CHESS JPM

CHESS: JPM (Joint Probability Method) Engine

Source files: CHESS_Library/Extreme_Value_Analysis/CHESS_JPM.m (entry point), CHESS_Library/Extreme_Value_Analysis/JPM/CHESS_JPM_integration.m (the integral itself), CHESS_Library/Extreme_Value_Analysis/JPM/CHESS_JPM_plot.m (plotting only). Called from: CHESS_Library/Extreme_Value_Analysis/call_hazard_curve_builder.m, itself invoked from the TC branch of chess_pros.m — see CHESS PROS for the workflow context (config.workflow ∈ {1,4}, config.storm_sampling selecting 'TC'/'CC'). Sibling engine: CHESS PST — the Stochastic Simulation Technique (SST) that plays the same hazard-curve-building role for extratropical (XC) storms. This page assumes familiarity with the "two storm populations, two statistical engines" framing already introduced on CHESS PROS.

What JPM is

Tropical cyclones are comparatively rare at any single site: a hurricane making a damaging, close pass to a given project location may only have happened a handful of times — or never — in the actual historical/hindcast record. That sparsity is exactly what makes SST's resampling approach (build a longer synthetic record by bootstrapping the historical one) impractical for TC. There simply isn't enough historical signal at a point to resample from.

JPM solves this differently. Instead of resampling a short observed record, it works from a large synthetic storm suite that spans the plausible TC parameter space directly — track, central pressure deficit, radius of maximum winds, translation speed, heading, and so on — as already constructed upstream by CHS's storm climatology and hydrodynamic modeling. Each synthetic storm in that suite is assigned a discrete probability-of-occurrence weight representing its climatological annual rate, rather than being treated as one draw from an observed sequence. The hazard curve is then built by integrating response magnitude against those weights directly, rather than by counting how often a given magnitude was empirically observed.

Theory: the JPM integral

The published form of the integral JPM implements (Melby et al. 2021, Appendix A/G — see References) is:

P(R>r) = ∫∫ P[R(θ)+ε>r | i,θ] f(θ) p(θ) dθ

where θ is the storm parameter vector (distance to the reference location, central pressure deficit, radius of maximum winds, translation speed, heading direction, and secondary parameters such as Holland B and tide), f(θ) is the joint probability density of those parameters, p(θ) is the annual rate of storm occurrence, and ε is a residual/uncertainty term. Conceptually, f(θ)p(θ)dθ says: take a slice of the storm-parameter space, weight it by how likely a storm with those characteristics is to occur and by how often storms occur at all, and integrate the response over every such slice.

CHESS implements this as a discrete integral rather than a continuous one: the synthetic storm suite already samples the parameter space at a finite set of points, so each synthetic storm stands in for one slice dθ of that space, carrying its own precomputed weight instead of a density function to be integrated analytically. "Synthetic storm suite weighted by climatological annual rate" means exactly this — each storm is not a historical event being counted, but a representative point in parameter space whose associated weight already encodes how often storms with roughly those characteristics occur per year.

Algorithm walkthrough

The discrete summation

CHESS_JPM_integration.m receives a response population and a matching per-storm (or per-storm-replicate) probability weight, and builds the exceedance curve as a direct discrete summation — no bootstrap, no Monte Carlo resampling at this stage:

pass_indx = ~isnan(Resp) & Resp>0;
Resp = Resp(pass_indx,1); ProbMass = ProbMass(pass_indx,1);
[y, I] = sort(Resp, 'descend');            % sort response high to low
x = cumsum(ProbMass(I,1));                 % x = P(R>r): cumulative discrete storm weight

Sorting every storm/replicate's response from largest to smallest and cumulatively summing its probability weight directly produces P(R>r) at every observed response level — the discretized analog of ∫∫f(θ)p(θ)dθ. The result is then log-transformed, de-duplicated, and interpolated onto standard AEF/AEP grids (HC_plt_x, HC_tbl_x) to produce the plotting curve and summary hazard table.

