Coverage for src/gwtransport/gamma.py: 0%
63 statements
« prev ^ index » next coverage.py v7.15.3, created at 2026-08-04 20:54 +0000
1"""
2Gamma Distribution Utilities for Aquifer Pore Volume Heterogeneity.
4This module provides utilities for working with gamma distributions to model heterogeneous
5aquifer pore volumes in groundwater transport analysis. The gamma distribution offers a
6flexible three-parameter model (shape, scale, location) for representing the natural
7variability in flow path lengths and residence times within aquifer systems. In
8heterogeneous aquifers, water travels through multiple flow paths with different pore
9volumes; the location parameter additionally represents a guaranteed minimum pore volume
10(for example, immobile porosity or a geometric minimum travel distance).
12Parameterizations
13-----------------
14Two equivalent parameterizations are supported, each optionally with a location shift:
16- **(mean, std, loc)** — physically intuitive. ``mean`` is the total expected value,
17 ``std`` is the spread (invariant under shift), and ``loc`` is the lower bound of
18 support. Constraint: ``0 <= loc < mean``.
19- **(alpha, beta, loc)** — scipy-style. ``alpha`` is shape, ``beta`` is scale, and
20 ``loc`` is the lower bound of support. Constraint: ``alpha > 0``, ``beta > 0``,
21 ``loc >= 0``.
23Conversion formulas (with constraint ``mean > loc``):
25 alpha = ((mean - loc) / std) ** 2
26 beta = std ** 2 / (mean - loc)
27 mean = alpha * beta + loc
28 std = sqrt(alpha) * beta
30When ``loc == 0`` the three-parameter model reduces to the standard two-parameter
31gamma distribution.
33Streamtube discretization by :func:`bins` is always into bins of **equal probability mass**.
35Available functions:
37- :func:`parse_parameters` - Resolve either ``(mean, std)`` or ``(alpha, beta)``, together with an
38 optional ``loc``, into the validated ``(alpha, beta, loc)`` triple the rest of the package works with.
39 Exactly one of the two pairs must be supplied; positivity, ``0 <= loc < mean`` and finiteness are
40 enforced, otherwise ``ValueError`` is raised.
42- :func:`mean_std_loc_to_alpha_beta` - Convert the physically intuitive ``(mean, std, loc)`` parameters
43 to gamma shape and scale from the excess-over-``loc`` moments,
44 ``alpha = ((mean - loc) / std) ** 2`` and ``beta = std ** 2 / (mean - loc)``.
46- :func:`bins` - Split the (shifted) gamma distribution at the ``n_bins + 1`` uniform quantile edges, so
47 every bin (streamtube) carries probability mass ``1 / n_bins``. Returns a dict of arrays: the bin edges
48 (``lower_bound``, ``upper_bound``, ``edges``; the first lower bound is ``loc`` and the last upper bound
49 is infinite), the per-bin conditional mean pore volume (``expected_values``) that serves as the
50 streamtube pore volume in transport calculations, and the per-bin ``probability_mass``.
52This file is part of gwtransport which is released under AGPL-3.0 license.
53See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
54"""
56import numpy as np
57import numpy.typing as npt
58from scipy.stats import gamma as gamma_dist
60# Numerical-envelope guard threshold for bins(): when alpha*eps exceeds this fraction of the
61# equal-mass quantile gap 1/n_bins, alpha + 1 == alpha to machine precision and the per-bin
62# expected values degrade to noise.
63_HUGE_ALPHA_GAP_FRACTION = 0.01
66def parse_parameters(
67 *,
68 mean: float | None = None,
69 std: float | None = None,
70 loc: float = 0.0,
71 alpha: float | None = None,
72 beta: float | None = None,
73) -> tuple[float, float, float]:
74 """
75 Parse parameters for gamma distribution.
77 Either ``(mean, std)`` or ``(alpha, beta)`` must be provided. ``loc`` is optional
78 and defaults to 0, which recovers the standard two-parameter gamma distribution.
80 Parameters
81 ----------
82 mean : float, optional
83 Mean of the gamma distribution. Must be strictly greater than ``loc``.
84 std : float, optional
85 Standard deviation of the gamma distribution. Must be positive. See
86 :ref:`concept-dispersion-scales` for what std represents depending
87 on APVD source. ``std`` is invariant under the ``loc`` shift.
88 loc : float, optional
89 Location (horizontal shift) of the gamma distribution; the lower bound of
90 support. Must satisfy ``loc >= 0`` and, when ``mean`` is supplied,
91 ``loc < mean``. Default is ``0.0``.
92 alpha : float, optional
93 Shape parameter of gamma distribution (must be > 0).
94 beta : float, optional
95 Scale parameter of gamma distribution (must be > 0).
97 Returns
98 -------
99 alpha : float
100 Shape parameter of gamma distribution.
101 beta : float
102 Scale parameter of gamma distribution.
103 loc : float
104 Location parameter of gamma distribution.
106 Raises
107 ------
108 ValueError
109 If neither ``(mean, std)`` nor ``(alpha, beta)`` is provided, if both pairs
110 are provided, if only one of a pair is provided, if ``alpha`` or ``beta`` are
111 not positive, if ``loc`` is negative, if the resolved ``alpha``, ``beta``, or
112 ``loc`` is not finite, or if ``mean <= loc``.
113 """
114 if loc < 0:
115 msg = "loc must be non-negative"
116 raise ValueError(msg)
118 if (alpha is None) != (beta is None):
119 msg = "alpha and beta must both be provided or both be None."
120 raise ValueError(msg)
122 if alpha is not None and (mean is not None or std is not None):
123 msg = "Provide either (alpha, beta) or (mean, std), not both."
124 raise ValueError(msg)
126 if (mean is None) != (std is None):
127 msg = "mean and std must both be provided or both be None."
128 raise ValueError(msg)
130 # The ``or beta is None`` is redundant at runtime (the check above pairs them) but lets the
131 # type checker narrow ``beta`` to a float on the fall-through return.
132 if alpha is None or beta is None:
133 if mean is None or std is None:
134 msg = "Either (alpha, beta) or (mean, std) must be provided."
135 raise ValueError(msg)
136 # mean_std_loc_to_alpha_beta enforces std>0 and mean>loc, which together with
137 # loc>=0 guarantee alpha=(mean-loc)**2/std**2 > 0 and beta=std**2/(mean-loc) > 0.
138 alpha, beta = mean_std_loc_to_alpha_beta(mean=mean, std=std, loc=loc)
139 elif alpha <= 0 or beta <= 0:
140 msg = "Alpha and beta must be positive"
141 raise ValueError(msg)
143 # A non-finite alpha/beta/loc slips past the comparisons above (``nan <= 0`` and ``nan < 0`` are
144 # both False), producing an all-NaN distribution instead of a clear error. Reject it loudly.
145 if not (np.isfinite(alpha) and np.isfinite(beta) and np.isfinite(loc)):
146 msg = "alpha, beta, and loc must be finite."
147 raise ValueError(msg)
149 return alpha, beta, loc
152def mean_std_loc_to_alpha_beta(*, mean: float, std: float, loc: float = 0.0) -> tuple[float, float]:
153 """
154 Convert mean, standard deviation, and location of gamma distribution to shape/scale.
156 The two-parameter shape/scale representation (``alpha``, ``beta``) is derived from
157 the excess-over-``loc`` moments: ``mean_excess = mean - loc``, ``std_excess = std``.
159 Parameters
160 ----------
161 mean : float
162 Mean of the gamma distribution. Must be strictly greater than ``loc``.
163 std : float
164 Standard deviation of the gamma distribution. Must be positive. See
165 :ref:`concept-dispersion-scales` for what std represents depending
166 on APVD source. ``std`` is invariant under the ``loc`` shift.
167 loc : float, optional
168 Location (horizontal shift) of the gamma distribution. Must satisfy
169 ``0 <= loc < mean``. Default is ``0.0``.
171 Returns
172 -------
173 alpha : float
174 Shape parameter of gamma distribution.
175 beta : float
176 Scale parameter of gamma distribution.
178 Raises
179 ------
180 ValueError
181 If ``std`` is not positive, if ``loc`` is negative, or if ``mean <= loc``.
183 See Also
184 --------
185 parse_parameters : Parse and validate gamma distribution parameters.
187 Examples
188 --------
189 >>> from gwtransport.gamma import mean_std_loc_to_alpha_beta
190 >>> mean_pore_volume = 30000.0 # m³
191 >>> std_pore_volume = 8100.0 # m³
192 >>> alpha, beta = mean_std_loc_to_alpha_beta(
193 ... mean=mean_pore_volume, std=std_pore_volume
194 ... )
195 >>> print(f"Shape parameter (alpha): {alpha:.2f}")
196 Shape parameter (alpha): 13.72
197 >>> print(f"Scale parameter (beta): {beta:.2f}")
198 Scale parameter (beta): 2187.00
200 With a 5000 m³ minimum pore volume:
202 >>> alpha, beta = mean_std_loc_to_alpha_beta(mean=30000.0, std=8100.0, loc=5000.0)
203 >>> print(f"Shape parameter (alpha): {alpha:.2f}")
204 Shape parameter (alpha): 9.53
205 >>> print(f"Scale parameter (beta): {beta:.2f}")
206 Scale parameter (beta): 2624.40
207 """
208 if std <= 0:
209 msg = "std must be positive"
210 raise ValueError(msg)
211 if loc < 0:
212 msg = "loc must be non-negative"
213 raise ValueError(msg)
214 if mean <= loc:
215 msg = "mean must be strictly greater than loc"
216 raise ValueError(msg)
218 mean_excess = mean - loc
219 alpha = mean_excess**2 / std**2
220 beta = std**2 / mean_excess
221 return alpha, beta
224def bins(
225 *,
226 mean: float | None = None,
227 std: float | None = None,
228 loc: float = 0.0,
229 alpha: float | None = None,
230 beta: float | None = None,
231 n_bins: int = 100,
232) -> dict[str, npt.NDArray[np.floating]]:
233 """
234 Divide a (shifted) gamma distribution into equal-probability-mass bins and compute bin properties.
236 The distribution is split at the ``n_bins + 1`` uniform quantile edges, so every bin
237 (streamtube) carries probability mass ``1 / n_bins``.
239 Parameters
240 ----------
241 mean : float, optional
242 Mean of the gamma distribution. Must be strictly greater than ``loc``.
243 std : float, optional
244 Standard deviation of the gamma distribution. Must be positive.
245 loc : float, optional
246 Location (horizontal shift) of the gamma distribution; the lower bound of
247 support. Must satisfy ``0 <= loc < mean`` (or ``loc >= 0`` when using
248 alpha/beta). Default is ``0.0``.
249 alpha : float, optional
250 Shape parameter of gamma distribution (must be > 0).
251 beta : float, optional
252 Scale parameter of gamma distribution (must be > 0).
253 n_bins : int, optional
254 Number of bins to divide the gamma distribution (must be >= 2). Default is 100.
256 Returns
257 -------
258 dict
259 Dictionary with keys of type str and values of type numpy.ndarray:
261 - ``lower_bound``: lower bounds of bins (first one equals ``loc``)
262 - ``upper_bound``: upper bounds of bins (last one is inf)
263 - ``edges``: bin edges (lower_bound[0], upper_bound[0], ..., upper_bound[-1])
264 - ``expected_values``: expected values in bins. Is what you would expect to
265 observe if you repeatedly sampled from the probability distribution, but only
266 considered samples that fall within that particular bin.
267 - ``probability_mass``: probability mass in bins (invariant under ``loc`` shift).
269 Raises
270 ------
271 ValueError
272 If ``n_bins`` is not greater than 1, or if parameter validation in
273 :func:`parse_parameters` fails. Also raised for numerically-degenerate requests that
274 would otherwise return silently-wrong structure: an ``alpha`` so large that
275 ``alpha + 1 == alpha`` in float64 relative to the ``1 / n_bins`` quantile gap (the
276 distribution is numerically a point mass), or a bin whose expected value underflows
277 to ``loc``.
279 See Also
280 --------
281 mean_std_loc_to_alpha_beta : Convert mean/std/loc to alpha/beta parameters.
282 gwtransport.advection.gamma_infiltration_to_extraction : Use bins for transport modeling.
283 :ref:`concept-gamma-distribution` : Two-parameter pore volume model.
284 :ref:`concept-gamma-loc` : Shifted gamma with minimum pore volume.
285 :ref:`concept-dispersion-scales` : What ``std`` represents (macrodispersion vs total spreading).
286 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate.
288 Examples
289 --------
290 Create equal-mass bins for a gamma distribution:
292 >>> from gwtransport.gamma import bins
293 >>> result = bins(mean=30000.0, std=8100.0, n_bins=5)
294 >>> print(f"Number of bins: {len(result['probability_mass'])}")
295 Number of bins: 5
297 With a location parameter representing a minimum pore volume:
299 >>> result = bins(mean=30000.0, std=8100.0, loc=5000.0, n_bins=5)
300 >>> float(result["edges"][0])
301 5000.0
302 """
303 alpha, beta, loc = parse_parameters(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta)
305 if n_bins <= 1:
306 # Validate before np.linspace: a negative n_bins would otherwise surface as
307 # numpy's opaque "Number of samples ... must be non-negative" error.
308 msg = "Number of bins must be greater than 1"
309 raise ValueError(msg)
311 quantile_edges = np.linspace(0, 1, n_bins + 1)
313 # Guard the numerical cliff of the closed-form conditional mean below: an alpha so large that
314 # alpha*eps exceeds ~1% of the equal-mass quantile gap 1/n_bins makes alpha+1 == alpha to machine
315 # precision, so the conditional means degrade to noise. Only reachable with a near-degenerate
316 # (delta-like) distribution.
317 if alpha * np.finfo(float).eps * n_bins > _HUGE_ALPHA_GAP_FRACTION:
318 msg = (
319 f"alpha ({alpha:.3g}) is too large for float64 bin resolution: alpha*eps exceeds 1% of "
320 f"the equal-mass quantile gap (1/{n_bins}), so alpha+1 == alpha to machine precision and "
321 "the per-bin expected values are numerical noise. The distribution is effectively a point "
322 "mass at alpha*beta + loc; use a larger std/(mean-loc) or fewer bins."
323 )
324 raise ValueError(msg)
326 # Unshifted bin edges for the standard Gamma(alpha, beta) distribution, then shift
327 unshifted_edges = gamma_dist.ppf(quantile_edges, alpha, scale=beta)
328 bin_edges = unshifted_edges + loc
329 probability_mass = np.diff(quantile_edges) # probability mass for each bin
331 # Conditional mean within each bin for the unshifted distribution, then shift by loc.
332 # E[X | a <= X < b] for X ~ Gamma(alpha, beta) uses the identity
333 # E[X * 1_{a<=X<b}] = alpha * beta * (F_{alpha+1}(b) - F_{alpha+1}(a))
334 # where F_{alpha+1} is the CDF of Gamma(alpha+1, scale=beta) (equivalently the regularized
335 # lower incomplete gamma P(alpha+1, b/beta) - P(alpha+1, a/beta)).
336 cdf_alpha_plus_1 = gamma_dist.cdf(unshifted_edges, alpha + 1, scale=beta)
337 diff_alpha_plus_1 = np.diff(cdf_alpha_plus_1)
339 # Pre-shift conditional mean of the excess over loc. Every positive-mass bin of a Gamma(alpha, beta)
340 # has a strictly positive conditional mean, so this must be > 0; a value <= 0 signals an underflow /
341 # cancellation of the CDF difference (very small alpha or extremely fine bins) that would emit a
342 # numerically-zero or negative pore volume. Test this pre-shift quantity rather than the shifted
343 # expected value: for loc > 0 a benign ``loc + tiny_positive_excess == loc`` rounding would otherwise
344 # be misread as underflow and reject correct, usable output for shifted heterogeneous APVDs.
345 cond_mean_excess = beta * alpha * diff_alpha_plus_1 / probability_mass
346 if np.any(cond_mean_excess <= 0.0):
347 msg = (
348 "A bin's conditional expected value underflowed to loc (its excess conditional mean over loc "
349 "is not strictly positive). This happens for very small alpha or extremely fine bins where "
350 "the CDF difference underflows. Use fewer bins or a larger alpha."
351 )
352 raise ValueError(msg)
354 expected_values = cond_mean_excess + loc
356 return {
357 "lower_bound": bin_edges[:-1],
358 "upper_bound": bin_edges[1:],
359 "edges": bin_edges,
360 "expected_values": expected_values,
361 "probability_mass": probability_mass,
362 }