Coverage for src/gwtransport/gamma.py: 97%
67 statements
« prev ^ index » next coverage.py v7.15.0, created at 2026-07-07 21:17 +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.
33Available functions:
35- :func:`parse_parameters` - Parse and validate gamma distribution parameters from either
36 (mean, std, loc) or (alpha, beta, loc). Requires exactly one parameter pair and raises
37 ``ValueError`` if both are supplied; validates positivity and ordering constraints.
39- :func:`mean_std_loc_to_alpha_beta` - Convert physically intuitive (mean, std, loc) parameters
40 to gamma shape/scale parameters.
42- :func:`alpha_beta_loc_to_mean_std` - Convert gamma (alpha, beta, loc) parameters back to
43 (mean, std) for physical interpretation.
45- :func:`bins` - Primary function for transport modeling. Creates discrete probability bins
46 from the (optionally shifted) gamma distribution with equal-probability bins (default) or
47 custom quantile edges. Returns bin edges, expected values (mean pore volume within each
48 bin), and probability masses (weight in transport calculations).
50This file is part of gwtransport which is released under AGPL-3.0 license.
51See the ./LICENSE file or go to https://github.com/gwtransport/gwtransport/blob/main/LICENSE for full license details.
52"""
54import numpy as np
55import numpy.typing as npt
56from scipy.stats import gamma as gamma_dist
59def parse_parameters(
60 *,
61 mean: float | None = None,
62 std: float | None = None,
63 loc: float = 0.0,
64 alpha: float | None = None,
65 beta: float | None = None,
66) -> tuple[float, float, float]:
67 """
68 Parse parameters for gamma distribution.
70 Either ``(mean, std)`` or ``(alpha, beta)`` must be provided. ``loc`` is optional
71 and defaults to 0, which recovers the standard two-parameter gamma distribution.
73 Parameters
74 ----------
75 mean : float, optional
76 Mean of the gamma distribution. Must be strictly greater than ``loc``.
77 std : float, optional
78 Standard deviation of the gamma distribution. Must be positive. See
79 :ref:`concept-dispersion-scales` for what std represents depending
80 on APVD source. ``std`` is invariant under the ``loc`` shift.
81 loc : float, optional
82 Location (horizontal shift) of the gamma distribution; the lower bound of
83 support. Must satisfy ``loc >= 0`` and, when ``mean`` is supplied,
84 ``loc < mean``. Default is ``0.0``.
85 alpha : float, optional
86 Shape parameter of gamma distribution (must be > 0).
87 beta : float, optional
88 Scale parameter of gamma distribution (must be > 0).
90 Returns
91 -------
92 alpha : float
93 Shape parameter of gamma distribution.
94 beta : float
95 Scale parameter of gamma distribution.
96 loc : float
97 Location parameter of gamma distribution.
99 Raises
100 ------
101 ValueError
102 If neither ``(mean, std)`` nor ``(alpha, beta)`` is provided, if both pairs
103 are provided, if only one of a pair is provided, if ``alpha`` or ``beta`` are
104 not positive, if ``loc`` is negative, or if ``mean <= loc``.
105 """
106 if loc < 0:
107 msg = "loc must be non-negative"
108 raise ValueError(msg)
110 if (alpha is None) != (beta is None):
111 msg = "alpha and beta must both be provided or both be None."
112 raise ValueError(msg)
114 if alpha is not None and (mean is not None or std is not None):
115 msg = "Provide either (alpha, beta) or (mean, std), not both."
116 raise ValueError(msg)
118 # The ``or beta is None`` is redundant at runtime (the check above pairs them) but lets the
119 # type checker narrow ``beta`` to a float on the fall-through return.
120 if alpha is None or beta is None:
121 if mean is None or std is None:
122 msg = "Either (alpha, beta) or (mean, std) must be provided."
123 raise ValueError(msg)
124 # mean_std_loc_to_alpha_beta enforces std>0 and mean>loc, which together with
125 # loc>=0 guarantee alpha=(mean-loc)**2/std**2 > 0 and beta=std**2/(mean-loc) > 0.
126 alpha, beta = mean_std_loc_to_alpha_beta(mean=mean, std=std, loc=loc)
127 elif alpha <= 0 or beta <= 0:
128 msg = "Alpha and beta must be positive"
129 raise ValueError(msg)
131 return alpha, beta, loc
134def mean_std_loc_to_alpha_beta(*, mean: float, std: float, loc: float = 0.0) -> tuple[float, float]:
135 """
136 Convert mean, standard deviation, and location of gamma distribution to shape/scale.
138 The two-parameter shape/scale representation (``alpha``, ``beta``) is derived from
139 the excess-over-``loc`` moments: ``mean_excess = mean - loc``, ``std_excess = std``.
141 Parameters
142 ----------
143 mean : float
144 Mean of the gamma distribution. Must be strictly greater than ``loc``.
145 std : float
146 Standard deviation of the gamma distribution. Must be positive. See
147 :ref:`concept-dispersion-scales` for what std represents depending
148 on APVD source. ``std`` is invariant under the ``loc`` shift.
149 loc : float, optional
150 Location (horizontal shift) of the gamma distribution. Must satisfy
151 ``0 <= loc < mean``. Default is ``0.0``.
153 Returns
154 -------
155 alpha : float
156 Shape parameter of gamma distribution.
157 beta : float
158 Scale parameter of gamma distribution.
160 Raises
161 ------
162 ValueError
163 If ``std`` is not positive, if ``loc`` is negative, or if ``mean <= loc``.
165 See Also
166 --------
167 alpha_beta_loc_to_mean_std : Convert shape/scale/loc parameters to mean and std.
168 parse_parameters : Parse and validate gamma distribution parameters.
170 Examples
171 --------
172 >>> from gwtransport.gamma import mean_std_loc_to_alpha_beta
173 >>> mean_pore_volume = 30000.0 # m³
174 >>> std_pore_volume = 8100.0 # m³
175 >>> alpha, beta = mean_std_loc_to_alpha_beta(
176 ... mean=mean_pore_volume, std=std_pore_volume
177 ... )
178 >>> print(f"Shape parameter (alpha): {alpha:.2f}")
179 Shape parameter (alpha): 13.72
180 >>> print(f"Scale parameter (beta): {beta:.2f}")
181 Scale parameter (beta): 2187.00
183 With a 5000 m³ minimum pore volume:
185 >>> alpha, beta = mean_std_loc_to_alpha_beta(mean=30000.0, std=8100.0, loc=5000.0)
186 >>> print(f"Shape parameter (alpha): {alpha:.2f}")
187 Shape parameter (alpha): 9.53
188 >>> print(f"Scale parameter (beta): {beta:.2f}")
189 Scale parameter (beta): 2624.40
190 """
191 if std <= 0:
192 msg = "std must be positive"
193 raise ValueError(msg)
194 if loc < 0:
195 msg = "loc must be non-negative"
196 raise ValueError(msg)
197 if mean <= loc:
198 msg = "mean must be strictly greater than loc"
199 raise ValueError(msg)
201 mean_excess = mean - loc
202 alpha = mean_excess**2 / std**2
203 beta = std**2 / mean_excess
204 return alpha, beta
207def alpha_beta_loc_to_mean_std(*, alpha: float, beta: float, loc: float = 0.0) -> tuple[float, float]:
208 """
209 Convert shape, scale, and location of gamma distribution to mean and standard deviation.
211 Parameters are validated via :func:`parse_parameters`, which raises ``ValueError`` if
212 ``alpha`` or ``beta`` are non-positive or ``loc`` is negative.
214 Parameters
215 ----------
216 alpha : float
217 Shape parameter of the gamma distribution. Must be positive.
218 beta : float
219 Scale parameter of the gamma distribution. Must be positive.
220 loc : float, optional
221 Location (horizontal shift) of the gamma distribution. Must be non-negative.
222 Default is ``0.0``.
224 Returns
225 -------
226 mean : float
227 Mean of the gamma distribution, equal to ``alpha * beta + loc``.
228 std : float
229 Standard deviation of the gamma distribution, equal to ``sqrt(alpha) * beta``.
230 ``std`` is invariant under the ``loc`` shift.
232 See Also
233 --------
234 mean_std_loc_to_alpha_beta : Convert mean/std/loc to shape and scale parameters.
235 parse_parameters : Parse and validate gamma distribution parameters.
237 Examples
238 --------
239 >>> from gwtransport.gamma import alpha_beta_loc_to_mean_std
240 >>> alpha = 13.72 # shape parameter
241 >>> beta = 2187.0 # scale parameter
242 >>> mean, std = alpha_beta_loc_to_mean_std(alpha=alpha, beta=beta)
243 >>> print(f"Mean pore volume: {mean:.0f} m³")
244 Mean pore volume: 30006 m³
245 >>> print(f"Std pore volume: {std:.0f} m³")
246 Std pore volume: 8101 m³
247 """
248 parse_parameters(alpha=alpha, beta=beta, loc=loc)
249 return alpha * beta + loc, np.sqrt(alpha) * beta
252def bins(
253 *,
254 mean: float | None = None,
255 std: float | None = None,
256 loc: float = 0.0,
257 alpha: float | None = None,
258 beta: float | None = None,
259 n_bins: int = 100,
260 quantile_edges: npt.ArrayLike | None = None,
261) -> dict[str, npt.NDArray[np.floating]]:
262 """
263 Divide a (shifted) gamma distribution into bins and compute bin properties.
265 If ``n_bins`` is provided, the gamma distribution is divided into ``n_bins``
266 equal-mass bins. If ``quantile_edges`` is provided, the distribution is divided
267 into bins defined by those quantile edges. The quantile edges must be a strictly
268 increasing 1-D array of at least 3 entries (>= 2 bins) in ``[0, 1]``, with the
269 first and last entries exactly 0 and 1; ``n_bins`` is then ignored.
271 Parameters
272 ----------
273 mean : float, optional
274 Mean of the gamma distribution. Must be strictly greater than ``loc``.
275 std : float, optional
276 Standard deviation of the gamma distribution. Must be positive.
277 loc : float, optional
278 Location (horizontal shift) of the gamma distribution; the lower bound of
279 support. Must satisfy ``0 <= loc < mean`` (or ``loc >= 0`` when using
280 alpha/beta). Default is ``0.0``.
281 alpha : float, optional
282 Shape parameter of gamma distribution (must be > 0).
283 beta : float, optional
284 Scale parameter of gamma distribution (must be > 0).
285 n_bins : int, optional
286 Number of bins to divide the gamma distribution (must be >= 2). Default is 100.
287 quantile_edges : array-like, optional
288 Quantile edges for binning. Must be a strictly increasing 1-D array of at least
289 3 entries (>= 2 bins), all in ``[0, 1]``, with the first and last entries exactly
290 0 and 1. If provided, ``n_bins`` is ignored.
292 Returns
293 -------
294 dict
295 Dictionary with keys of type str and values of type numpy.ndarray:
297 - ``lower_bound``: lower bounds of bins (first one equals ``loc``)
298 - ``upper_bound``: upper bounds of bins (last one is inf)
299 - ``edges``: bin edges (lower_bound[0], upper_bound[0], ..., upper_bound[-1])
300 - ``expected_values``: expected values in bins. Is what you would expect to
301 observe if you repeatedly sampled from the probability distribution, but only
302 considered samples that fall within that particular bin.
303 - ``probability_mass``: probability mass in bins (invariant under ``loc`` shift).
305 Raises
306 ------
307 ValueError
308 If ``n_bins`` is not greater than 1, if ``quantile_edges`` is not a strictly
309 increasing 1-D array in ``[0, 1]`` with endpoints exactly 0 and 1, or if
310 parameter validation in :func:`parse_parameters` fails.
312 See Also
313 --------
314 mean_std_loc_to_alpha_beta : Convert mean/std/loc to alpha/beta parameters.
315 gwtransport.advection.gamma_infiltration_to_extraction : Use bins for transport modeling.
316 :ref:`concept-gamma-distribution` : Two-parameter pore volume model.
317 :ref:`concept-gamma-loc` : Shifted gamma with minimum pore volume.
318 :ref:`concept-dispersion-scales` : What ``std`` represents (macrodispersion vs total spreading).
319 :ref:`assumption-gamma-distribution` : When gamma distribution is adequate.
321 Examples
322 --------
323 Create equal-mass bins for a gamma distribution:
325 >>> from gwtransport.gamma import bins
326 >>> result = bins(mean=30000.0, std=8100.0, n_bins=5)
328 With a location parameter representing a minimum pore volume:
330 >>> result = bins(mean=30000.0, std=8100.0, loc=5000.0, n_bins=5)
331 >>> float(result["edges"][0])
332 5000.0
334 Create bins with custom quantile edges:
336 >>> import numpy as np
337 >>> quantiles = np.array([0.0, 0.25, 0.5, 0.75, 1.0])
338 >>> result = bins(mean=30000.0, std=8100.0, quantile_edges=quantiles)
339 >>> print(f"Number of bins: {len(result['probability_mass'])}")
340 Number of bins: 4
341 """
342 alpha, beta, loc = parse_parameters(mean=mean, std=std, loc=loc, alpha=alpha, beta=beta)
344 if quantile_edges is not None:
345 quantile_edges = np.asarray(quantile_edges, dtype=float)
346 if quantile_edges.ndim != 1:
347 msg = "quantile_edges must be a 1-D array."
348 raise ValueError(msg)
349 if quantile_edges.size == 0:
350 msg = "quantile_edges must not be empty."
351 raise ValueError(msg)
352 if not np.all(np.diff(quantile_edges) > 0):
353 msg = "quantile_edges must be strictly increasing."
354 raise ValueError(msg)
355 if quantile_edges[0] != 0.0 or quantile_edges[-1] != 1.0:
356 msg = "quantile_edges must start at 0 and end at 1."
357 raise ValueError(msg)
358 n_bins = len(quantile_edges) - 1
359 else:
360 quantile_edges = np.linspace(0, 1, n_bins + 1)
362 if n_bins <= 1:
363 msg = "Number of bins must be greater than 1"
364 raise ValueError(msg)
366 # Unshifted bin edges for the standard Gamma(alpha, beta) distribution, then shift
367 unshifted_edges = gamma_dist.ppf(quantile_edges, alpha, scale=beta)
368 bin_edges = unshifted_edges + loc
369 probability_mass = np.diff(quantile_edges) # probability mass for each bin
371 # Conditional mean within each bin for the unshifted distribution, then shift by loc.
372 # E[X | a <= X < b] for X ~ Gamma(alpha, beta) uses the identity
373 # E[X * 1_{a<=X<b}] = alpha * beta * (F_{alpha+1}(b/beta) - F_{alpha+1}(a/beta))
374 # where F_{alpha+1} is the CDF of Gamma(alpha+1, beta).
375 cdf_alpha_plus_1 = gamma_dist.cdf(unshifted_edges, alpha + 1, scale=beta)
376 diff_alpha_plus_1 = np.diff(cdf_alpha_plus_1)
377 expected_values = beta * alpha * diff_alpha_plus_1 / probability_mass + loc
379 return {
380 "lower_bound": bin_edges[:-1],
381 "upper_bound": bin_edges[1:],
382 "edges": bin_edges,
383 "expected_values": expected_values,
384 "probability_mass": probability_mass,
385 }