Where the weights come from. The probability weight consumed here — ProbMass = response_data.data(:,4), documented in CHESS_JPM.m as "TC event Discrete Storm Weights (DSW)" — is not an empirical count of how often a storm was observed. It traces back to chs_probability_mass_loader.m, which loads a regional CHS *_ITCS_DSW_600km.mat file's DSW_ITCS field (or Freq(Nsvpt).TC for the NACCS region): a genuine climatological annual-frequency table (storms/year) computed upstream as part of CHS's PCHA framework. That value flows through chess_pros_formater.m as project_forcing.TC.TC_Prob, into call_hazard_curve_builder.m (split evenly across uncertainty-replicate columns via tc_dsw = repmat(TC_Prob./dcols, 1, dcols) when the response has multiple replicates), and finally into column 4 of response_data.data.

This is the key structural asymmetry with the sibling PST/SST engine, worth understanding explicitly. Because CHS's DSW/TC_Freq values are already annual rates, JPM never needs a record-length ("Nyrs") normalization step anywhere in its code. PST, by contrast, works from a finite historical-plus-synthetic record and has to convert an empirical or bootstrapped event count into a rate using config.Nyrs_XC (see CHESS PST) — a real, confirmed methodological difference between the two engines, not an implementation detail that happens to differ.

The two-layer uncertainty mechanism

JPM uses a 444-point standard-normal discretization table (discrete_norm_444.txt, spanning roughly −3.00 to +3.00) as a deterministic stand-in for Monte Carlo sampling of assigned uncertainty. Critically, which of two structurally different places that discretization is applied depends on which response variable is being processed — this is call-site-dependent, not a single fixed mechanism, and is the most important subtlety in how JPM actually works:

1. Derived structure-response variables (R2%, Dn50, q, P1, etc.). For these, the 444-point discretization happens upstream of JPM entirely, inside the forcing itself. chess_pros_formater.m replicates each storm's SWL/Hm0/Tp into 444 columns, and the shared uncertainty engine (call_uncertainty_engine.m / apply_uncertainty.m / pcha_forcing_uncertainty.m) applies the same 444-length z-score vector identically across every storm — an epistemic/system-wide uncertainty model, not per-storm aleatory noise — perturbing SWL/Hm0 before the structure-response equations ever run. By the time this reaches JPM, call_hazard_curve_builder.m sets eva_options.integration_method = 1 ("PCHA ATCS"), and CHESS_JPM_integration.m's own internal discretization block (integration_method == 2) is skipped entirely. JPM simply treats the already-444-perturbed, flattened response population as final discrete atoms and runs the plain discrete summation described above.

2. The three raw forcing variables (SWL, Hm0, Tp) themselves. Their own hazard curves bypass the upstream uncertainty engine — they're built from chess_pros.m's un-replicated, per-storm-max _no_rep fields — so call_hazard_curve_builder.m instead forces integration_method = 2 ("Let JPM Manage Replicates," a literal comment in the code), and CHESS_JPM_integration.m performs the discretization itself: it reads discrete_norm_444.txt, replicates it 444× per storm, and applies a variance-partition trick — a fixed 0.1-unit (absolute) / 0.1-fraction (relative) "first partition" (p1_a/p1_r) is carved out of the total assigned uncertainty (chs_swl_u_a/u_r, chs_hm0_u_a/u_r) and applied directly to perturb the response across the 444 replicate states, shaping the central/mean curve. The residual uncertainty left after subtracting that 0.1 partition is what feeds the analytic confidence-limit formula below instead of the discrete perturbation.

The practical takeaway: the answer to "what does the 444-point discretization actually perturb?" is call-site-dependent. For derived responses it perturbs the forcing (SWL/Hm0) upstream of JPM; for the three raw forcing variables' own hazard curves it perturbs the response value directly inside JPM via the partition trick. Both paths draw from the same 444-value table but combine it with the uncertainty differently.

Confidence limits: analytic, not bootstrap

Once the discrete hazard curve y (response) vs. x (cumulative exceedance) is built, JPM computes confidence-limit bands with a closed-form Gaussian-percentile formula applied directly to that curve — CHESS_JPM_integration.m, selected by jpm_options.uncertainty_treatment:

• 'absolute': y + z·√(U_a² + U_t²)
• 'relative': y + z·y·√((y·U_r)² + U_t²)
• 'combined': y + z/√(1/U_a² + 1/(y·U_r)² [+ 1/U_t²])

