import numpy as np
import corehydropy as ch25. Estimation methods
bulletin17c_analysis() and univariate_analysis() are convenience wrappers: pick a model, pick a method, get a frequency curve. Underneath them sits a smaller estimation layer that fits a model directly and hands back the fit itself: fit_mle(), fit_map(), fit_bayesian(), and fit_gmm(). This example fits the same model three different ways over the same censored record as example 23, and compares what each estimator returns.
What you’ll learn
- Fit one model three ways, by maximum likelihood, maximum a posteriori, and a Bayesian MCMC chain, and see how an informative prior moves the maximum a posteriori estimate.
- Compare AIC, BIC, and DIC across estimators, and why a method-of-moments fit like
fit_gmm()reports none of the three. - Set profile-likelihood confidence intervals against posterior credible intervals.
- Read a Bayesian fit’s R-hat and effective sample size.
Setup
The censored record
The same gauge record and historical information as example 23: forty-eight years of annual peaks, two pre-gauge floods known only within a range, and a fifty-year perception threshold that makes those two floods usable at all. See that example for the full derivation; here the frame is assembled directly.
peak_flows = [
6290, 2700, 13100, 16900, 14600, 9600, 7740, 8490, 8130, 12000,
17200, 15000, 12400, 6960, 6500, 5840, 10400, 18800, 21400, 22600,
14200, 11000, 12800, 15700, 4740, 6950, 11800, 12100, 20600, 14600,
14600, 8900, 10600, 14200, 14100, 14100, 12500, 7530, 13400, 17600,
13400, 19200, 16900, 15500, 14500, 21900, 10400, 7460,
]
n = len(peak_flows)
historical = {
"index": [n, n + 1],
"lower": [30000, 26000],
"value": [34000, 29000],
"upper": [38000, 32000],
}
perception = {
"start_index": n, "end_index": n + 49, "value": 25000, "number_above": 2,
}
flood_data = ch.analysis_data(exact=peak_flows, interval=historical, threshold=perception)An informative prior on skew
model_univariate() defaults to flat priors within each parameter’s bounds, and the mode of a flat-times-likelihood posterior is just the likelihood’s own mode. A fit_map() fit with no explicit prior would then land exactly on fit_mle()’s estimate, leaving nothing to compare. Log-Pearson III’s skew is exactly where a Bulletin 17C analysis leans on outside information (a regional skew estimate, in practice), so this model gives it a weakly informative prior instead of a flat one.
skew_prior = ch.Distribution("Normal", [-0.3, 0.5])
m = ch.model_univariate(
"LogPearsonTypeIII", flood_data,
parameters=ch.model_parameter("skew", prior=skew_prior),
)Three fits of the same model
f_mle = ch.fit_mle(m)
f_map = ch.fit_map(m)
f_bayes = ch.fit_bayesian(
m,
sampler="DEMCzs", iterations=400, output_length=2000,
seed=20250811, thinning_interval=1,
)
print(f"{'':<10}{'mean (log)':>14}{'sd (log)':>12}{'skew':>10}")
for label, fit in (("MLE", f_mle), ("MAP", f_map), ("Bayesian", f_bayes)):
v = list(fit.parameters.values())
print(f"{label:<10}{v[0]:>14.4f}{v[1]:>12.4f}{v[2]:>10.4f}") mean (log) sd (log) skew
MLE 4.0843 0.1991 -0.4657
MAP 4.0852 0.1962 -0.4189
Bayesian 4.0770 0.2063 -0.4263Every Fit, regardless of method, carries a .converged flag and a .status string. Here both are structurally True/"Success": the estimator raises on a genuine failure to estimate rather than handing back a failed fit, so these two fields record that the run completed, not that it necessarily found a good optimum.
The skew estimate moves from -0.4657 at MLE toward the prior’s mean of -0.3 once a prior is in play: -0.4189 at MAP, and -0.4263 for the Bayesian fit’s own point estimate, which is the posterior mean reported below in .posterior_summary (bit-identical, since fit_bayesian()’s default point_estimator is "PosteriorMean"; pass point_estimator="PosteriorMode" to report the MAP-like point instead). Both parameter priors and censored data pull in the same direction here, pulling the log-Pearson III skew away from the strongly negative value example 23 found for the systematic-only record.
Goodness of fit: AIC, BIC, DIC
print(f"{'method':<10}{'aic':>12}{'bic':>12}{'dic':>12}")
for label, fit, dic in (("MLE", f_mle, None), ("MAP", f_map, None), ("Bayesian", f_bayes, f_bayes.dic)):
dic_str = f"{dic:.4f}" if dic is not None else "NA"
print(f"{label:<10}{fit.aic:>12.4f}{fit.bic:>12.4f}{dic_str:>12}")method aic bic dic
MLE 990.1401 995.9356 NA
MAP 992.3973 998.1928 NA
Bayesian 992.6791 998.4746 989.7603AIC and BIC are defined at a single point estimate, so all three fits report them; DIC is built from the posterior draws themselves and only fit_bayesian() has those, hence the two NAs. MAP’s AIC is worse than MLE’s because AIC is built from whatever value the estimator maximized: the data log-likelihood for MLE, but the log posterior, data log-likelihood plus the prior’s own log-density, for MAP. Moving the skew toward the prior’s mean costs some data log-likelihood, and the prior’s log-density there adds nothing AIC counts as a credit, so the value AIC sees comes in below the data optimum. DIC, which does account for the posterior’s effective complexity, is close to but not equal to the Bayesian fit’s own AIC/BIC.
Not shown in this table: a fit_gmm() fit’s .aic/.bic/.log_likelihood come back None too, for a different reason. GMM is method-of-moments; it never had a likelihood surface to compute them from in the first place, so they are not “missing” the way a Bayesian-only DIC on an MLE fit is. They are structurally absent. See the GMM section below.
Profile intervals versus credible intervals
Fit.confint() dispatches on the fit’s own method: MLE and MAP get profile-likelihood confidence intervals, and a Bayesian fit gets its posterior credible interval.
ci_mle = f_mle.confint()
ci_map = f_map.confint()
ci_bayes = f_bayes.confint()
skew = "Skew (of log) (γ)"
print(f"{'':<22}{'lower':>10}{'upper':>10}")
for label, ci in (
("MLE (profile)", ci_mle), ("MAP (profile)", ci_map), ("Bayesian (credible)", ci_bayes)
):
print(f"{label:<22}{ci['lower'][skew]:>10.4f}{ci['upper'][skew]:>10.4f}") lower upper
MLE (profile) -0.8216 0.0835
MAP (profile) -0.7725 0.0531
Bayesian (credible) -0.8679 0.0730The Bayesian credible interval turns out to be the widest of the three, about 4% wider than the MLE profile interval, and MAP’s is the narrowest, pulled in by the same prior that shifted its point estimate. All three straddle zero: even with the historical record and, for MAP, an informative prior, this fit cannot rule out a skew near zero. That ordering is not a sign the Bayesian fit is less certain than the others. A profile interval traces the log-likelihood (or log-posterior) outward from the best fit until it drops by a chi-squared threshold, so it need not be symmetric around the point estimate: for the MLE and MAP fits above, each interval’s distance below its point estimate differs from its distance above. A credible interval, by contrast, is read directly off the spread of the sampled posterior draws, and the two need not land on the same width. They also answer different questions to begin with: a probability statement about the parameter given the data and prior, versus a repeated-sampling statement about the estimator. Close numeric agreement between them, in either direction, was never guaranteed.
R-hat and effective sample size
for nm, r, e in zip(f_bayes.parameter_names, f_bayes.posterior_summary["rhat"], f_bayes.posterior_summary["ess"]):
print(f"{nm:<22}{r:>8.3f}{e:>10.3f}")Mean (of log) (µ) 1.011 126.176
Std Dev (of log) (σ) 1.127 142.537
Skew (of log) (γ) 1.083 120.766R-hat compares between-chain and within-chain variance across the four parallel chains fit_bayesian() ran; values near 1 indicate the chains agree on where the posterior mass is. Effective sample size (ESS) is smaller than the raw output_length=2000 retained draws because consecutive MCMC draws are autocorrelated. An ESS in the low hundreds means each parameter’s posterior is effectively summarized by that many independent draws, not 2000.
A different estimator entirely: GMM
fit_gmm() only fits a model_bulletin17c() model, since IGMMModel has exactly one implementation, so it is not a fourth entry in the table above. It is, however, the estimator bulletin17c_analysis() itself builds on.
f_gmm = ch.fit_gmm(ch.model_bulletin17c(flood_data))
for label, v in zip(("mean (log10)", "sd (log10)", "skew"), f_gmm.parameters.values()):
print(f"{label:<14}{v:.4f}")
print(f_gmm.summary())mean (log10) 4.0862
sd (log10) 0.2032
skew -0.6863
<Fit GMM (Success): 3 parameters>
parameters:
p1: 4.08623
p2: 0.203167
p3: -0.686259
j-statistic: not interpretable at 0 degrees of freedom gmm iterations: 2
converged: True
standard errors:
p1: 0.028255
p2: 0.0233435
p3: 0.486138The J-statistic is GMM’s own goodness-of-fit diagnostic, an overidentification test, and the method-of-moments analogue of a likelihood-based AIC/BIC/DIC rather than a value on the same scale as any of them. Bulletin 17C’s moment conditions exactly match its parameter count, so this fit has zero over-identifying degrees of freedom and .j_stat_pval is None. .j_stat itself carries no information there either, which is why summary() names the reason instead of showing a number: the residual covariance the statistic is scaled by is theoretically zero at zero degrees of freedom, so what comes back is whatever inverting a numerically singular matrix happened to give. The attribute is still on the fit for anyone who wants it.
Writing more moment conditions than parameters with fit_gmm_moments() restores the p-value. It does not make the statistic itself trustworthy: the same residual covariance has rank q - p for any q and p, so it is singular whether or not the fit is over-identified, and inverting it amplifies the optimizer’s convergence tolerance. docs/upstream-csharp-issues.md carries the measurement.
quantile_variance() gives the delta-method variance of a fitted quantile off a GMM fit’s sandwich covariance, useful directly without going through bulletin17c_analysis()’s full uncertainty-quantification machinery:
ch.quantile_variance(f_gmm, 0.01)0.00432876400509739Diagnostics off the MAP fit
fit_diagnostics() computes Cook’s distance and leverage at a point-estimate fit’s optimum (MAP or GMM), or PSIS-LOO Pareto-k and prior influence off a Bayesian fit’s posterior.
diag = ch.fit_diagnostics(f_map)
for i in range(5):
print(f"{i + 1:>3} cooks_distance={diag['cooks_distance'][i]:.4f} leverage={diag['leverage'][i]:.4f}")
idx = int(np.argmax(diag["cooks_distance"]))
print(f"largest Cook's distance: observation {idx + 1}, a {peak_flows[idx]:,} cfs peak") 1 cooks_distance=0.0169 leverage=0.0331
2 cooks_distance=0.2842 leverage=0.4024
3 cooks_distance=0.0028 leverage=0.0066
4 cooks_distance=0.0060 leverage=0.0142
5 cooks_distance=0.0036 leverage=0.0084
largest Cook's distance: observation 2, a 2,700 cfs peakThe single largest Cook’s distance in the record belongs to the second systematic year, 2,700 cfs, an unusually small peak that pulls a log-space fit the way example 23’s Multiple Grubbs-Beck low-outlier screen exists to guard against.
Reproduction check
# The pipeline is deterministic (MLE/MAP by a deterministic optimizer; the Bayesian fit is
# seeded), so these literals are what this port computes -- and the R version asserts exactly
# the same ones, which is what proves the cross-language identity. Both packages run the same
# compiled core with a bit-exact Mersenne Twister, so the match is exact here rather than merely
# close.
mean_key, sd_key, skew_key = f_mle.parameter_names
assert f_mle.parameters[mean_key] == 4.084346130939519
assert f_mle.parameters[sd_key] == 0.19914759083044403
assert f_mle.parameters[skew_key] == -0.4657240567252959
assert f_mle.aic == 990.140077715221
assert f_mle.bic == 995.935554613394
assert f_map.parameters[mean_key] == 4.085236124852605
assert f_map.parameters[sd_key] == 0.19616835592097398
assert f_map.parameters[skew_key] == -0.4189188626281035
assert f_map.aic == 992.3973115540235
assert f_map.bic == 998.1927884521965
assert f_bayes.parameters[mean_key] == 4.076977143876432
assert f_bayes.parameters[sd_key] == 0.2063000968835454
assert f_bayes.parameters[skew_key] == -0.4263428595238764
assert f_bayes.aic == 992.679135891769
assert f_bayes.bic == 998.474612789942
assert f_bayes.dic == 989.7603249256322
# Profile versus credible intervals on skew.
assert ci_mle["lower"][skew] == -0.821569546770356
assert ci_mle["upper"][skew] == 0.08353932473268728
assert ci_map["lower"][skew] == -0.7725382794403562
assert ci_map["upper"][skew] == 0.05314301785784643
assert ci_bayes["lower"][skew] == -0.8679286484920807
assert ci_bayes["upper"][skew] == 0.0730369777667995
# R-hat and ESS.
assert f_bayes.posterior_summary["rhat"].tolist() == [
1.0112377055613953, 1.1271063062090425, 1.0830639682776986
]
assert f_bayes.posterior_summary["ess"].tolist() == [
126.17628141940655, 142.53659604889492, 120.76568406444153
]
# GMM: the same censored record, a different estimator.
assert list(f_gmm.parameters.values()) == [
4.086226682554578, 0.2031671188484615, -0.686258748960609
]
# .j_stat is deliberately NOT pinned: at zero degrees of freedom it is the inverse of a
# theoretically zero matrix, and it is not reproducible across optimizers or across languages.
assert f_gmm.j_stat_pval is None
assert ch.quantile_variance(f_gmm, 0.01) == 0.00432876400509739
# Diagnostics off the MAP fit.
assert diag["cooks_distance"][1] == 0.2841569426904578
assert diag["leverage"][1] == 0.40243088779094255
# Internal consistency: the informative prior moves MAP away from MLE and toward the prior
# mean, every confidence/credible interval brackets its own point estimate, and GMM's
# log-likelihood-based fields are structurally None.
assert abs(f_map.parameters[skew_key] - (-0.3)) < abs(f_mle.parameters[skew_key] - (-0.3))
assert ci_mle["lower"][skew] < f_mle.parameters[skew_key] < ci_mle["upper"][skew]
assert ci_map["lower"][skew] < f_map.parameters[skew_key] < ci_map["upper"][skew]
assert f_gmm.aic is None and f_gmm.bic is None and f_gmm.log_likelihood is None
print("All reproduction checks passed.")All reproduction checks passed.