REVIEW 3 major objections 4 minor 12 references
Joint modeling of low and high extremes using a multivariate extended generalized Pareto distribution
T0 review · 3 major / 4 minor · reviewed 2026-08-05 · deepseek-v4-flash
Pith's one-line read The paper proposes a multivariate extended generalized Pareto distribution that models the entire range of a random vector—low, moderate, and high—while keeping joint lower and upper tails compliant with extreme-value theory and governed by
desk verdict New threshold-free multivariate eGPD construction is a real contribution, but the application's neural credible intervals run outside the training prior and need fixing. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The object that carries the argument is the stochastic construction Y = R([1 − ω{F_R(R)}]L + ω{F_R(R)}U). R is a univariate eGPD with lower-tail shape κ and upper-tail shape ξ; L and U are positive random vectors on the L1 simplex, independent of R; ω is a cumulative distribution function on [0,1] acting as a smooth switch between L and U as the radius R moves from 0 to ∞. The identity ||Y||_1 = R makes the radius exactly the univariate eGPD, and the tail limits ω→1 and ω→0 transfer the angular dependence to U and L respectively. This is what breaks the dependence structure into two separately controllable pieces while preserving a smooth bulk.
What would settle it
Simulate from the model with known parameters and a weight function that approaches 1 slowly, such as ω(u) = 1 − (1 − u)^β with small β, then compute the multivariate angular measure or extremal coefficient at very high radial thresholds; if the observed tail dependence differs systematically from U's angular distribution, the claimed reduction of the upper tail to RU has failed. Analytically, check whether E[{R(1 − ω(F_R(R)))}^α L_j^α] is negligible relative to E[{R ω(F_R(R))}^α U_j^α] for α = 1/ξ; if the ratio does not vanish as the tail threshold increases, the separation does not hold.
Extended reading notes
Core claim
The central construction is Y = R([1 − ω{F_R(R)}]L + ω{F_R(R)}U), with R a univariate eGPD, L and U independent simplex-valued random vectors, and ω a cumulative distribution function on [0,1]. Because ||Y||_1 = R, the radius is exactly the univariate eGPD. As the radius tends to infinity, ω approaches 1 and Y behaves like RU; as it tends to zero, ω approaches 0 and Y behaves like RL. Using Breiman's lemma and its multivariate extension, the paper shows the upper joint tail is regularly varying with index 1/ξ and limiting angular measure determined by U, while the lower joint tail has index κ (or, in a special case, the index of L^{-1}) with angular measure determined by L. Margins are heavy
Load-bearing premise
The paper asserts that in the far upper tail the term [1 − ω{F_R(R)}]L becomes negligible compared with ω{F_R(R)}U, so Y behaves like RU, and symmetrically at low values Y behaves like RL; if the weight function does not reach its limits fast enough relative to R's tail, that separation of tail dependence into U and L may fail.
Editorial extensions
If this is right
- A single fitted model replaces the usual three-piece procedure (lower-tail GP, bulk, upper-tail GP) and avoids the threshold selection step in multivariate settings, where such choices are especially fragile.
- Joint exceedance probabilities in both tails are EVT-compliant: high-level joint probabilities follow a regular variation limit whose angular measure is fixed by U, and low-level joint probabilities follow the analogous limit set by L.
- Because both limiting angular measures are controlled by free parameters in L and U, practitioners can model drought-type joint lows and flood-type joint highs with different strengths of dependence in one model.
- Fast simulation makes simulation-based inference practical; the paper's neural estimators provide point estimates or posteriors orders of magnitude faster than the hybrid likelihood-moment benchmark.
- For rainfall, fits from the model reproduce observed χ-measures and quantile-quantile behavior from the lowest to highest values, matching classical bivariate GP fits in the tails while also covering the bulk.
Reading between the lines
- Because L and U can be any positive simplex-valued vectors, the construction offers a general template for full-range dependence modeling; a natural next step is to choose L or U so that one or both tails exhibit asymptotic independence, a case the paper explicitly leaves open.
- The radial-angular separation is only asymptotic: for moderate values of R, both L and U contribute to dependence. If a user needs exact control at intermediate quantiles, the weight function ω would need to be chosen to shape that transition, and inference on θ_ω may be the weakest link (the paper's simulations show lower posterior contraction for θ_ω).
- The paper notes that L and U contaminate marginal scales; when sites are not exchangeable, the component-specific scale extension it sketches would be needed. A testable extension is to fit that extended model and check whether the estimated L/U angular measures shift materially.
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. The paper proposes a new multivariate extended generalized Pareto distribution (eGPD) for jointly modeling low, moderate, and high values, with the stated aim of avoiding threshold selection while remaining compliant with multivariate extreme-value theory in both tails. The model is defined through the stochastic representation Y = R([1-ω{F_R(R)}]L + ω{F_R(R)}U), where R is a univariate eGPD radial variable, L and U are independent simplex-valued angular vectors, and ω is a distribution function on [0,1] that transitions from L to U as R moves from small to large values. The authors claim that the upper joint tail is governed by U, the lower joint tail by L, and that the margins behave like univariate eGPDs. Inference is performed using amortized neural Bayes estimators and neural posterior estimators, after also developing a hybrid likelihood--moment estimator. The methodology is illustrated on bivariate Dutch daily rainfall data from three stations. The central contribution is a generative full-range multivariate model with separately controlled lower- and upper-tail dependence.
Significance. If the mathematical claims hold, this is a useful contribution to the growing literature on bulk--tail models: it provides a simple, fast-to-simulate multivariate construction that avoids threshold selection and separates lower- and upper-tail dependence parameters. The authors are to be credited for making the code publicly available, for using amortized simulation-based inference, and for providing a fairly extensive simulation study with calibration diagnostics. The paper also honestly lists several limitations of the modeling framework, including the fact that L and U affect margins and that asymptotic independence is not covered. The main value is conceptual: a clean stochastic representation whose upper and lower extremal dependence can be tuned independently. However, as detailed below, the current manuscript contains a definitional ambiguity in the weight function, a missing proof of the central tail-equivalence assertion, and, most importantly, an application whose reported neural posterior estimates lie outside the training prior support and therefore cannot be taken at face value.
major comments (3)
- [Table 2 and §4.1] The NPE was trained with θ_U ∼ Unif(0.1,20) and θ_ω ∈ (0,0.5), but Table 2 reports for Ammerzoden–Zaltbommel a posterior median θ_U = 21.043 with 95% CI (19.358, 22.980), and for Giersbergen–Zaltbommel a θ_ω CI (0.417, 0.504) whose upper endpoint exceeds the parameter-space bound 0.5. Under the stated prior, the exact posterior has zero mass outside (0.1,20) and (0,0.5), respectively. These results therefore indicate that the normalizing flow is extrapolating outside its training support, so the reported posterior intervals and the conclusions drawn from them are not currently valid. The authors should either widen the priors (and retrain) or explicitly constrain the flow to respect the parameter bounds, and then reassess the application.
- [§4.1] The weight function is not well-defined as written. The text defines ω(u) as a Beta(3,3) CDF evaluated at (u−θ_ω)/(1−2θ_ω) with θ_ω ∈ (0,0.5). For u∈[0,1], this argument takes values below 0 (e.g., at u=0, it equals −θ_ω/(1−2θ_ω)) and above 1 (e.g., at u=1, it equals (1−θ_ω)/(1−2θ_ω) > 1). Since the Beta CDF is only defined on [0,1], the model (5) is not fully specified. The authors need to clarify whether the argument is clipped, renormalized, or replaced by a different transition function, and verify that the claimed limits ω(u)→0 as u→0 and ω(u)→1 as u→1 still hold.
- [§2.3.2] The central tail-separation result is asserted rather than proved. The text states that the upper tail of Y can be studied 'equivalently' through RU 'since ω(u)↑1 as u↑1', and similarly for the lower tail. This is a load-bearing step: it is needed to conclude that the joint upper (lower) tail of Y is fully governed by U (L). The convergence ω{F_R(R)}→1 in probability as R→∞ is plausible and likely sufficient, but the paper should supply a precise argument, e.g., a lemma showing that the remainder R(1−ω{F_R(R)})L (and the analogous lower-tail term) is asymptotically negligible in the relevant tail probability. Without such a proof, the main mathematical claim rests on an unverified heuristic.
minor comments (4)
- [§4.1] The variance-stabilizing transformation is written as 'h(x) = sign(x) log(1+|x|)−1', which is ambiguous. Please clarify the intended formula.
- [§2.3.2] Equation (6) and (7): it would be helpful to state explicitly that the moment formulas are only valid when the relevant moments of R exist, and to define ω_R before using it in (6).
- [§4.2] The SBC and z-score diagnostics are computed on test parameters drawn from the training prior. This is appropriate for checking calibration within the prior, but it does not validate behavior outside the prior support; a brief sentence acknowledging this limitation would be useful.
- [§5] The plots in Figure 5 and Appendix A are informative, but the lower-tail χ-measure appears to show systematic discrepancy at small levels; the caption mentions discretization of rainfall, but the paper would benefit from a more quantitative discussion of the magnitude of this lack-of-fit.
Circularity Check
No significant circularity: the model's tail properties follow from its explicit stochastic construction and are validated on independent simulations; the NPE prior-support issue is a validity concern, not circularity.
full rationale
The derivation chain is self-contained. The tail properties in Section 2.3.2 are consequences of the explicit stochastic representation (5), the exact identity ||Y||_1 = R in (8), and the assumed limits of the weight function ω; no fitted parameter is renamed as a prediction. The upper-tail equivalence with RU is asserted without a full proof, but it is a standard regular-variation/Slutsky consequence and does not assume the conclusion, so it is a mathematical gap rather than circularity. The neural estimators are trained on simulated data and validated on a separate test set (Section 4.2) with simulation-based calibration diagnostics; the application uses the trained NPE to produce posterior estimates, not predictions forced by a fitted input. The only caveat is non-circular: Table 2 reports credible intervals for θ_U and θ_ω that partially lie outside the training prior support (θ_U ~ Unif(0.1,20), θ_ω ~ Unif(0,0.5)), which undermines the validity of those NPE estimates but does not constitute circularity. Minor self-citations (Naveau et al. 2016; Sainsbury-Dale et al. 2024) provide background methods and are not load-bearing for the new model's mathematical construction.
Assumptions & free parameters
free parameters (6)
- kappa (lower tail shape of R) =
1.116 (Ammerzoden-Zaltbommel)
- sigma (scale of R) =
1.381 (Ammerzoden-Zaltbommel)
- xi (upper tail shape of R) =
0.193 (Ammerzoden-Zaltbommel)
- theta_L (lower dependence parameter) =
4.082 (Ammerzoden-Zaltbommel)
- theta_U (upper dependence parameter) =
21.043 (Ammerzoden-Zaltbommel)
- theta_omega (bulk transition parameter) =
0.237 (Ammerzoden-Zaltbommel)
assumptions (4)
- domain assumption R, L, and U are mutually independent in (5).
- domain assumption Asymptotic dependence in both upper and lower tails.
- ad hoc to paper Tail equivalence: Y behaves as RU in the upper tail and as RL in the lower tail.
- standard math Breiman's lemma and the multivariate extension of Fougeres and Mercadier (2012).
invented entities (1)
-
Latent angular vectors L and U (on the L1 simplex)
Cite this review
Pith. "Pith review of Joint modeling of low and high extremes using a multivariate extended generalized Pareto distribution." pith.science (2026). https://pith.science/paper/VSDCWEEP
@misc{pith2026250905982,
author = {Pith},
title = {Pith review of: Joint modeling of low and high extremes using a multivariate extended generalized Pareto distribution},
year = {2026},
howpublished = {\url{https://pith.science/paper/VSDCWEEP}},
note = {Machine review of arXiv:2509.05982}
}
read the original abstract
In most risk assessment studies, it is important to accurately capture the entire distribution of the multivariate random vector of interest from low to high values. For example, in climate sciences, low precipitation events may lead to droughts, while heavy rainfall may generate large floods, and both of these extreme scenarios can have major impacts on the safety of people and infrastructure, as well as agricultural or other economic sectors. In the univariate case, the extended generalized Pareto distribution (eGPD) was specifically developed to accurately model low, moderate, and high precipitation intensities, while bypassing the threshold selection procedure usually conducted in extreme-value analyses. In this work, we extend this approach to the multivariate case. The proposed multivariate eGPD has the following appealing properties: (1) its marginal distributions behave like univariate eGPDs; (2) its lower and upper joint tails comply with multivariate extreme-value theory, with key parameters separately controlling dependence in each joint tail; and (3) the model allows for fast simulation and is thus amenable to simulation-based inference. We propose estimating model parameters by leveraging modern neural approaches, where a neural network, once trained, can provide point estimates, credible intervals, or full posterior approximations in a fraction of a second. Our new methodology is illustrated by application to daily rainfall times series data from the Netherlands. The proposed model is shown to provide satisfactory marginal and dependence fits from low to high quantiles.
Figures
Figures from the paper (3 more)
Reference graph
Works this paper leans on
-
[2]
ArXiv preprint arXiv:2503.06255
Gaussian mixture copulas for flexible dependence modelling in the body and tails of joint distributions. ArXiv preprint arXiv:2503.06255. André, L., Wadsworth, J., O’Hagan, A.,
-
[4]
Implications of Artificial Intelligence on Health Data Privacy and Confidentiality
Fast likelihood-free parameter estimation for Lévy processes. ArXiv preprint arXiv.2501.01639. Davison, A.C., Huser, R.,
-
[6]
ArXiv preprint arXiv.2404.10899
A variational neural Bayes framework for inference on intractable posterior distributions. ArXiv preprint arXiv.2404.10899. Mackay, E., Jonathan, P.,
-
[7]
ArXiv preprint arXiv:2310.12711
Modelling multivariate extremes through angular-radial decomposition of the density function. ArXiv preprint arXiv:2310.12711. McFadden, D.,
-
[8]
Theoretical guarantees for neural estimators in parametric statistics. ArXiv:2506.18508. Sainsbury-Dale, M.,
-
[9]
R package version 0.1.2,https://CRAN.R-project.org/package= NeuralEstimators
NeuralEstimators: Likelihood-Free Parameter Estimation us- ing Neural Networks. R package version 0.1.2,https://CRAN.R-project.org/package= NeuralEstimators. Sainsbury-Dale, M., Richards, J., Zammit-Mangion, A., Huser, R., 2025a. Neural Bayes es- timators for irregular spatial data using graph neural networks. Journal of Computational and Graphical Statis...
-
[10]
ArXiv preprint arXiv:2412.07957
Spatial scale-aware tail dependence modeling for high-dimensional spatial extremes. ArXiv preprint arXiv:2412.07957. Smith, R.L.,
-
[2018]
ArXiv preprint arXiv:1804.06788
Validating Bayesian inference algorithms with simulation-based calibration. ArXiv preprint arXiv:1804.06788. Tencaliec, P., Favre, A.C., Naveau, P., Prieur, C., Nicolet, G.,
Show all 12 references
-
[2019]
secondary
Natural catastrophes and man-made disasters in 2018: “secondary” perils on the frontline. (Tech. Rep. No. 2). Säilynoja, T., Bürkner, P., Vehtari, A.,
2018
-
[2023]
(Eds.), Contribution of Working Groups I, II, and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change
Climate Change 2023: Synthesis Report, in: Core Writing Team, H.L., Romero, J. (Eds.), Contribution of Working Groups I, II, and III to the Sixth Assessment Report of the Intergovernmental Panel on Climate Change. IPCC, Geneva, Switzerland, pp. 35–115. 33 Kakampakou, L., Simps...
2023
-
[2024]
Computational Statistics & Data Analysis 189, 107841
Joint modelling of the body and tail of bivariate data. Computational Statistics & Data Analysis 189, 107841. Ardizzone, L., Kruse, J., Rother, C., Köthe, U., 2019a. Analyzing inverse problems with invertibleneuralnetworks, in: Proceedingsofthe7thInternationalConferenceonLearn...
2019 arXiv
-
[2025]
ArXiv preprint arXiv:2503.23156
Neural Bayes inference for complex bivariate extremal dependence models. ArXiv preprint arXiv:2503.23156. André, L.M., Tawn, J.A.,
Reviewed August 5, 2026 · model on record in the stance chip above.
Discussion (0). Sign in to comment.