where z = norminv(prc/100) and U_t is tide uncertainty (see "Tide handling" below).

This is not resampling. Contrast this explicitly with PST/SST, which derives its confidence limits from genuine bootstrap resampling (on the order of 100 replicates, empirical-record resampling with jitter — see CHESS PST). JPM's percentile bands come from a single analytic formula evaluated at each point on the already-computed best-estimate curve, with no repeated resampling involved.

Because both engines are built off the same shared options-struct constructor in call_hazard_curve_builder.m, JPM is passed — and silently ignores — several fields that are meaningful only to PST: bootstrap_sims, GPD_TH_crit, apply_GPD_to_SS, Nyrs, DataType. This is harmless (each engine only reads the fields relevant to its own math) but worth knowing: one shared struct constructor serves two very differently-implemented statistical engines.

The "50" label caveat. call_hazard_curve_builder.m labels JPM's output columns Output.CL = [50, prc]. The "50" here is just a label attached to JPM's raw, discretely-integrated central curve — it is not a literal 50th percentile computed via the analytic CL formula above. Only the columns for the percentiles in prc are true CL_unc-derived values. (The PST page covers this same "50" labeling nuance in more depth, since it applies identically there.)

Tide handling (present but dormant)

CHESS_JPM_integration.m implements a complete tide_application ∈ {0,1,2} branch: 0 = no tide, 1 = tide folded into the confidence-limit formula as an added variance term, 2 = tide added directly to the response (as a random draw, a discrete-Gaussian draw, or literal skew-tide values, depending on ind_Skew). In practice, however, call_hazard_curve_builder.m always hardcodes tide_application: 0 — this branch is real, complete, tested-looking code, but it is never exercised by the current, sole caller. Random-tide sampling for PROS instead happens earlier in the pipeline, via config.apply_random_tides in call_project_forcing_adjuster.m, entirely outside JPM.

Config fields

JPM itself takes no config directly — call_hazard_curve_builder.m translates a handful of config fields into the jpm_options/response_data structs described above. Full field-by-field definitions live in config_reference; this table only maps which fields feed JPM and how.

Field	Effect
config.storm_sampling	Gates whether the 'TC'/'CC' branch (and therefore JPM) runs at all
config.pros_use_aep	AEP vs. AEF convention for the hazard curve
config.project_CLs	Percentile list used by the analytic CL formula (max 4, plus the "50" central curve)
config.chs_swl_u_a/u_r, config.chs_hm0_u_a/u_r	Feed U_a/U_r for the SWL/Hm0/Tp forcing hazard curves, and (via the response library headers) several derived responses
config.r2p_u, config.q_u, config.dn50_u, config.p1_u	Feed U_a/U_r for the respective structure responses

