Coverage for /usr/lib/python3/dist-packages/scipy/stats/_qmvnt.py: 6%
250 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1# Integration of multivariate normal and t distributions.
3# Adapted from the MATLAB original implementations by Dr. Alan Genz.
5# http://www.math.wsu.edu/faculty/genz/software/software.html
7# Copyright (C) 2013, Alan Genz, All rights reserved.
8# Python implementation is copyright (C) 2022, Robert Kern, All rights
9# reserved.
11# Redistribution and use in source and binary forms, with or without
12# modification, are permitted provided the following conditions are met:
13# 1. Redistributions of source code must retain the above copyright
14# notice, this list of conditions and the following disclaimer.
15# 2. Redistributions in binary form must reproduce the above copyright
16# notice, this list of conditions and the following disclaimer in
17# the documentation and/or other materials provided with the
18# distribution.
19# 3. The contributor name(s) may not be used to endorse or promote
20# products derived from this software without specific prior
21# written permission.
22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
25# FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
26# COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
27# INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
28# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS
29# OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
30# ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR
31# TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF USE
32# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
35import numpy as np
37from scipy.fft import fft, ifft
38from scipy.special import gammaincinv, ndtr, ndtri
39from scipy.stats._qmc import primes_from_2_to
42phi = ndtr
43phinv = ndtri
46def _factorize_int(n):
47 """Return a sorted list of the unique prime factors of a positive integer.
48 """
49 # NOTE: There are lots faster ways to do this, but this isn't terrible.
50 factors = set()
51 for p in primes_from_2_to(int(np.sqrt(n)) + 1):
52 while not (n % p):
53 factors.add(p)
54 n //= p
55 if n == 1:
56 break
57 if n != 1:
58 factors.add(n)
59 return sorted(factors)
62def _primitive_root(p):
63 """Compute a primitive root of the prime number `p`.
65 Used in the CBC lattice construction.
67 References
68 ----------
69 .. [1] https://en.wikipedia.org/wiki/Primitive_root_modulo_n
70 """
71 # p is prime
72 pm = p - 1
73 factors = _factorize_int(pm)
74 n = len(factors)
75 r = 2
76 k = 0
77 while k < n:
78 d = pm // factors[k]
79 # pow() doesn't like numpy scalar types.
80 rd = pow(int(r), int(d), int(p))
81 if rd == 1:
82 r += 1
83 k = 0
84 else:
85 k += 1
86 return r
89def _cbc_lattice(n_dim, n_qmc_samples):
90 """Compute a QMC lattice generator using a Fast CBC construction.
92 Parameters
93 ----------
94 n_dim : int > 0
95 The number of dimensions for the lattice.
96 n_qmc_samples : int > 0
97 The desired number of QMC samples. This will be rounded down to the
98 nearest prime to enable the CBC construction.
100 Returns
101 -------
102 q : float array : shape=(n_dim,)
103 The lattice generator vector. All values are in the open interval
104 `(0, 1)`.
105 actual_n_qmc_samples : int
106 The prime number of QMC samples that must be used with this lattice,
107 no more, no less.
109 References
110 ----------
111 .. [1] Nuyens, D. and Cools, R. "Fast Component-by-Component Construction,
112 a Reprise for Different Kernels", In H. Niederreiter and D. Talay,
113 editors, Monte-Carlo and Quasi-Monte Carlo Methods 2004,
114 Springer-Verlag, 2006, 371-385.
115 """
116 # Round down to the nearest prime number.
117 primes = primes_from_2_to(n_qmc_samples + 1)
118 n_qmc_samples = primes[-1]
120 bt = np.ones(n_dim)
121 gm = np.hstack([1.0, 0.8 ** np.arange(n_dim - 1)])
122 q = 1
123 w = 0
124 z = np.arange(1, n_dim + 1)
125 m = (n_qmc_samples - 1) // 2
126 g = _primitive_root(n_qmc_samples)
127 # Slightly faster way to compute perm[j] = pow(g, j, n_qmc_samples)
128 # Shame that we don't have modulo pow() implemented as a ufunc.
129 perm = np.ones(m, dtype=int)
130 for j in range(m - 1):
131 perm[j + 1] = (g * perm[j]) % n_qmc_samples
132 perm = np.minimum(n_qmc_samples - perm, perm)
133 pn = perm / n_qmc_samples
134 c = pn * pn - pn + 1.0 / 6
135 fc = fft(c)
136 for s in range(1, n_dim):
137 reordered = np.hstack([
138 c[:w+1][::-1],
139 c[w+1:m][::-1],
140 ])
141 q = q * (bt[s-1] + gm[s-1] * reordered)
142 w = ifft(fc * fft(q)).real.argmin()
143 z[s] = perm[w]
144 q = z / n_qmc_samples
145 return q, n_qmc_samples
148# Note: this function is not currently used or tested by any SciPy code. It is
149# included in this file to facilitate the development of a parameter for users
150# to set the desired CDF accuracy, but must be reviewed and tested before use.
151def _qauto(func, covar, low, high, rng, error=1e-3, limit=10_000, **kwds):
152 """Automatically rerun the integration to get the required error bound.
154 Parameters
155 ----------
156 func : callable
157 Either :func:`_qmvn` or :func:`_qmvt`.
158 covar, low, high : array
159 As specified in :func:`_qmvn` and :func:`_qmvt`.
160 rng : Generator, optional
161 default_rng(), yada, yada
162 error : float > 0
163 The desired error bound.
164 limit : int > 0:
165 The rough limit of the number of integration points to consider. The
166 integration will stop looping once this limit has been *exceeded*.
167 **kwds :
168 Other keyword arguments to pass to `func`. When using :func:`_qmvt`, be
169 sure to include ``nu=`` as one of these.
171 Returns
172 -------
173 prob : float
174 The estimated probability mass within the bounds.
175 est_error : float
176 3 times the standard error of the batch estimates.
177 n_samples : int
178 The number of integration points actually used.
179 """
180 n = len(covar)
181 n_samples = 0
182 if n == 1:
183 prob = phi(high) - phi(low)
184 # More or less
185 est_error = 1e-15
186 else:
187 mi = min(limit, n * 1000)
188 prob = 0.0
189 est_error = 1.0
190 ei = 0.0
191 while est_error > error and n_samples < limit:
192 mi = round(np.sqrt(2) * mi)
193 pi, ei, ni = func(mi, covar, low, high, rng=rng, **kwds)
194 n_samples += ni
195 wt = 1.0 / (1 + (ei / est_error)**2)
196 prob += wt * (pi - prob)
197 est_error = np.sqrt(wt) * ei
198 return prob, est_error, n_samples
201# Note: this function is not currently used or tested by any SciPy code. It is
202# included in this file to facilitate the resolution of gh-8367, gh-16142, and
203# possibly gh-14286, but must be reviewed and tested before use.
204def _qmvn(m, covar, low, high, rng, lattice='cbc', n_batches=10):
205 """Multivariate normal integration over box bounds.
207 Parameters
208 ----------
209 m : int > n_batches
210 The number of points to sample. This number will be divided into
211 `n_batches` batches that apply random offsets of the sampling lattice
212 for each batch in order to estimate the error.
213 covar : (n, n) float array
214 Possibly singular, positive semidefinite symmetric covariance matrix.
215 low, high : (n,) float array
216 The low and high integration bounds.
217 rng : Generator, optional
218 default_rng(), yada, yada
219 lattice : 'cbc' or callable
220 The type of lattice rule to use to construct the integration points.
221 n_batches : int > 0, optional
222 The number of QMC batches to apply.
224 Returns
225 -------
226 prob : float
227 The estimated probability mass within the bounds.
228 est_error : float
229 3 times the standard error of the batch estimates.
230 """
231 cho, lo, hi = _permuted_cholesky(covar, low, high)
232 n = cho.shape[0]
233 ct = cho[0, 0]
234 c = phi(lo[0] / ct)
235 d = phi(hi[0] / ct)
236 ci = c
237 dci = d - ci
238 prob = 0.0
239 error_var = 0.0
240 q, n_qmc_samples = _cbc_lattice(n - 1, max(m // n_batches, 1))
241 y = np.zeros((n - 1, n_qmc_samples))
242 i_samples = np.arange(n_qmc_samples) + 1
243 for j in range(n_batches):
244 c = np.full(n_qmc_samples, ci)
245 dc = np.full(n_qmc_samples, dci)
246 pv = dc.copy()
247 for i in range(1, n):
248 # Pseudorandomly-shifted lattice coordinate.
249 z = q[i - 1] * i_samples + rng.random()
250 # Fast remainder(z, 1.0)
251 z -= z.astype(int)
252 # Tent periodization transform.
253 x = abs(2 * z - 1)
254 y[i - 1, :] = phinv(c + x * dc)
255 s = cho[i, :i] @ y[:i, :]
256 ct = cho[i, i]
257 c = phi((lo[i] - s) / ct)
258 d = phi((hi[i] - s) / ct)
259 dc = d - c
260 pv = pv * dc
261 # Accumulate the mean and error variances with online formulations.
262 d = (pv.mean() - prob) / (j + 1)
263 prob += d
264 error_var = (j - 1) * error_var / (j + 1) + d * d
265 # Error bounds are 3 times the standard error of the estimates.
266 est_error = 3 * np.sqrt(error_var)
267 n_samples = n_qmc_samples * n_batches
268 return prob, est_error, n_samples
271# Note: this function is not currently used or tested by any SciPy code. It is
272# included in this file to facilitate the resolution of gh-8367, gh-16142, and
273# possibly gh-14286, but must be reviewed and tested before use.
274def _mvn_qmc_integrand(covar, low, high, use_tent=False):
275 """Transform the multivariate normal integration into a QMC integrand over
276 a unit hypercube.
278 The dimensionality of the resulting hypercube integration domain is one
279 less than the dimensionality of the original integrand. Note that this
280 transformation subsumes the integration bounds in order to account for
281 infinite bounds. The QMC integration one does with the returned integrand
282 should be on the unit hypercube.
284 Parameters
285 ----------
286 covar : (n, n) float array
287 Possibly singular, positive semidefinite symmetric covariance matrix.
288 low, high : (n,) float array
289 The low and high integration bounds.
290 use_tent : bool, optional
291 If True, then use tent periodization. Only helpful for lattice rules.
293 Returns
294 -------
295 integrand : Callable[[NDArray], NDArray]
296 The QMC-integrable integrand. It takes an
297 ``(n_qmc_samples, ndim_integrand)`` array of QMC samples in the unit
298 hypercube and returns the ``(n_qmc_samples,)`` evaluations of at these
299 QMC points.
300 ndim_integrand : int
301 The dimensionality of the integrand. Equal to ``n-1``.
302 """
303 cho, lo, hi = _permuted_cholesky(covar, low, high)
304 n = cho.shape[0]
305 ndim_integrand = n - 1
306 ct = cho[0, 0]
307 c = phi(lo[0] / ct)
308 d = phi(hi[0] / ct)
309 ci = c
310 dci = d - ci
312 def integrand(*zs):
313 ndim_qmc = len(zs)
314 n_qmc_samples = len(np.atleast_1d(zs[0]))
315 assert ndim_qmc == ndim_integrand
316 y = np.zeros((ndim_qmc, n_qmc_samples))
317 c = np.full(n_qmc_samples, ci)
318 dc = np.full(n_qmc_samples, dci)
319 pv = dc.copy()
320 for i in range(1, n):
321 if use_tent:
322 # Tent periodization transform.
323 x = abs(2 * zs[i-1] - 1)
324 else:
325 x = zs[i-1]
326 y[i - 1, :] = phinv(c + x * dc)
327 s = cho[i, :i] @ y[:i, :]
328 ct = cho[i, i]
329 c = phi((lo[i] - s) / ct)
330 d = phi((hi[i] - s) / ct)
331 dc = d - c
332 pv = pv * dc
333 return pv
335 return integrand, ndim_integrand
338def _qmvt(m, nu, covar, low, high, rng, lattice='cbc', n_batches=10):
339 """Multivariate t integration over box bounds.
341 Parameters
342 ----------
343 m : int > n_batches
344 The number of points to sample. This number will be divided into
345 `n_batches` batches that apply random offsets of the sampling lattice
346 for each batch in order to estimate the error.
347 nu : float >= 0
348 The shape parameter of the multivariate t distribution.
349 covar : (n, n) float array
350 Possibly singular, positive semidefinite symmetric covariance matrix.
351 low, high : (n,) float array
352 The low and high integration bounds.
353 rng : Generator, optional
354 default_rng(), yada, yada
355 lattice : 'cbc' or callable
356 The type of lattice rule to use to construct the integration points.
357 n_batches : int > 0, optional
358 The number of QMC batches to apply.
360 Returns
361 -------
362 prob : float
363 The estimated probability mass within the bounds.
364 est_error : float
365 3 times the standard error of the batch estimates.
366 n_samples : int
367 The number of samples actually used.
368 """
369 sn = max(1.0, np.sqrt(nu))
370 low = np.asarray(low, dtype=np.float64)
371 high = np.asarray(high, dtype=np.float64)
372 cho, lo, hi = _permuted_cholesky(covar, low / sn, high / sn)
373 n = cho.shape[0]
374 prob = 0.0
375 error_var = 0.0
376 q, n_qmc_samples = _cbc_lattice(n, max(m // n_batches, 1))
377 i_samples = np.arange(n_qmc_samples) + 1
378 for j in range(n_batches):
379 pv = np.ones(n_qmc_samples)
380 s = np.zeros((n, n_qmc_samples))
381 for i in range(n):
382 # Pseudorandomly-shifted lattice coordinate.
383 z = q[i] * i_samples + rng.random()
384 # Fast remainder(z, 1.0)
385 z -= z.astype(int)
386 # Tent periodization transform.
387 x = abs(2 * z - 1)
388 # FIXME: Lift the i==0 case out of the loop to make the logic
389 # easier to follow.
390 if i == 0:
391 # We'll use one of the QR variates to pull out the
392 # t-distribution scaling.
393 if nu > 0:
394 r = np.sqrt(2 * gammaincinv(nu / 2, x))
395 else:
396 r = np.ones_like(x)
397 else:
398 y = phinv(c + x * dc) # noqa: F821
399 with np.errstate(invalid='ignore'):
400 s[i:, :] += cho[i:, i - 1][:, np.newaxis] * y
401 si = s[i, :]
403 c = np.ones(n_qmc_samples)
404 d = np.ones(n_qmc_samples)
405 with np.errstate(invalid='ignore'):
406 lois = lo[i] * r - si
407 hiis = hi[i] * r - si
408 c[lois < -9] = 0.0
409 d[hiis < -9] = 0.0
410 lo_mask = abs(lois) < 9
411 hi_mask = abs(hiis) < 9
412 c[lo_mask] = phi(lois[lo_mask])
413 d[hi_mask] = phi(hiis[hi_mask])
415 dc = d - c
416 pv *= dc
418 # Accumulate the mean and error variances with online formulations.
419 d = (pv.mean() - prob) / (j + 1)
420 prob += d
421 error_var = (j - 1) * error_var / (j + 1) + d * d
422 # Error bounds are 3 times the standard error of the estimates.
423 est_error = 3 * np.sqrt(error_var)
424 n_samples = n_qmc_samples * n_batches
425 return prob, est_error, n_samples
428def _permuted_cholesky(covar, low, high, tol=1e-10):
429 """Compute a scaled, permuted Cholesky factor, with integration bounds.
431 The scaling and permuting of the dimensions accomplishes part of the
432 transformation of the original integration problem into a more numerically
433 tractable form. The lower-triangular Cholesky factor will then be used in
434 the subsequent integration. The integration bounds will be scaled and
435 permuted as well.
437 Parameters
438 ----------
439 covar : (n, n) float array
440 Possibly singular, positive semidefinite symmetric covariance matrix.
441 low, high : (n,) float array
442 The low and high integration bounds.
443 tol : float, optional
444 The singularity tolerance.
446 Returns
447 -------
448 cho : (n, n) float array
449 Lower Cholesky factor, scaled and permuted.
450 new_low, new_high : (n,) float array
451 The scaled and permuted low and high integration bounds.
452 """
453 # Make copies for outputting.
454 cho = np.array(covar, dtype=np.float64)
455 new_lo = np.array(low, dtype=np.float64)
456 new_hi = np.array(high, dtype=np.float64)
457 n = cho.shape[0]
458 if cho.shape != (n, n):
459 raise ValueError("expected a square symmetric array")
460 if new_lo.shape != (n,) or new_hi.shape != (n,):
461 raise ValueError(
462 "expected integration boundaries the same dimensions "
463 "as the covariance matrix"
464 )
465 # Scale by the sqrt of the diagonal.
466 dc = np.sqrt(np.maximum(np.diag(cho), 0.0))
467 # But don't divide by 0.
468 dc[dc == 0.0] = 1.0
469 new_lo /= dc
470 new_hi /= dc
471 cho /= dc
472 cho /= dc[:, np.newaxis]
474 y = np.zeros(n)
475 sqtp = np.sqrt(2 * np.pi)
476 for k in range(n):
477 epk = (k + 1) * tol
478 im = k
479 ck = 0.0
480 dem = 1.0
481 s = 0.0
482 lo_m = 0.0
483 hi_m = 0.0
484 for i in range(k, n):
485 if cho[i, i] > tol:
486 ci = np.sqrt(cho[i, i])
487 if i > 0:
488 s = cho[i, :k] @ y[:k]
489 lo_i = (new_lo[i] - s) / ci
490 hi_i = (new_hi[i] - s) / ci
491 de = phi(hi_i) - phi(lo_i)
492 if de <= dem:
493 ck = ci
494 dem = de
495 lo_m = lo_i
496 hi_m = hi_i
497 im = i
498 if im > k:
499 # Swap im and k
500 cho[im, im] = cho[k, k]
501 _swap_slices(cho, np.s_[im, :k], np.s_[k, :k])
502 _swap_slices(cho, np.s_[im + 1:, im], np.s_[im + 1:, k])
503 _swap_slices(cho, np.s_[k + 1:im, k], np.s_[im, k + 1:im])
504 _swap_slices(new_lo, k, im)
505 _swap_slices(new_hi, k, im)
506 if ck > epk:
507 cho[k, k] = ck
508 cho[k, k + 1:] = 0.0
509 for i in range(k + 1, n):
510 cho[i, k] /= ck
511 cho[i, k + 1:i + 1] -= cho[i, k] * cho[k + 1:i + 1, k]
512 if abs(dem) > tol:
513 y[k] = ((np.exp(-lo_m * lo_m / 2) - np.exp(-hi_m * hi_m / 2)) /
514 (sqtp * dem))
515 else:
516 y[k] = (lo_m + hi_m) / 2
517 if lo_m < -10:
518 y[k] = hi_m
519 elif hi_m > 10:
520 y[k] = lo_m
521 cho[k, :k + 1] /= ck
522 new_lo[k] /= ck
523 new_hi[k] /= ck
524 else:
525 cho[k:, k] = 0.0
526 y[k] = (new_lo[k] + new_hi[k]) / 2
527 return cho, new_lo, new_hi
530def _swap_slices(x, slc1, slc2):
531 t = x[slc1].copy()
532 x[slc1] = x[slc2].copy()
533 x[slc2] = t