Not configurable — hardcoded in source: the 444-point discretization table itself (discrete_norm_444.txt), the 0.1 variance-partition constants (p1_a, p1_r), and the integration_method selection logic (JPM's caller decides 1 vs. 2 based on which variable is being processed, not a user-facing switch).

Implementation reference

• Signature: [JPM_output] = CHESS_JPM(response_data, jpm_options, plot_options).
• response_data.data must be an N×4 matrix: [time, response, skew_tides, dsw]. Column 1 ("time") is always zeroed by the caller and unused by the integration itself.
• jpm_options fields: integration_method (1 = "PCHA ATCS" — Probabilistic Coastal Hazard Analysis with Gaussian-process-metamodel-Augmented Tropical Cyclone Suite; 2 = "PCHA ITCS" — standard suite used as-is), uncertainty_treatment ('absolute'|'relative'|'combined'), tide_application (0/1/2, dormant per above), ind_Skew, use_AEP, prc (max 4 percentiles; defaults to [2.28 15.87 84.13 97.72] — the ±1σ/±2σ Gaussian quantiles — if left empty), stat_print.
• Output: mapped by call_hazard_curve_builder.m onto the same Output.x_plot/y_plot/x_table/y_table/x_table_ARI/tbl_rsp_x/tbl_rsp_y/CL struct contract used by PST — sourced from JPM's own HC_plt_x/HC_plt/HC_tbl_x/HC_tbl/HC_tbl_rsp_x/HC_tbl_rsp_y fields, so downstream consumers (documented on CHESS PROS) don't need to distinguish which engine produced a given result.
• Entry chain: chess_pros.m (TC branch) → call_hazard_curve_builder.m → CHESS_JPM.m → CHESS_JPM_integration.m.
• Two dormant/dead-branch issues, informational only: inside CHESS_JPM.m's create_plots==1 branch (never reached, since the sole caller hardcodes create_plots=0), there is a reference to an undefined path_out variable and a save-path concatenation missing a separator. Both are harmless today but would break immediately if a future caller ever enabled in-function plotting.

Naming and terminology notes

• Stale header comment. CHESS_JPM.m's first line reads %% CHESS_JPM_Tool.m — a name that doesn't match the actual filename CHESS_JPM.m.
• Copy-paste docstring artifact. The same file's docstring opens with "SST Input Data Structures:" — clearly copied from the sibling PST/SST tool's docstring and never updated for the JPM copy. This is a documentation artifact, not a sign of shared code.
• Alpha-version banner. The runtime banner reads *** CHESS-JPM Tool Alpha Version 0.3 *** Release 1 - 20210831 *** FOR TESTING ONLY *** — still labeled alpha/testing-only despite being live in the production PROS pipeline. The banner never actually prints in production (it's gated behind stat_print=1, which the caller always hardcodes to 0), but the alpha-version status is worth surfacing regardless.
• integration_method naming is real, deliberate terminology — not an inconsistency. Unlike the items above, "PCHA ATCS" and "PCHA ITCS" are genuine, documented methodology names given verbatim in CHESS_JPM.m's own header, each citing real prior USACE studies as provenance (ATCS: South Atlantic Coast Study, Louisiana Coast Protection and Restoration Study, Coastal Texas Study revision; ITCS: North Atlantic Coast Comprehensive Study, Coastal Texas Study initial study). Treat this as a legitimate, deliberate toggle between two established methodologies, not documentation drift.
• An orphaned, unrelated second implementation. CHESS_Library/CHS/CHS_Data_Prep/itcs_hazard_remap.m implements a structurally similar discrete-integral pattern (sort/cumsum/interp1) internally labeled "JPM INTEGRATION," using entirely different variable names and terminology (DSW/ATCS/ITCS), with zero in-repo callers and no call relationship to CHESS_JPM.m. Don't confuse this dead script with the active JPM pipeline documented on this page.
• Open question: combining JPM and PST curves by percentile label. When PROS combines the XC (PST/bootstrap-derived) and TC (JPM/analytic-Gaussian-derived) hazard curves into the "CC" combined curve, call_hazard_curve_combiner.m treats both engines' percentile-labeled output columns (CL = [50, prc]) as directly comparable/combinable by label alone — even though the two curves were derived by genuinely different statistical mechanisms (bootstrap resampling vs. a closed-form Gaussian formula). Whether this is a real accuracy concern or an accepted simplification within the CHS/PCHA methodology could not be determined from the code alone; a full trace of combine_hazard_curves.m's internals was outside this investigation's scope. This is flagged here as an open question, not an asserted bug.

References

• Nadal-Caraballo, N.C., Campbell, M.O., Gonzalez, V.M., Torres, M.J., Melby, J.A., and Taflanidis, A.A. (2020). "Coastal Hazards System: A Probabilistic Coastal Hazard Analysis Framework." Journal of Coastal Research, Vol. 95, No. sp1. DOI: 10.2112/SI95-235.1.
• Melby, J.A., Massey, T.C., Stehno, A.L., Nadal-Caraballo, N.C., Misra, S., and Gonzalez, V.M. (2021). "Sabine Pass to Galveston Bay, TX Pre-construction, Engineering and Design (PED): Coastal Storm Surge and Wave Hazard Assessment. Report 1 – Background and Approach." ERDC/CHL TR-21-15, Report 1. US Army Engineer Research and Development Center. Appendix G, "Stochastic Structure Response Simulations," pp. 190-206. DOI: 10.21079/11681/41820.