Coverage for /usr/lib/python3/dist-packages/scipy/stats/_multivariate.py: 27%
1787 statements
« prev ^ index » next coverage.py v7.9.1, created at 2025-06-14 15:55 +0200
1#
2# Author: Joris Vankerschaver 2013
3#
4import math
5import numpy as np
6from numpy import asarray_chkfinite, asarray
7from numpy.lib import NumpyVersion
8import scipy.linalg
9from scipy._lib import doccer
10from scipy.special import (gammaln, psi, multigammaln, xlogy, entr, betaln,
11 ive, loggamma)
12from scipy._lib._util import check_random_state
13from scipy.linalg.blas import drot
14from scipy.linalg._misc import LinAlgError
15from scipy.linalg.lapack import get_lapack_funcs
16from ._discrete_distns import binom
17from . import _mvn, _covariance, _rcont
18from ._qmvnt import _qmvt
19from ._morestats import directional_stats
20from scipy.optimize import root_scalar
22__all__ = ['multivariate_normal',
23 'matrix_normal',
24 'dirichlet',
25 'dirichlet_multinomial',
26 'wishart',
27 'invwishart',
28 'multinomial',
29 'special_ortho_group',
30 'ortho_group',
31 'random_correlation',
32 'unitary_group',
33 'multivariate_t',
34 'multivariate_hypergeom',
35 'random_table',
36 'uniform_direction',
37 'vonmises_fisher']
39_LOG_2PI = np.log(2 * np.pi)
40_LOG_2 = np.log(2)
41_LOG_PI = np.log(np.pi)
44_doc_random_state = """\
45seed : {None, int, np.random.RandomState, np.random.Generator}, optional
46 Used for drawing random variates.
47 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
48 If `seed` is an int, a new ``RandomState`` instance is used, seeded
49 with seed.
50 If `seed` is already a ``RandomState`` or ``Generator`` instance,
51 then that object is used.
52 Default is `None`.
53"""
56def _squeeze_output(out):
57 """
58 Remove single-dimensional entries from array and convert to scalar,
59 if necessary.
60 """
61 out = out.squeeze()
62 if out.ndim == 0:
63 out = out[()]
64 return out
67def _eigvalsh_to_eps(spectrum, cond=None, rcond=None):
68 """Determine which eigenvalues are "small" given the spectrum.
70 This is for compatibility across various linear algebra functions
71 that should agree about whether or not a Hermitian matrix is numerically
72 singular and what is its numerical matrix rank.
73 This is designed to be compatible with scipy.linalg.pinvh.
75 Parameters
76 ----------
77 spectrum : 1d ndarray
78 Array of eigenvalues of a Hermitian matrix.
79 cond, rcond : float, optional
80 Cutoff for small eigenvalues.
81 Singular values smaller than rcond * largest_eigenvalue are
82 considered zero.
83 If None or -1, suitable machine precision is used.
85 Returns
86 -------
87 eps : float
88 Magnitude cutoff for numerical negligibility.
90 """
91 if rcond is not None:
92 cond = rcond
93 if cond in [None, -1]:
94 t = spectrum.dtype.char.lower()
95 factor = {'f': 1E3, 'd': 1E6}
96 cond = factor[t] * np.finfo(t).eps
97 eps = cond * np.max(abs(spectrum))
98 return eps
101def _pinv_1d(v, eps=1e-5):
102 """A helper function for computing the pseudoinverse.
104 Parameters
105 ----------
106 v : iterable of numbers
107 This may be thought of as a vector of eigenvalues or singular values.
108 eps : float
109 Values with magnitude no greater than eps are considered negligible.
111 Returns
112 -------
113 v_pinv : 1d float ndarray
114 A vector of pseudo-inverted numbers.
116 """
117 return np.array([0 if abs(x) <= eps else 1/x for x in v], dtype=float)
120class _PSD:
121 """
122 Compute coordinated functions of a symmetric positive semidefinite matrix.
124 This class addresses two issues. Firstly it allows the pseudoinverse,
125 the logarithm of the pseudo-determinant, and the rank of the matrix
126 to be computed using one call to eigh instead of three.
127 Secondly it allows these functions to be computed in a way
128 that gives mutually compatible results.
129 All of the functions are computed with a common understanding as to
130 which of the eigenvalues are to be considered negligibly small.
131 The functions are designed to coordinate with scipy.linalg.pinvh()
132 but not necessarily with np.linalg.det() or with np.linalg.matrix_rank().
134 Parameters
135 ----------
136 M : array_like
137 Symmetric positive semidefinite matrix (2-D).
138 cond, rcond : float, optional
139 Cutoff for small eigenvalues.
140 Singular values smaller than rcond * largest_eigenvalue are
141 considered zero.
142 If None or -1, suitable machine precision is used.
143 lower : bool, optional
144 Whether the pertinent array data is taken from the lower
145 or upper triangle of M. (Default: lower)
146 check_finite : bool, optional
147 Whether to check that the input matrices contain only finite
148 numbers. Disabling may give a performance gain, but may result
149 in problems (crashes, non-termination) if the inputs do contain
150 infinities or NaNs.
151 allow_singular : bool, optional
152 Whether to allow a singular matrix. (Default: True)
154 Notes
155 -----
156 The arguments are similar to those of scipy.linalg.pinvh().
158 """
160 def __init__(self, M, cond=None, rcond=None, lower=True,
161 check_finite=True, allow_singular=True):
162 self._M = np.asarray(M)
164 # Compute the symmetric eigendecomposition.
165 # Note that eigh takes care of array conversion, chkfinite,
166 # and assertion that the matrix is square.
167 s, u = scipy.linalg.eigh(M, lower=lower, check_finite=check_finite)
169 eps = _eigvalsh_to_eps(s, cond, rcond)
170 if np.min(s) < -eps:
171 msg = "The input matrix must be symmetric positive semidefinite."
172 raise ValueError(msg)
173 d = s[s > eps]
174 if len(d) < len(s) and not allow_singular:
175 msg = ("When `allow_singular is False`, the input matrix must be "
176 "symmetric positive definite.")
177 raise np.linalg.LinAlgError(msg)
178 s_pinv = _pinv_1d(s, eps)
179 U = np.multiply(u, np.sqrt(s_pinv))
181 # Save the eigenvector basis, and tolerance for testing support
182 self.eps = 1e3*eps
183 self.V = u[:, s <= eps]
185 # Initialize the eagerly precomputed attributes.
186 self.rank = len(d)
187 self.U = U
188 self.log_pdet = np.sum(np.log(d))
190 # Initialize attributes to be lazily computed.
191 self._pinv = None
193 def _support_mask(self, x):
194 """
195 Check whether x lies in the support of the distribution.
196 """
197 residual = np.linalg.norm(x @ self.V, axis=-1)
198 in_support = residual < self.eps
199 return in_support
201 @property
202 def pinv(self):
203 if self._pinv is None:
204 self._pinv = np.dot(self.U, self.U.T)
205 return self._pinv
208class multi_rv_generic:
209 """
210 Class which encapsulates common functionality between all multivariate
211 distributions.
212 """
213 def __init__(self, seed=None):
214 super().__init__()
215 self._random_state = check_random_state(seed)
217 @property
218 def random_state(self):
219 """ Get or set the Generator object for generating random variates.
221 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
222 singleton is used.
223 If `seed` is an int, a new ``RandomState`` instance is used,
224 seeded with `seed`.
225 If `seed` is already a ``Generator`` or ``RandomState`` instance then
226 that instance is used.
228 """
229 return self._random_state
231 @random_state.setter
232 def random_state(self, seed):
233 self._random_state = check_random_state(seed)
235 def _get_random_state(self, random_state):
236 if random_state is not None:
237 return check_random_state(random_state)
238 else:
239 return self._random_state
242class multi_rv_frozen:
243 """
244 Class which encapsulates common functionality between all frozen
245 multivariate distributions.
246 """
247 @property
248 def random_state(self):
249 return self._dist._random_state
251 @random_state.setter
252 def random_state(self, seed):
253 self._dist._random_state = check_random_state(seed)
256_mvn_doc_default_callparams = """\
257mean : array_like, default: ``[0]``
258 Mean of the distribution.
259cov : array_like or `Covariance`, default: ``[1]``
260 Symmetric positive (semi)definite covariance matrix of the distribution.
261allow_singular : bool, default: ``False``
262 Whether to allow a singular covariance matrix. This is ignored if `cov` is
263 a `Covariance` object.
264"""
266_mvn_doc_callparams_note = """\
267Setting the parameter `mean` to `None` is equivalent to having `mean`
268be the zero-vector. The parameter `cov` can be a scalar, in which case
269the covariance matrix is the identity times that value, a vector of
270diagonal entries for the covariance matrix, a two-dimensional array_like,
271or a `Covariance` object.
272"""
274_mvn_doc_frozen_callparams = ""
276_mvn_doc_frozen_callparams_note = """\
277See class definition for a detailed description of parameters."""
279mvn_docdict_params = {
280 '_mvn_doc_default_callparams': _mvn_doc_default_callparams,
281 '_mvn_doc_callparams_note': _mvn_doc_callparams_note,
282 '_doc_random_state': _doc_random_state
283}
285mvn_docdict_noparams = {
286 '_mvn_doc_default_callparams': _mvn_doc_frozen_callparams,
287 '_mvn_doc_callparams_note': _mvn_doc_frozen_callparams_note,
288 '_doc_random_state': _doc_random_state
289}
292class multivariate_normal_gen(multi_rv_generic):
293 r"""A multivariate normal random variable.
295 The `mean` keyword specifies the mean. The `cov` keyword specifies the
296 covariance matrix.
298 Methods
299 -------
300 pdf(x, mean=None, cov=1, allow_singular=False)
301 Probability density function.
302 logpdf(x, mean=None, cov=1, allow_singular=False)
303 Log of the probability density function.
304 cdf(x, mean=None, cov=1, allow_singular=False, maxpts=1000000*dim, abseps=1e-5, releps=1e-5, lower_limit=None) # noqa
305 Cumulative distribution function.
306 logcdf(x, mean=None, cov=1, allow_singular=False, maxpts=1000000*dim, abseps=1e-5, releps=1e-5)
307 Log of the cumulative distribution function.
308 rvs(mean=None, cov=1, size=1, random_state=None)
309 Draw random samples from a multivariate normal distribution.
310 entropy()
311 Compute the differential entropy of the multivariate normal.
313 Parameters
314 ----------
315 %(_mvn_doc_default_callparams)s
316 %(_doc_random_state)s
318 Notes
319 -----
320 %(_mvn_doc_callparams_note)s
322 The covariance matrix `cov` may be an instance of a subclass of
323 `Covariance`, e.g. `scipy.stats.CovViaPrecision`. If so, `allow_singular`
324 is ignored.
326 Otherwise, `cov` must be a symmetric positive semidefinite
327 matrix when `allow_singular` is True; it must be (strictly) positive
328 definite when `allow_singular` is False.
329 Symmetry is not checked; only the lower triangular portion is used.
330 The determinant and inverse of `cov` are computed
331 as the pseudo-determinant and pseudo-inverse, respectively, so
332 that `cov` does not need to have full rank.
334 The probability density function for `multivariate_normal` is
336 .. math::
338 f(x) = \frac{1}{\sqrt{(2 \pi)^k \det \Sigma}}
339 \exp\left( -\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right),
341 where :math:`\mu` is the mean, :math:`\Sigma` the covariance matrix,
342 :math:`k` the rank of :math:`\Sigma`. In case of singular :math:`\Sigma`,
343 SciPy extends this definition according to [1]_.
345 .. versionadded:: 0.14.0
347 References
348 ----------
349 .. [1] Multivariate Normal Distribution - Degenerate Case, Wikipedia,
350 https://en.wikipedia.org/wiki/Multivariate_normal_distribution#Degenerate_case
352 Examples
353 --------
354 >>> import numpy as np
355 >>> import matplotlib.pyplot as plt
356 >>> from scipy.stats import multivariate_normal
358 >>> x = np.linspace(0, 5, 10, endpoint=False)
359 >>> y = multivariate_normal.pdf(x, mean=2.5, cov=0.5); y
360 array([ 0.00108914, 0.01033349, 0.05946514, 0.20755375, 0.43939129,
361 0.56418958, 0.43939129, 0.20755375, 0.05946514, 0.01033349])
362 >>> fig1 = plt.figure()
363 >>> ax = fig1.add_subplot(111)
364 >>> ax.plot(x, y)
365 >>> plt.show()
367 Alternatively, the object may be called (as a function) to fix the mean
368 and covariance parameters, returning a "frozen" multivariate normal
369 random variable:
371 >>> rv = multivariate_normal(mean=None, cov=1, allow_singular=False)
372 >>> # Frozen object with the same methods but holding the given
373 >>> # mean and covariance fixed.
375 The input quantiles can be any shape of array, as long as the last
376 axis labels the components. This allows us for instance to
377 display the frozen pdf for a non-isotropic random variable in 2D as
378 follows:
380 >>> x, y = np.mgrid[-1:1:.01, -1:1:.01]
381 >>> pos = np.dstack((x, y))
382 >>> rv = multivariate_normal([0.5, -0.2], [[2.0, 0.3], [0.3, 0.5]])
383 >>> fig2 = plt.figure()
384 >>> ax2 = fig2.add_subplot(111)
385 >>> ax2.contourf(x, y, rv.pdf(pos))
387 """
389 def __init__(self, seed=None):
390 super().__init__(seed)
391 self.__doc__ = doccer.docformat(self.__doc__, mvn_docdict_params)
393 def __call__(self, mean=None, cov=1, allow_singular=False, seed=None):
394 """Create a frozen multivariate normal distribution.
396 See `multivariate_normal_frozen` for more information.
397 """
398 return multivariate_normal_frozen(mean, cov,
399 allow_singular=allow_singular,
400 seed=seed)
402 def _process_parameters(self, mean, cov, allow_singular=True):
403 """
404 Infer dimensionality from mean or covariance matrix, ensure that
405 mean and covariance are full vector resp. matrix.
406 """
407 if isinstance(cov, _covariance.Covariance):
408 return self._process_parameters_Covariance(mean, cov)
409 else:
410 # Before `Covariance` classes were introduced,
411 # `multivariate_normal` accepted plain arrays as `cov` and used the
412 # following input validation. To avoid disturbing the behavior of
413 # `multivariate_normal` when plain arrays are used, we use the
414 # original input validation here.
415 dim, mean, cov = self._process_parameters_psd(None, mean, cov)
416 # After input validation, some methods then processed the arrays
417 # with a `_PSD` object and used that to perform computation.
418 # To avoid branching statements in each method depending on whether
419 # `cov` is an array or `Covariance` object, we always process the
420 # array with `_PSD`, and then use wrapper that satisfies the
421 # `Covariance` interface, `CovViaPSD`.
422 psd = _PSD(cov, allow_singular=allow_singular)
423 cov_object = _covariance.CovViaPSD(psd)
424 return dim, mean, cov_object
426 def _process_parameters_Covariance(self, mean, cov):
427 dim = cov.shape[-1]
428 mean = np.array([0.]) if mean is None else mean
429 message = (f"`cov` represents a covariance matrix in {dim} dimensions,"
430 f"and so `mean` must be broadcastable to shape {(dim,)}")
431 try:
432 mean = np.broadcast_to(mean, dim)
433 except ValueError as e:
434 raise ValueError(message) from e
435 return dim, mean, cov
437 def _process_parameters_psd(self, dim, mean, cov):
438 # Try to infer dimensionality
439 if dim is None:
440 if mean is None:
441 if cov is None:
442 dim = 1
443 else:
444 cov = np.asarray(cov, dtype=float)
445 if cov.ndim < 2:
446 dim = 1
447 else:
448 dim = cov.shape[0]
449 else:
450 mean = np.asarray(mean, dtype=float)
451 dim = mean.size
452 else:
453 if not np.isscalar(dim):
454 raise ValueError("Dimension of random variable must be "
455 "a scalar.")
457 # Check input sizes and return full arrays for mean and cov if
458 # necessary
459 if mean is None:
460 mean = np.zeros(dim)
461 mean = np.asarray(mean, dtype=float)
463 if cov is None:
464 cov = 1.0
465 cov = np.asarray(cov, dtype=float)
467 if dim == 1:
468 mean = mean.reshape(1)
469 cov = cov.reshape(1, 1)
471 if mean.ndim != 1 or mean.shape[0] != dim:
472 raise ValueError("Array 'mean' must be a vector of length %d." %
473 dim)
474 if cov.ndim == 0:
475 cov = cov * np.eye(dim)
476 elif cov.ndim == 1:
477 cov = np.diag(cov)
478 elif cov.ndim == 2 and cov.shape != (dim, dim):
479 rows, cols = cov.shape
480 if rows != cols:
481 msg = ("Array 'cov' must be square if it is two dimensional,"
482 " but cov.shape = %s." % str(cov.shape))
483 else:
484 msg = ("Dimension mismatch: array 'cov' is of shape %s,"
485 " but 'mean' is a vector of length %d.")
486 msg = msg % (str(cov.shape), len(mean))
487 raise ValueError(msg)
488 elif cov.ndim > 2:
489 raise ValueError("Array 'cov' must be at most two-dimensional,"
490 " but cov.ndim = %d" % cov.ndim)
492 return dim, mean, cov
494 def _process_quantiles(self, x, dim):
495 """
496 Adjust quantiles array so that last axis labels the components of
497 each data point.
498 """
499 x = np.asarray(x, dtype=float)
501 if x.ndim == 0:
502 x = x[np.newaxis]
503 elif x.ndim == 1:
504 if dim == 1:
505 x = x[:, np.newaxis]
506 else:
507 x = x[np.newaxis, :]
509 return x
511 def _logpdf(self, x, mean, cov_object):
512 """Log of the multivariate normal probability density function.
514 Parameters
515 ----------
516 x : ndarray
517 Points at which to evaluate the log of the probability
518 density function
519 mean : ndarray
520 Mean of the distribution
521 cov_object : Covariance
522 An object representing the Covariance matrix
524 Notes
525 -----
526 As this function does no argument checking, it should not be
527 called directly; use 'logpdf' instead.
529 """
530 log_det_cov, rank = cov_object.log_pdet, cov_object.rank
531 dev = x - mean
532 if dev.ndim > 1:
533 log_det_cov = log_det_cov[..., np.newaxis]
534 rank = rank[..., np.newaxis]
535 maha = np.sum(np.square(cov_object.whiten(dev)), axis=-1)
536 return -0.5 * (rank * _LOG_2PI + log_det_cov + maha)
538 def logpdf(self, x, mean=None, cov=1, allow_singular=False):
539 """Log of the multivariate normal probability density function.
541 Parameters
542 ----------
543 x : array_like
544 Quantiles, with the last axis of `x` denoting the components.
545 %(_mvn_doc_default_callparams)s
547 Returns
548 -------
549 pdf : ndarray or scalar
550 Log of the probability density function evaluated at `x`
552 Notes
553 -----
554 %(_mvn_doc_callparams_note)s
556 """
557 params = self._process_parameters(mean, cov, allow_singular)
558 dim, mean, cov_object = params
559 x = self._process_quantiles(x, dim)
560 out = self._logpdf(x, mean, cov_object)
561 if np.any(cov_object.rank < dim):
562 out_of_bounds = ~cov_object._support_mask(x-mean)
563 out[out_of_bounds] = -np.inf
564 return _squeeze_output(out)
566 def pdf(self, x, mean=None, cov=1, allow_singular=False):
567 """Multivariate normal probability density function.
569 Parameters
570 ----------
571 x : array_like
572 Quantiles, with the last axis of `x` denoting the components.
573 %(_mvn_doc_default_callparams)s
575 Returns
576 -------
577 pdf : ndarray or scalar
578 Probability density function evaluated at `x`
580 Notes
581 -----
582 %(_mvn_doc_callparams_note)s
584 """
585 params = self._process_parameters(mean, cov, allow_singular)
586 dim, mean, cov_object = params
587 x = self._process_quantiles(x, dim)
588 out = np.exp(self._logpdf(x, mean, cov_object))
589 if np.any(cov_object.rank < dim):
590 out_of_bounds = ~cov_object._support_mask(x-mean)
591 out[out_of_bounds] = 0.0
592 return _squeeze_output(out)
594 def _cdf(self, x, mean, cov, maxpts, abseps, releps, lower_limit):
595 """Multivariate normal cumulative distribution function.
597 Parameters
598 ----------
599 x : ndarray
600 Points at which to evaluate the cumulative distribution function.
601 mean : ndarray
602 Mean of the distribution
603 cov : array_like
604 Covariance matrix of the distribution
605 maxpts : integer
606 The maximum number of points to use for integration
607 abseps : float
608 Absolute error tolerance
609 releps : float
610 Relative error tolerance
611 lower_limit : array_like, optional
612 Lower limit of integration of the cumulative distribution function.
613 Default is negative infinity. Must be broadcastable with `x`.
615 Notes
616 -----
617 As this function does no argument checking, it should not be
618 called directly; use 'cdf' instead.
621 .. versionadded:: 1.0.0
623 """
624 lower = (np.full(mean.shape, -np.inf)
625 if lower_limit is None else lower_limit)
626 # In 2d, _mvn.mvnun accepts input in which `lower` bound elements
627 # are greater than `x`. Not so in other dimensions. Fix this by
628 # ensuring that lower bounds are indeed lower when passed, then
629 # set signs of resulting CDF manually.
630 b, a = np.broadcast_arrays(x, lower)
631 i_swap = b < a
632 signs = (-1)**(i_swap.sum(axis=-1)) # odd # of swaps -> negative
633 a, b = a.copy(), b.copy()
634 a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
635 n = x.shape[-1]
636 limits = np.concatenate((a, b), axis=-1)
638 # mvnun expects 1-d arguments, so process points sequentially
639 def func1d(limits):
640 return _mvn.mvnun(limits[:n], limits[n:], mean, cov,
641 maxpts, abseps, releps)[0]
643 out = np.apply_along_axis(func1d, -1, limits) * signs
644 return _squeeze_output(out)
646 def logcdf(self, x, mean=None, cov=1, allow_singular=False, maxpts=None,
647 abseps=1e-5, releps=1e-5, *, lower_limit=None):
648 """Log of the multivariate normal cumulative distribution function.
650 Parameters
651 ----------
652 x : array_like
653 Quantiles, with the last axis of `x` denoting the components.
654 %(_mvn_doc_default_callparams)s
655 maxpts : integer, optional
656 The maximum number of points to use for integration
657 (default `1000000*dim`)
658 abseps : float, optional
659 Absolute error tolerance (default 1e-5)
660 releps : float, optional
661 Relative error tolerance (default 1e-5)
662 lower_limit : array_like, optional
663 Lower limit of integration of the cumulative distribution function.
664 Default is negative infinity. Must be broadcastable with `x`.
666 Returns
667 -------
668 cdf : ndarray or scalar
669 Log of the cumulative distribution function evaluated at `x`
671 Notes
672 -----
673 %(_mvn_doc_callparams_note)s
675 .. versionadded:: 1.0.0
677 """
678 params = self._process_parameters(mean, cov, allow_singular)
679 dim, mean, cov_object = params
680 cov = cov_object.covariance
681 x = self._process_quantiles(x, dim)
682 if not maxpts:
683 maxpts = 1000000 * dim
684 cdf = self._cdf(x, mean, cov, maxpts, abseps, releps, lower_limit)
685 # the log of a negative real is complex, and cdf can be negative
686 # if lower limit is greater than upper limit
687 cdf = cdf + 0j if np.any(cdf < 0) else cdf
688 out = np.log(cdf)
689 return out
691 def cdf(self, x, mean=None, cov=1, allow_singular=False, maxpts=None,
692 abseps=1e-5, releps=1e-5, *, lower_limit=None):
693 """Multivariate normal cumulative distribution function.
695 Parameters
696 ----------
697 x : array_like
698 Quantiles, with the last axis of `x` denoting the components.
699 %(_mvn_doc_default_callparams)s
700 maxpts : integer, optional
701 The maximum number of points to use for integration
702 (default `1000000*dim`)
703 abseps : float, optional
704 Absolute error tolerance (default 1e-5)
705 releps : float, optional
706 Relative error tolerance (default 1e-5)
707 lower_limit : array_like, optional
708 Lower limit of integration of the cumulative distribution function.
709 Default is negative infinity. Must be broadcastable with `x`.
711 Returns
712 -------
713 cdf : ndarray or scalar
714 Cumulative distribution function evaluated at `x`
716 Notes
717 -----
718 %(_mvn_doc_callparams_note)s
720 .. versionadded:: 1.0.0
722 """
723 params = self._process_parameters(mean, cov, allow_singular)
724 dim, mean, cov_object = params
725 cov = cov_object.covariance
726 x = self._process_quantiles(x, dim)
727 if not maxpts:
728 maxpts = 1000000 * dim
729 out = self._cdf(x, mean, cov, maxpts, abseps, releps, lower_limit)
730 return out
732 def rvs(self, mean=None, cov=1, size=1, random_state=None):
733 """Draw random samples from a multivariate normal distribution.
735 Parameters
736 ----------
737 %(_mvn_doc_default_callparams)s
738 size : integer, optional
739 Number of samples to draw (default 1).
740 %(_doc_random_state)s
742 Returns
743 -------
744 rvs : ndarray or scalar
745 Random variates of size (`size`, `N`), where `N` is the
746 dimension of the random variable.
748 Notes
749 -----
750 %(_mvn_doc_callparams_note)s
752 """
753 dim, mean, cov_object = self._process_parameters(mean, cov)
754 random_state = self._get_random_state(random_state)
756 if isinstance(cov_object, _covariance.CovViaPSD):
757 cov = cov_object.covariance
758 out = random_state.multivariate_normal(mean, cov, size)
759 out = _squeeze_output(out)
760 else:
761 size = size or tuple()
762 if not np.iterable(size):
763 size = (size,)
764 shape = tuple(size) + (cov_object.shape[-1],)
765 x = random_state.normal(size=shape)
766 out = mean + cov_object.colorize(x)
767 return out
769 def entropy(self, mean=None, cov=1):
770 """Compute the differential entropy of the multivariate normal.
772 Parameters
773 ----------
774 %(_mvn_doc_default_callparams)s
776 Returns
777 -------
778 h : scalar
779 Entropy of the multivariate normal distribution
781 Notes
782 -----
783 %(_mvn_doc_callparams_note)s
785 """
786 dim, mean, cov_object = self._process_parameters(mean, cov)
787 return 0.5 * (cov_object.rank * (_LOG_2PI + 1) + cov_object.log_pdet)
790multivariate_normal = multivariate_normal_gen()
793class multivariate_normal_frozen(multi_rv_frozen):
794 def __init__(self, mean=None, cov=1, allow_singular=False, seed=None,
795 maxpts=None, abseps=1e-5, releps=1e-5):
796 """Create a frozen multivariate normal distribution.
798 Parameters
799 ----------
800 mean : array_like, default: ``[0]``
801 Mean of the distribution.
802 cov : array_like, default: ``[1]``
803 Symmetric positive (semi)definite covariance matrix of the
804 distribution.
805 allow_singular : bool, default: ``False``
806 Whether to allow a singular covariance matrix.
807 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
808 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
809 singleton is used.
810 If `seed` is an int, a new ``RandomState`` instance is used,
811 seeded with `seed`.
812 If `seed` is already a ``Generator`` or ``RandomState`` instance
813 then that instance is used.
814 maxpts : integer, optional
815 The maximum number of points to use for integration of the
816 cumulative distribution function (default `1000000*dim`)
817 abseps : float, optional
818 Absolute error tolerance for the cumulative distribution function
819 (default 1e-5)
820 releps : float, optional
821 Relative error tolerance for the cumulative distribution function
822 (default 1e-5)
824 Examples
825 --------
826 When called with the default parameters, this will create a 1D random
827 variable with mean 0 and covariance 1:
829 >>> from scipy.stats import multivariate_normal
830 >>> r = multivariate_normal()
831 >>> r.mean
832 array([ 0.])
833 >>> r.cov
834 array([[1.]])
836 """
837 self._dist = multivariate_normal_gen(seed)
838 self.dim, self.mean, self.cov_object = (
839 self._dist._process_parameters(mean, cov, allow_singular))
840 self.allow_singular = allow_singular or self.cov_object._allow_singular
841 if not maxpts:
842 maxpts = 1000000 * self.dim
843 self.maxpts = maxpts
844 self.abseps = abseps
845 self.releps = releps
847 @property
848 def cov(self):
849 return self.cov_object.covariance
851 def logpdf(self, x):
852 x = self._dist._process_quantiles(x, self.dim)
853 out = self._dist._logpdf(x, self.mean, self.cov_object)
854 if np.any(self.cov_object.rank < self.dim):
855 out_of_bounds = ~self.cov_object._support_mask(x-self.mean)
856 out[out_of_bounds] = -np.inf
857 return _squeeze_output(out)
859 def pdf(self, x):
860 return np.exp(self.logpdf(x))
862 def logcdf(self, x, *, lower_limit=None):
863 cdf = self.cdf(x, lower_limit=lower_limit)
864 # the log of a negative real is complex, and cdf can be negative
865 # if lower limit is greater than upper limit
866 cdf = cdf + 0j if np.any(cdf < 0) else cdf
867 out = np.log(cdf)
868 return out
870 def cdf(self, x, *, lower_limit=None):
871 x = self._dist._process_quantiles(x, self.dim)
872 out = self._dist._cdf(x, self.mean, self.cov_object.covariance,
873 self.maxpts, self.abseps, self.releps,
874 lower_limit)
875 return _squeeze_output(out)
877 def rvs(self, size=1, random_state=None):
878 return self._dist.rvs(self.mean, self.cov_object, size, random_state)
880 def entropy(self):
881 """Computes the differential entropy of the multivariate normal.
883 Returns
884 -------
885 h : scalar
886 Entropy of the multivariate normal distribution
888 """
889 log_pdet = self.cov_object.log_pdet
890 rank = self.cov_object.rank
891 return 0.5 * (rank * (_LOG_2PI + 1) + log_pdet)
894# Set frozen generator docstrings from corresponding docstrings in
895# multivariate_normal_gen and fill in default strings in class docstrings
896for name in ['logpdf', 'pdf', 'logcdf', 'cdf', 'rvs']:
897 method = multivariate_normal_gen.__dict__[name]
898 method_frozen = multivariate_normal_frozen.__dict__[name]
899 method_frozen.__doc__ = doccer.docformat(method.__doc__,
900 mvn_docdict_noparams)
901 method.__doc__ = doccer.docformat(method.__doc__, mvn_docdict_params)
903_matnorm_doc_default_callparams = """\
904mean : array_like, optional
905 Mean of the distribution (default: `None`)
906rowcov : array_like, optional
907 Among-row covariance matrix of the distribution (default: `1`)
908colcov : array_like, optional
909 Among-column covariance matrix of the distribution (default: `1`)
910"""
912_matnorm_doc_callparams_note = """\
913If `mean` is set to `None` then a matrix of zeros is used for the mean.
914The dimensions of this matrix are inferred from the shape of `rowcov` and
915`colcov`, if these are provided, or set to `1` if ambiguous.
917`rowcov` and `colcov` can be two-dimensional array_likes specifying the
918covariance matrices directly. Alternatively, a one-dimensional array will
919be be interpreted as the entries of a diagonal matrix, and a scalar or
920zero-dimensional array will be interpreted as this value times the
921identity matrix.
922"""
924_matnorm_doc_frozen_callparams = ""
926_matnorm_doc_frozen_callparams_note = """\
927See class definition for a detailed description of parameters."""
929matnorm_docdict_params = {
930 '_matnorm_doc_default_callparams': _matnorm_doc_default_callparams,
931 '_matnorm_doc_callparams_note': _matnorm_doc_callparams_note,
932 '_doc_random_state': _doc_random_state
933}
935matnorm_docdict_noparams = {
936 '_matnorm_doc_default_callparams': _matnorm_doc_frozen_callparams,
937 '_matnorm_doc_callparams_note': _matnorm_doc_frozen_callparams_note,
938 '_doc_random_state': _doc_random_state
939}
942class matrix_normal_gen(multi_rv_generic):
943 r"""A matrix normal random variable.
945 The `mean` keyword specifies the mean. The `rowcov` keyword specifies the
946 among-row covariance matrix. The 'colcov' keyword specifies the
947 among-column covariance matrix.
949 Methods
950 -------
951 pdf(X, mean=None, rowcov=1, colcov=1)
952 Probability density function.
953 logpdf(X, mean=None, rowcov=1, colcov=1)
954 Log of the probability density function.
955 rvs(mean=None, rowcov=1, colcov=1, size=1, random_state=None)
956 Draw random samples.
957 entropy(rowcol=1, colcov=1)
958 Differential entropy.
960 Parameters
961 ----------
962 %(_matnorm_doc_default_callparams)s
963 %(_doc_random_state)s
965 Notes
966 -----
967 %(_matnorm_doc_callparams_note)s
969 The covariance matrices specified by `rowcov` and `colcov` must be
970 (symmetric) positive definite. If the samples in `X` are
971 :math:`m \times n`, then `rowcov` must be :math:`m \times m` and
972 `colcov` must be :math:`n \times n`. `mean` must be the same shape as `X`.
974 The probability density function for `matrix_normal` is
976 .. math::
978 f(X) = (2 \pi)^{-\frac{mn}{2}}|U|^{-\frac{n}{2}} |V|^{-\frac{m}{2}}
979 \exp\left( -\frac{1}{2} \mathrm{Tr}\left[ U^{-1} (X-M) V^{-1}
980 (X-M)^T \right] \right),
982 where :math:`M` is the mean, :math:`U` the among-row covariance matrix,
983 :math:`V` the among-column covariance matrix.
985 The `allow_singular` behaviour of the `multivariate_normal`
986 distribution is not currently supported. Covariance matrices must be
987 full rank.
989 The `matrix_normal` distribution is closely related to the
990 `multivariate_normal` distribution. Specifically, :math:`\mathrm{Vec}(X)`
991 (the vector formed by concatenating the columns of :math:`X`) has a
992 multivariate normal distribution with mean :math:`\mathrm{Vec}(M)`
993 and covariance :math:`V \otimes U` (where :math:`\otimes` is the Kronecker
994 product). Sampling and pdf evaluation are
995 :math:`\mathcal{O}(m^3 + n^3 + m^2 n + m n^2)` for the matrix normal, but
996 :math:`\mathcal{O}(m^3 n^3)` for the equivalent multivariate normal,
997 making this equivalent form algorithmically inefficient.
999 .. versionadded:: 0.17.0
1001 Examples
1002 --------
1004 >>> import numpy as np
1005 >>> from scipy.stats import matrix_normal
1007 >>> M = np.arange(6).reshape(3,2); M
1008 array([[0, 1],
1009 [2, 3],
1010 [4, 5]])
1011 >>> U = np.diag([1,2,3]); U
1012 array([[1, 0, 0],
1013 [0, 2, 0],
1014 [0, 0, 3]])
1015 >>> V = 0.3*np.identity(2); V
1016 array([[ 0.3, 0. ],
1017 [ 0. , 0.3]])
1018 >>> X = M + 0.1; X
1019 array([[ 0.1, 1.1],
1020 [ 2.1, 3.1],
1021 [ 4.1, 5.1]])
1022 >>> matrix_normal.pdf(X, mean=M, rowcov=U, colcov=V)
1023 0.023410202050005054
1025 >>> # Equivalent multivariate normal
1026 >>> from scipy.stats import multivariate_normal
1027 >>> vectorised_X = X.T.flatten()
1028 >>> equiv_mean = M.T.flatten()
1029 >>> equiv_cov = np.kron(V,U)
1030 >>> multivariate_normal.pdf(vectorised_X, mean=equiv_mean, cov=equiv_cov)
1031 0.023410202050005054
1033 Alternatively, the object may be called (as a function) to fix the mean
1034 and covariance parameters, returning a "frozen" matrix normal
1035 random variable:
1037 >>> rv = matrix_normal(mean=None, rowcov=1, colcov=1)
1038 >>> # Frozen object with the same methods but holding the given
1039 >>> # mean and covariance fixed.
1041 """
1043 def __init__(self, seed=None):
1044 super().__init__(seed)
1045 self.__doc__ = doccer.docformat(self.__doc__, matnorm_docdict_params)
1047 def __call__(self, mean=None, rowcov=1, colcov=1, seed=None):
1048 """Create a frozen matrix normal distribution.
1050 See `matrix_normal_frozen` for more information.
1052 """
1053 return matrix_normal_frozen(mean, rowcov, colcov, seed=seed)
1055 def _process_parameters(self, mean, rowcov, colcov):
1056 """
1057 Infer dimensionality from mean or covariance matrices. Handle
1058 defaults. Ensure compatible dimensions.
1059 """
1061 # Process mean
1062 if mean is not None:
1063 mean = np.asarray(mean, dtype=float)
1064 meanshape = mean.shape
1065 if len(meanshape) != 2:
1066 raise ValueError("Array `mean` must be two dimensional.")
1067 if np.any(meanshape == 0):
1068 raise ValueError("Array `mean` has invalid shape.")
1070 # Process among-row covariance
1071 rowcov = np.asarray(rowcov, dtype=float)
1072 if rowcov.ndim == 0:
1073 if mean is not None:
1074 rowcov = rowcov * np.identity(meanshape[0])
1075 else:
1076 rowcov = rowcov * np.identity(1)
1077 elif rowcov.ndim == 1:
1078 rowcov = np.diag(rowcov)
1079 rowshape = rowcov.shape
1080 if len(rowshape) != 2:
1081 raise ValueError("`rowcov` must be a scalar or a 2D array.")
1082 if rowshape[0] != rowshape[1]:
1083 raise ValueError("Array `rowcov` must be square.")
1084 if rowshape[0] == 0:
1085 raise ValueError("Array `rowcov` has invalid shape.")
1086 numrows = rowshape[0]
1088 # Process among-column covariance
1089 colcov = np.asarray(colcov, dtype=float)
1090 if colcov.ndim == 0:
1091 if mean is not None:
1092 colcov = colcov * np.identity(meanshape[1])
1093 else:
1094 colcov = colcov * np.identity(1)
1095 elif colcov.ndim == 1:
1096 colcov = np.diag(colcov)
1097 colshape = colcov.shape
1098 if len(colshape) != 2:
1099 raise ValueError("`colcov` must be a scalar or a 2D array.")
1100 if colshape[0] != colshape[1]:
1101 raise ValueError("Array `colcov` must be square.")
1102 if colshape[0] == 0:
1103 raise ValueError("Array `colcov` has invalid shape.")
1104 numcols = colshape[0]
1106 # Ensure mean and covariances compatible
1107 if mean is not None:
1108 if meanshape[0] != numrows:
1109 raise ValueError("Arrays `mean` and `rowcov` must have the "
1110 "same number of rows.")
1111 if meanshape[1] != numcols:
1112 raise ValueError("Arrays `mean` and `colcov` must have the "
1113 "same number of columns.")
1114 else:
1115 mean = np.zeros((numrows, numcols))
1117 dims = (numrows, numcols)
1119 return dims, mean, rowcov, colcov
1121 def _process_quantiles(self, X, dims):
1122 """
1123 Adjust quantiles array so that last two axes labels the components of
1124 each data point.
1125 """
1126 X = np.asarray(X, dtype=float)
1127 if X.ndim == 2:
1128 X = X[np.newaxis, :]
1129 if X.shape[-2:] != dims:
1130 raise ValueError("The shape of array `X` is not compatible "
1131 "with the distribution parameters.")
1132 return X
1134 def _logpdf(self, dims, X, mean, row_prec_rt, log_det_rowcov,
1135 col_prec_rt, log_det_colcov):
1136 """Log of the matrix normal probability density function.
1138 Parameters
1139 ----------
1140 dims : tuple
1141 Dimensions of the matrix variates
1142 X : ndarray
1143 Points at which to evaluate the log of the probability
1144 density function
1145 mean : ndarray
1146 Mean of the distribution
1147 row_prec_rt : ndarray
1148 A decomposition such that np.dot(row_prec_rt, row_prec_rt.T)
1149 is the inverse of the among-row covariance matrix
1150 log_det_rowcov : float
1151 Logarithm of the determinant of the among-row covariance matrix
1152 col_prec_rt : ndarray
1153 A decomposition such that np.dot(col_prec_rt, col_prec_rt.T)
1154 is the inverse of the among-column covariance matrix
1155 log_det_colcov : float
1156 Logarithm of the determinant of the among-column covariance matrix
1158 Notes
1159 -----
1160 As this function does no argument checking, it should not be
1161 called directly; use 'logpdf' instead.
1163 """
1164 numrows, numcols = dims
1165 roll_dev = np.moveaxis(X-mean, -1, 0)
1166 scale_dev = np.tensordot(col_prec_rt.T,
1167 np.dot(roll_dev, row_prec_rt), 1)
1168 maha = np.sum(np.sum(np.square(scale_dev), axis=-1), axis=0)
1169 return -0.5 * (numrows*numcols*_LOG_2PI + numcols*log_det_rowcov
1170 + numrows*log_det_colcov + maha)
1172 def logpdf(self, X, mean=None, rowcov=1, colcov=1):
1173 """Log of the matrix normal probability density function.
1175 Parameters
1176 ----------
1177 X : array_like
1178 Quantiles, with the last two axes of `X` denoting the components.
1179 %(_matnorm_doc_default_callparams)s
1181 Returns
1182 -------
1183 logpdf : ndarray
1184 Log of the probability density function evaluated at `X`
1186 Notes
1187 -----
1188 %(_matnorm_doc_callparams_note)s
1190 """
1191 dims, mean, rowcov, colcov = self._process_parameters(mean, rowcov,
1192 colcov)
1193 X = self._process_quantiles(X, dims)
1194 rowpsd = _PSD(rowcov, allow_singular=False)
1195 colpsd = _PSD(colcov, allow_singular=False)
1196 out = self._logpdf(dims, X, mean, rowpsd.U, rowpsd.log_pdet, colpsd.U,
1197 colpsd.log_pdet)
1198 return _squeeze_output(out)
1200 def pdf(self, X, mean=None, rowcov=1, colcov=1):
1201 """Matrix normal probability density function.
1203 Parameters
1204 ----------
1205 X : array_like
1206 Quantiles, with the last two axes of `X` denoting the components.
1207 %(_matnorm_doc_default_callparams)s
1209 Returns
1210 -------
1211 pdf : ndarray
1212 Probability density function evaluated at `X`
1214 Notes
1215 -----
1216 %(_matnorm_doc_callparams_note)s
1218 """
1219 return np.exp(self.logpdf(X, mean, rowcov, colcov))
1221 def rvs(self, mean=None, rowcov=1, colcov=1, size=1, random_state=None):
1222 """Draw random samples from a matrix normal distribution.
1224 Parameters
1225 ----------
1226 %(_matnorm_doc_default_callparams)s
1227 size : integer, optional
1228 Number of samples to draw (default 1).
1229 %(_doc_random_state)s
1231 Returns
1232 -------
1233 rvs : ndarray or scalar
1234 Random variates of size (`size`, `dims`), where `dims` is the
1235 dimension of the random matrices.
1237 Notes
1238 -----
1239 %(_matnorm_doc_callparams_note)s
1241 """
1242 size = int(size)
1243 dims, mean, rowcov, colcov = self._process_parameters(mean, rowcov,
1244 colcov)
1245 rowchol = scipy.linalg.cholesky(rowcov, lower=True)
1246 colchol = scipy.linalg.cholesky(colcov, lower=True)
1247 random_state = self._get_random_state(random_state)
1248 # We aren't generating standard normal variates with size=(size,
1249 # dims[0], dims[1]) directly to ensure random variates remain backwards
1250 # compatible. See https://github.com/scipy/scipy/pull/12312 for more
1251 # details.
1252 std_norm = random_state.standard_normal(
1253 size=(dims[1], size, dims[0])
1254 ).transpose(1, 2, 0)
1255 out = mean + np.einsum('jp,ipq,kq->ijk',
1256 rowchol, std_norm, colchol,
1257 optimize=True)
1258 if size == 1:
1259 out = out.reshape(mean.shape)
1260 return out
1262 def entropy(self, rowcov=1, colcov=1):
1263 """Log of the matrix normal probability density function.
1265 Parameters
1266 ----------
1267 rowcov : array_like, optional
1268 Among-row covariance matrix of the distribution (default: `1`)
1269 colcov : array_like, optional
1270 Among-column covariance matrix of the distribution (default: `1`)
1272 Returns
1273 -------
1274 entropy : float
1275 Entropy of the distribution
1277 Notes
1278 -----
1279 %(_matnorm_doc_callparams_note)s
1281 """
1282 dummy_mean = np.zeros((rowcov.shape[0], colcov.shape[0]))
1283 dims, _, rowcov, colcov = self._process_parameters(dummy_mean,
1284 rowcov,
1285 colcov)
1286 rowpsd = _PSD(rowcov, allow_singular=False)
1287 colpsd = _PSD(colcov, allow_singular=False)
1289 return self._entropy(dims, rowpsd.log_pdet, colpsd.log_pdet)
1291 def _entropy(self, dims, row_cov_logdet, col_cov_logdet):
1292 n, p = dims
1293 return (0.5 * n * p * (1 + _LOG_2PI) + 0.5 * p * row_cov_logdet +
1294 0.5 * n * col_cov_logdet)
1297matrix_normal = matrix_normal_gen()
1300class matrix_normal_frozen(multi_rv_frozen):
1301 """
1302 Create a frozen matrix normal distribution.
1304 Parameters
1305 ----------
1306 %(_matnorm_doc_default_callparams)s
1307 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
1308 If `seed` is `None` the `~np.random.RandomState` singleton is used.
1309 If `seed` is an int, a new ``RandomState`` instance is used, seeded
1310 with seed.
1311 If `seed` is already a ``RandomState`` or ``Generator`` instance,
1312 then that object is used.
1313 Default is `None`.
1315 Examples
1316 --------
1317 >>> import numpy as np
1318 >>> from scipy.stats import matrix_normal
1320 >>> distn = matrix_normal(mean=np.zeros((3,3)))
1321 >>> X = distn.rvs(); X
1322 array([[-0.02976962, 0.93339138, -0.09663178],
1323 [ 0.67405524, 0.28250467, -0.93308929],
1324 [-0.31144782, 0.74535536, 1.30412916]])
1325 >>> distn.pdf(X)
1326 2.5160642368346784e-05
1327 >>> distn.logpdf(X)
1328 -10.590229595124615
1329 """
1331 def __init__(self, mean=None, rowcov=1, colcov=1, seed=None):
1332 self._dist = matrix_normal_gen(seed)
1333 self.dims, self.mean, self.rowcov, self.colcov = \
1334 self._dist._process_parameters(mean, rowcov, colcov)
1335 self.rowpsd = _PSD(self.rowcov, allow_singular=False)
1336 self.colpsd = _PSD(self.colcov, allow_singular=False)
1338 def logpdf(self, X):
1339 X = self._dist._process_quantiles(X, self.dims)
1340 out = self._dist._logpdf(self.dims, X, self.mean, self.rowpsd.U,
1341 self.rowpsd.log_pdet, self.colpsd.U,
1342 self.colpsd.log_pdet)
1343 return _squeeze_output(out)
1345 def pdf(self, X):
1346 return np.exp(self.logpdf(X))
1348 def rvs(self, size=1, random_state=None):
1349 return self._dist.rvs(self.mean, self.rowcov, self.colcov, size,
1350 random_state)
1352 def entropy(self):
1353 return self._dist._entropy(self.dims, self.rowpsd.log_pdet,
1354 self.colpsd.log_pdet)
1357# Set frozen generator docstrings from corresponding docstrings in
1358# matrix_normal_gen and fill in default strings in class docstrings
1359for name in ['logpdf', 'pdf', 'rvs', 'entropy']:
1360 method = matrix_normal_gen.__dict__[name]
1361 method_frozen = matrix_normal_frozen.__dict__[name]
1362 method_frozen.__doc__ = doccer.docformat(method.__doc__,
1363 matnorm_docdict_noparams)
1364 method.__doc__ = doccer.docformat(method.__doc__, matnorm_docdict_params)
1366_dirichlet_doc_default_callparams = """\
1367alpha : array_like
1368 The concentration parameters. The number of entries determines the
1369 dimensionality of the distribution.
1370"""
1371_dirichlet_doc_frozen_callparams = ""
1373_dirichlet_doc_frozen_callparams_note = """\
1374See class definition for a detailed description of parameters."""
1376dirichlet_docdict_params = {
1377 '_dirichlet_doc_default_callparams': _dirichlet_doc_default_callparams,
1378 '_doc_random_state': _doc_random_state
1379}
1381dirichlet_docdict_noparams = {
1382 '_dirichlet_doc_default_callparams': _dirichlet_doc_frozen_callparams,
1383 '_doc_random_state': _doc_random_state
1384}
1387def _dirichlet_check_parameters(alpha):
1388 alpha = np.asarray(alpha)
1389 if np.min(alpha) <= 0:
1390 raise ValueError("All parameters must be greater than 0")
1391 elif alpha.ndim != 1:
1392 raise ValueError("Parameter vector 'a' must be one dimensional, "
1393 "but a.shape = {}.".format(alpha.shape))
1394 return alpha
1397def _dirichlet_check_input(alpha, x):
1398 x = np.asarray(x)
1400 if x.shape[0] + 1 != alpha.shape[0] and x.shape[0] != alpha.shape[0]:
1401 raise ValueError("Vector 'x' must have either the same number "
1402 "of entries as, or one entry fewer than, "
1403 "parameter vector 'a', but alpha.shape = {} "
1404 "and x.shape = {}.".format(alpha.shape, x.shape))
1406 if x.shape[0] != alpha.shape[0]:
1407 xk = np.array([1 - np.sum(x, 0)])
1408 if xk.ndim == 1:
1409 x = np.append(x, xk)
1410 elif xk.ndim == 2:
1411 x = np.vstack((x, xk))
1412 else:
1413 raise ValueError("The input must be one dimensional or a two "
1414 "dimensional matrix containing the entries.")
1416 if np.min(x) < 0:
1417 raise ValueError("Each entry in 'x' must be greater than or equal "
1418 "to zero.")
1420 if np.max(x) > 1:
1421 raise ValueError("Each entry in 'x' must be smaller or equal one.")
1423 # Check x_i > 0 or alpha_i > 1
1424 xeq0 = (x == 0)
1425 alphalt1 = (alpha < 1)
1426 if x.shape != alpha.shape:
1427 alphalt1 = np.repeat(alphalt1, x.shape[-1], axis=-1).reshape(x.shape)
1428 chk = np.logical_and(xeq0, alphalt1)
1430 if np.sum(chk):
1431 raise ValueError("Each entry in 'x' must be greater than zero if its "
1432 "alpha is less than one.")
1434 if (np.abs(np.sum(x, 0) - 1.0) > 10e-10).any():
1435 raise ValueError("The input vector 'x' must lie within the normal "
1436 "simplex. but np.sum(x, 0) = %s." % np.sum(x, 0))
1438 return x
1441def _lnB(alpha):
1442 r"""Internal helper function to compute the log of the useful quotient.
1444 .. math::
1446 B(\alpha) = \frac{\prod_{i=1}{K}\Gamma(\alpha_i)}
1447 {\Gamma\left(\sum_{i=1}^{K} \alpha_i \right)}
1449 Parameters
1450 ----------
1451 %(_dirichlet_doc_default_callparams)s
1453 Returns
1454 -------
1455 B : scalar
1456 Helper quotient, internal use only
1458 """
1459 return np.sum(gammaln(alpha)) - gammaln(np.sum(alpha))
1462class dirichlet_gen(multi_rv_generic):
1463 r"""A Dirichlet random variable.
1465 The ``alpha`` keyword specifies the concentration parameters of the
1466 distribution.
1468 .. versionadded:: 0.15.0
1470 Methods
1471 -------
1472 pdf(x, alpha)
1473 Probability density function.
1474 logpdf(x, alpha)
1475 Log of the probability density function.
1476 rvs(alpha, size=1, random_state=None)
1477 Draw random samples from a Dirichlet distribution.
1478 mean(alpha)
1479 The mean of the Dirichlet distribution
1480 var(alpha)
1481 The variance of the Dirichlet distribution
1482 entropy(alpha)
1483 Compute the differential entropy of the Dirichlet distribution.
1485 Parameters
1486 ----------
1487 %(_dirichlet_doc_default_callparams)s
1488 %(_doc_random_state)s
1490 Notes
1491 -----
1492 Each :math:`\alpha` entry must be positive. The distribution has only
1493 support on the simplex defined by
1495 .. math::
1496 \sum_{i=1}^{K} x_i = 1
1498 where :math:`0 < x_i < 1`.
1500 If the quantiles don't lie within the simplex, a ValueError is raised.
1502 The probability density function for `dirichlet` is
1504 .. math::
1506 f(x) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^K x_i^{\alpha_i - 1}
1508 where
1510 .. math::
1512 \mathrm{B}(\boldsymbol\alpha) = \frac{\prod_{i=1}^K \Gamma(\alpha_i)}
1513 {\Gamma\bigl(\sum_{i=1}^K \alpha_i\bigr)}
1515 and :math:`\boldsymbol\alpha=(\alpha_1,\ldots,\alpha_K)`, the
1516 concentration parameters and :math:`K` is the dimension of the space
1517 where :math:`x` takes values.
1519 Note that the `dirichlet` interface is somewhat inconsistent.
1520 The array returned by the rvs function is transposed
1521 with respect to the format expected by the pdf and logpdf.
1523 Examples
1524 --------
1525 >>> import numpy as np
1526 >>> from scipy.stats import dirichlet
1528 Generate a dirichlet random variable
1530 >>> quantiles = np.array([0.2, 0.2, 0.6]) # specify quantiles
1531 >>> alpha = np.array([0.4, 5, 15]) # specify concentration parameters
1532 >>> dirichlet.pdf(quantiles, alpha)
1533 0.2843831684937255
1535 The same PDF but following a log scale
1537 >>> dirichlet.logpdf(quantiles, alpha)
1538 -1.2574327653159187
1540 Once we specify the dirichlet distribution
1541 we can then calculate quantities of interest
1543 >>> dirichlet.mean(alpha) # get the mean of the distribution
1544 array([0.01960784, 0.24509804, 0.73529412])
1545 >>> dirichlet.var(alpha) # get variance
1546 array([0.00089829, 0.00864603, 0.00909517])
1547 >>> dirichlet.entropy(alpha) # calculate the differential entropy
1548 -4.3280162474082715
1550 We can also return random samples from the distribution
1552 >>> dirichlet.rvs(alpha, size=1, random_state=1)
1553 array([[0.00766178, 0.24670518, 0.74563305]])
1554 >>> dirichlet.rvs(alpha, size=2, random_state=2)
1555 array([[0.01639427, 0.1292273 , 0.85437844],
1556 [0.00156917, 0.19033695, 0.80809388]])
1558 Alternatively, the object may be called (as a function) to fix
1559 concentration parameters, returning a "frozen" Dirichlet
1560 random variable:
1562 >>> rv = dirichlet(alpha)
1563 >>> # Frozen object with the same methods but holding the given
1564 >>> # concentration parameters fixed.
1566 """
1568 def __init__(self, seed=None):
1569 super().__init__(seed)
1570 self.__doc__ = doccer.docformat(self.__doc__, dirichlet_docdict_params)
1572 def __call__(self, alpha, seed=None):
1573 return dirichlet_frozen(alpha, seed=seed)
1575 def _logpdf(self, x, alpha):
1576 """Log of the Dirichlet probability density function.
1578 Parameters
1579 ----------
1580 x : ndarray
1581 Points at which to evaluate the log of the probability
1582 density function
1583 %(_dirichlet_doc_default_callparams)s
1585 Notes
1586 -----
1587 As this function does no argument checking, it should not be
1588 called directly; use 'logpdf' instead.
1590 """
1591 lnB = _lnB(alpha)
1592 return - lnB + np.sum((xlogy(alpha - 1, x.T)).T, 0)
1594 def logpdf(self, x, alpha):
1595 """Log of the Dirichlet probability density function.
1597 Parameters
1598 ----------
1599 x : array_like
1600 Quantiles, with the last axis of `x` denoting the components.
1601 %(_dirichlet_doc_default_callparams)s
1603 Returns
1604 -------
1605 pdf : ndarray or scalar
1606 Log of the probability density function evaluated at `x`.
1608 """
1609 alpha = _dirichlet_check_parameters(alpha)
1610 x = _dirichlet_check_input(alpha, x)
1612 out = self._logpdf(x, alpha)
1613 return _squeeze_output(out)
1615 def pdf(self, x, alpha):
1616 """The Dirichlet probability density function.
1618 Parameters
1619 ----------
1620 x : array_like
1621 Quantiles, with the last axis of `x` denoting the components.
1622 %(_dirichlet_doc_default_callparams)s
1624 Returns
1625 -------
1626 pdf : ndarray or scalar
1627 The probability density function evaluated at `x`.
1629 """
1630 alpha = _dirichlet_check_parameters(alpha)
1631 x = _dirichlet_check_input(alpha, x)
1633 out = np.exp(self._logpdf(x, alpha))
1634 return _squeeze_output(out)
1636 def mean(self, alpha):
1637 """Mean of the Dirichlet distribution.
1639 Parameters
1640 ----------
1641 %(_dirichlet_doc_default_callparams)s
1643 Returns
1644 -------
1645 mu : ndarray or scalar
1646 Mean of the Dirichlet distribution.
1648 """
1649 alpha = _dirichlet_check_parameters(alpha)
1651 out = alpha / (np.sum(alpha))
1652 return _squeeze_output(out)
1654 def var(self, alpha):
1655 """Variance of the Dirichlet distribution.
1657 Parameters
1658 ----------
1659 %(_dirichlet_doc_default_callparams)s
1661 Returns
1662 -------
1663 v : ndarray or scalar
1664 Variance of the Dirichlet distribution.
1666 """
1668 alpha = _dirichlet_check_parameters(alpha)
1670 alpha0 = np.sum(alpha)
1671 out = (alpha * (alpha0 - alpha)) / ((alpha0 * alpha0) * (alpha0 + 1))
1672 return _squeeze_output(out)
1674 def entropy(self, alpha):
1675 """
1676 Differential entropy of the Dirichlet distribution.
1678 Parameters
1679 ----------
1680 %(_dirichlet_doc_default_callparams)s
1682 Returns
1683 -------
1684 h : scalar
1685 Entropy of the Dirichlet distribution
1687 """
1689 alpha = _dirichlet_check_parameters(alpha)
1691 alpha0 = np.sum(alpha)
1692 lnB = _lnB(alpha)
1693 K = alpha.shape[0]
1695 out = lnB + (alpha0 - K) * scipy.special.psi(alpha0) - np.sum(
1696 (alpha - 1) * scipy.special.psi(alpha))
1697 return _squeeze_output(out)
1699 def rvs(self, alpha, size=1, random_state=None):
1700 """
1701 Draw random samples from a Dirichlet distribution.
1703 Parameters
1704 ----------
1705 %(_dirichlet_doc_default_callparams)s
1706 size : int, optional
1707 Number of samples to draw (default 1).
1708 %(_doc_random_state)s
1710 Returns
1711 -------
1712 rvs : ndarray or scalar
1713 Random variates of size (`size`, `N`), where `N` is the
1714 dimension of the random variable.
1716 """
1717 alpha = _dirichlet_check_parameters(alpha)
1718 random_state = self._get_random_state(random_state)
1719 return random_state.dirichlet(alpha, size=size)
1722dirichlet = dirichlet_gen()
1725class dirichlet_frozen(multi_rv_frozen):
1726 def __init__(self, alpha, seed=None):
1727 self.alpha = _dirichlet_check_parameters(alpha)
1728 self._dist = dirichlet_gen(seed)
1730 def logpdf(self, x):
1731 return self._dist.logpdf(x, self.alpha)
1733 def pdf(self, x):
1734 return self._dist.pdf(x, self.alpha)
1736 def mean(self):
1737 return self._dist.mean(self.alpha)
1739 def var(self):
1740 return self._dist.var(self.alpha)
1742 def entropy(self):
1743 return self._dist.entropy(self.alpha)
1745 def rvs(self, size=1, random_state=None):
1746 return self._dist.rvs(self.alpha, size, random_state)
1749# Set frozen generator docstrings from corresponding docstrings in
1750# multivariate_normal_gen and fill in default strings in class docstrings
1751for name in ['logpdf', 'pdf', 'rvs', 'mean', 'var', 'entropy']:
1752 method = dirichlet_gen.__dict__[name]
1753 method_frozen = dirichlet_frozen.__dict__[name]
1754 method_frozen.__doc__ = doccer.docformat(
1755 method.__doc__, dirichlet_docdict_noparams)
1756 method.__doc__ = doccer.docformat(method.__doc__, dirichlet_docdict_params)
1759_wishart_doc_default_callparams = """\
1760df : int
1761 Degrees of freedom, must be greater than or equal to dimension of the
1762 scale matrix
1763scale : array_like
1764 Symmetric positive definite scale matrix of the distribution
1765"""
1767_wishart_doc_callparams_note = ""
1769_wishart_doc_frozen_callparams = ""
1771_wishart_doc_frozen_callparams_note = """\
1772See class definition for a detailed description of parameters."""
1774wishart_docdict_params = {
1775 '_doc_default_callparams': _wishart_doc_default_callparams,
1776 '_doc_callparams_note': _wishart_doc_callparams_note,
1777 '_doc_random_state': _doc_random_state
1778}
1780wishart_docdict_noparams = {
1781 '_doc_default_callparams': _wishart_doc_frozen_callparams,
1782 '_doc_callparams_note': _wishart_doc_frozen_callparams_note,
1783 '_doc_random_state': _doc_random_state
1784}
1787class wishart_gen(multi_rv_generic):
1788 r"""A Wishart random variable.
1790 The `df` keyword specifies the degrees of freedom. The `scale` keyword
1791 specifies the scale matrix, which must be symmetric and positive definite.
1792 In this context, the scale matrix is often interpreted in terms of a
1793 multivariate normal precision matrix (the inverse of the covariance
1794 matrix). These arguments must satisfy the relationship
1795 ``df > scale.ndim - 1``, but see notes on using the `rvs` method with
1796 ``df < scale.ndim``.
1798 Methods
1799 -------
1800 pdf(x, df, scale)
1801 Probability density function.
1802 logpdf(x, df, scale)
1803 Log of the probability density function.
1804 rvs(df, scale, size=1, random_state=None)
1805 Draw random samples from a Wishart distribution.
1806 entropy()
1807 Compute the differential entropy of the Wishart distribution.
1809 Parameters
1810 ----------
1811 %(_doc_default_callparams)s
1812 %(_doc_random_state)s
1814 Raises
1815 ------
1816 scipy.linalg.LinAlgError
1817 If the scale matrix `scale` is not positive definite.
1819 See Also
1820 --------
1821 invwishart, chi2
1823 Notes
1824 -----
1825 %(_doc_callparams_note)s
1827 The scale matrix `scale` must be a symmetric positive definite
1828 matrix. Singular matrices, including the symmetric positive semi-definite
1829 case, are not supported. Symmetry is not checked; only the lower triangular
1830 portion is used.
1832 The Wishart distribution is often denoted
1834 .. math::
1836 W_p(\nu, \Sigma)
1838 where :math:`\nu` is the degrees of freedom and :math:`\Sigma` is the
1839 :math:`p \times p` scale matrix.
1841 The probability density function for `wishart` has support over positive
1842 definite matrices :math:`S`; if :math:`S \sim W_p(\nu, \Sigma)`, then
1843 its PDF is given by:
1845 .. math::
1847 f(S) = \frac{|S|^{\frac{\nu - p - 1}{2}}}{2^{ \frac{\nu p}{2} }
1848 |\Sigma|^\frac{\nu}{2} \Gamma_p \left ( \frac{\nu}{2} \right )}
1849 \exp\left( -tr(\Sigma^{-1} S) / 2 \right)
1851 If :math:`S \sim W_p(\nu, \Sigma)` (Wishart) then
1852 :math:`S^{-1} \sim W_p^{-1}(\nu, \Sigma^{-1})` (inverse Wishart).
1854 If the scale matrix is 1-dimensional and equal to one, then the Wishart
1855 distribution :math:`W_1(\nu, 1)` collapses to the :math:`\chi^2(\nu)`
1856 distribution.
1858 The algorithm [2]_ implemented by the `rvs` method may
1859 produce numerically singular matrices with :math:`p - 1 < \nu < p`; the
1860 user may wish to check for this condition and generate replacement samples
1861 as necessary.
1864 .. versionadded:: 0.16.0
1866 References
1867 ----------
1868 .. [1] M.L. Eaton, "Multivariate Statistics: A Vector Space Approach",
1869 Wiley, 1983.
1870 .. [2] W.B. Smith and R.R. Hocking, "Algorithm AS 53: Wishart Variate
1871 Generator", Applied Statistics, vol. 21, pp. 341-345, 1972.
1873 Examples
1874 --------
1875 >>> import numpy as np
1876 >>> import matplotlib.pyplot as plt
1877 >>> from scipy.stats import wishart, chi2
1878 >>> x = np.linspace(1e-5, 8, 100)
1879 >>> w = wishart.pdf(x, df=3, scale=1); w[:5]
1880 array([ 0.00126156, 0.10892176, 0.14793434, 0.17400548, 0.1929669 ])
1881 >>> c = chi2.pdf(x, 3); c[:5]
1882 array([ 0.00126156, 0.10892176, 0.14793434, 0.17400548, 0.1929669 ])
1883 >>> plt.plot(x, w)
1884 >>> plt.show()
1886 The input quantiles can be any shape of array, as long as the last
1887 axis labels the components.
1889 Alternatively, the object may be called (as a function) to fix the degrees
1890 of freedom and scale parameters, returning a "frozen" Wishart random
1891 variable:
1893 >>> rv = wishart(df=1, scale=1)
1894 >>> # Frozen object with the same methods but holding the given
1895 >>> # degrees of freedom and scale fixed.
1897 """
1899 def __init__(self, seed=None):
1900 super().__init__(seed)
1901 self.__doc__ = doccer.docformat(self.__doc__, wishart_docdict_params)
1903 def __call__(self, df=None, scale=None, seed=None):
1904 """Create a frozen Wishart distribution.
1906 See `wishart_frozen` for more information.
1907 """
1908 return wishart_frozen(df, scale, seed)
1910 def _process_parameters(self, df, scale):
1911 if scale is None:
1912 scale = 1.0
1913 scale = np.asarray(scale, dtype=float)
1915 if scale.ndim == 0:
1916 scale = scale[np.newaxis, np.newaxis]
1917 elif scale.ndim == 1:
1918 scale = np.diag(scale)
1919 elif scale.ndim == 2 and not scale.shape[0] == scale.shape[1]:
1920 raise ValueError("Array 'scale' must be square if it is two"
1921 " dimensional, but scale.scale = %s."
1922 % str(scale.shape))
1923 elif scale.ndim > 2:
1924 raise ValueError("Array 'scale' must be at most two-dimensional,"
1925 " but scale.ndim = %d" % scale.ndim)
1927 dim = scale.shape[0]
1929 if df is None:
1930 df = dim
1931 elif not np.isscalar(df):
1932 raise ValueError("Degrees of freedom must be a scalar.")
1933 elif df <= dim - 1:
1934 raise ValueError("Degrees of freedom must be greater than the "
1935 "dimension of scale matrix minus 1.")
1937 return dim, df, scale
1939 def _process_quantiles(self, x, dim):
1940 """
1941 Adjust quantiles array so that last axis labels the components of
1942 each data point.
1943 """
1944 x = np.asarray(x, dtype=float)
1946 if x.ndim == 0:
1947 x = x * np.eye(dim)[:, :, np.newaxis]
1948 if x.ndim == 1:
1949 if dim == 1:
1950 x = x[np.newaxis, np.newaxis, :]
1951 else:
1952 x = np.diag(x)[:, :, np.newaxis]
1953 elif x.ndim == 2:
1954 if not x.shape[0] == x.shape[1]:
1955 raise ValueError("Quantiles must be square if they are two"
1956 " dimensional, but x.shape = %s."
1957 % str(x.shape))
1958 x = x[:, :, np.newaxis]
1959 elif x.ndim == 3:
1960 if not x.shape[0] == x.shape[1]:
1961 raise ValueError("Quantiles must be square in the first two"
1962 " dimensions if they are three dimensional"
1963 ", but x.shape = %s." % str(x.shape))
1964 elif x.ndim > 3:
1965 raise ValueError("Quantiles must be at most two-dimensional with"
1966 " an additional dimension for multiple"
1967 "components, but x.ndim = %d" % x.ndim)
1969 # Now we have 3-dim array; should have shape [dim, dim, *]
1970 if not x.shape[0:2] == (dim, dim):
1971 raise ValueError('Quantiles have incompatible dimensions: should'
1972 ' be {}, got {}.'.format((dim, dim), x.shape[0:2]))
1974 return x
1976 def _process_size(self, size):
1977 size = np.asarray(size)
1979 if size.ndim == 0:
1980 size = size[np.newaxis]
1981 elif size.ndim > 1:
1982 raise ValueError('Size must be an integer or tuple of integers;'
1983 ' thus must have dimension <= 1.'
1984 ' Got size.ndim = %s' % str(tuple(size)))
1985 n = size.prod()
1986 shape = tuple(size)
1988 return n, shape
1990 def _logpdf(self, x, dim, df, scale, log_det_scale, C):
1991 """Log of the Wishart probability density function.
1993 Parameters
1994 ----------
1995 x : ndarray
1996 Points at which to evaluate the log of the probability
1997 density function
1998 dim : int
1999 Dimension of the scale matrix
2000 df : int
2001 Degrees of freedom
2002 scale : ndarray
2003 Scale matrix
2004 log_det_scale : float
2005 Logarithm of the determinant of the scale matrix
2006 C : ndarray
2007 Cholesky factorization of the scale matrix, lower triagular.
2009 Notes
2010 -----
2011 As this function does no argument checking, it should not be
2012 called directly; use 'logpdf' instead.
2014 """
2015 # log determinant of x
2016 # Note: x has components along the last axis, so that x.T has
2017 # components alone the 0-th axis. Then since det(A) = det(A'), this
2018 # gives us a 1-dim vector of determinants
2020 # Retrieve tr(scale^{-1} x)
2021 log_det_x = np.empty(x.shape[-1])
2022 scale_inv_x = np.empty(x.shape)
2023 tr_scale_inv_x = np.empty(x.shape[-1])
2024 for i in range(x.shape[-1]):
2025 _, log_det_x[i] = self._cholesky_logdet(x[:, :, i])
2026 scale_inv_x[:, :, i] = scipy.linalg.cho_solve((C, True), x[:, :, i])
2027 tr_scale_inv_x[i] = scale_inv_x[:, :, i].trace()
2029 # Log PDF
2030 out = ((0.5 * (df - dim - 1) * log_det_x - 0.5 * tr_scale_inv_x) -
2031 (0.5 * df * dim * _LOG_2 + 0.5 * df * log_det_scale +
2032 multigammaln(0.5*df, dim)))
2034 return out
2036 def logpdf(self, x, df, scale):
2037 """Log of the Wishart probability density function.
2039 Parameters
2040 ----------
2041 x : array_like
2042 Quantiles, with the last axis of `x` denoting the components.
2043 Each quantile must be a symmetric positive definite matrix.
2044 %(_doc_default_callparams)s
2046 Returns
2047 -------
2048 pdf : ndarray
2049 Log of the probability density function evaluated at `x`
2051 Notes
2052 -----
2053 %(_doc_callparams_note)s
2055 """
2056 dim, df, scale = self._process_parameters(df, scale)
2057 x = self._process_quantiles(x, dim)
2059 # Cholesky decomposition of scale, get log(det(scale))
2060 C, log_det_scale = self._cholesky_logdet(scale)
2062 out = self._logpdf(x, dim, df, scale, log_det_scale, C)
2063 return _squeeze_output(out)
2065 def pdf(self, x, df, scale):
2066 """Wishart probability density function.
2068 Parameters
2069 ----------
2070 x : array_like
2071 Quantiles, with the last axis of `x` denoting the components.
2072 Each quantile must be a symmetric positive definite matrix.
2073 %(_doc_default_callparams)s
2075 Returns
2076 -------
2077 pdf : ndarray
2078 Probability density function evaluated at `x`
2080 Notes
2081 -----
2082 %(_doc_callparams_note)s
2084 """
2085 return np.exp(self.logpdf(x, df, scale))
2087 def _mean(self, dim, df, scale):
2088 """Mean of the Wishart distribution.
2090 Parameters
2091 ----------
2092 dim : int
2093 Dimension of the scale matrix
2094 %(_doc_default_callparams)s
2096 Notes
2097 -----
2098 As this function does no argument checking, it should not be
2099 called directly; use 'mean' instead.
2101 """
2102 return df * scale
2104 def mean(self, df, scale):
2105 """Mean of the Wishart distribution.
2107 Parameters
2108 ----------
2109 %(_doc_default_callparams)s
2111 Returns
2112 -------
2113 mean : float
2114 The mean of the distribution
2115 """
2116 dim, df, scale = self._process_parameters(df, scale)
2117 out = self._mean(dim, df, scale)
2118 return _squeeze_output(out)
2120 def _mode(self, dim, df, scale):
2121 """Mode of the Wishart distribution.
2123 Parameters
2124 ----------
2125 dim : int
2126 Dimension of the scale matrix
2127 %(_doc_default_callparams)s
2129 Notes
2130 -----
2131 As this function does no argument checking, it should not be
2132 called directly; use 'mode' instead.
2134 """
2135 if df >= dim + 1:
2136 out = (df-dim-1) * scale
2137 else:
2138 out = None
2139 return out
2141 def mode(self, df, scale):
2142 """Mode of the Wishart distribution
2144 Only valid if the degrees of freedom are greater than the dimension of
2145 the scale matrix.
2147 Parameters
2148 ----------
2149 %(_doc_default_callparams)s
2151 Returns
2152 -------
2153 mode : float or None
2154 The Mode of the distribution
2155 """
2156 dim, df, scale = self._process_parameters(df, scale)
2157 out = self._mode(dim, df, scale)
2158 return _squeeze_output(out) if out is not None else out
2160 def _var(self, dim, df, scale):
2161 """Variance of the Wishart distribution.
2163 Parameters
2164 ----------
2165 dim : int
2166 Dimension of the scale matrix
2167 %(_doc_default_callparams)s
2169 Notes
2170 -----
2171 As this function does no argument checking, it should not be
2172 called directly; use 'var' instead.
2174 """
2175 var = scale**2
2176 diag = scale.diagonal() # 1 x dim array
2177 var += np.outer(diag, diag)
2178 var *= df
2179 return var
2181 def var(self, df, scale):
2182 """Variance of the Wishart distribution.
2184 Parameters
2185 ----------
2186 %(_doc_default_callparams)s
2188 Returns
2189 -------
2190 var : float
2191 The variance of the distribution
2192 """
2193 dim, df, scale = self._process_parameters(df, scale)
2194 out = self._var(dim, df, scale)
2195 return _squeeze_output(out)
2197 def _standard_rvs(self, n, shape, dim, df, random_state):
2198 """
2199 Parameters
2200 ----------
2201 n : integer
2202 Number of variates to generate
2203 shape : iterable
2204 Shape of the variates to generate
2205 dim : int
2206 Dimension of the scale matrix
2207 df : int
2208 Degrees of freedom
2209 random_state : {None, int, `numpy.random.Generator`,
2210 `numpy.random.RandomState`}, optional
2212 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
2213 singleton is used.
2214 If `seed` is an int, a new ``RandomState`` instance is used,
2215 seeded with `seed`.
2216 If `seed` is already a ``Generator`` or ``RandomState`` instance
2217 then that instance is used.
2219 Notes
2220 -----
2221 As this function does no argument checking, it should not be
2222 called directly; use 'rvs' instead.
2224 """
2225 # Random normal variates for off-diagonal elements
2226 n_tril = dim * (dim-1) // 2
2227 covariances = random_state.normal(
2228 size=n*n_tril).reshape(shape+(n_tril,))
2230 # Random chi-square variates for diagonal elements
2231 variances = (np.r_[[random_state.chisquare(df-(i+1)+1, size=n)**0.5
2232 for i in range(dim)]].reshape((dim,) +
2233 shape[::-1]).T)
2235 # Create the A matri(ces) - lower triangular
2236 A = np.zeros(shape + (dim, dim))
2238 # Input the covariances
2239 size_idx = tuple([slice(None, None, None)]*len(shape))
2240 tril_idx = np.tril_indices(dim, k=-1)
2241 A[size_idx + tril_idx] = covariances
2243 # Input the variances
2244 diag_idx = np.diag_indices(dim)
2245 A[size_idx + diag_idx] = variances
2247 return A
2249 def _rvs(self, n, shape, dim, df, C, random_state):
2250 """Draw random samples from a Wishart distribution.
2252 Parameters
2253 ----------
2254 n : integer
2255 Number of variates to generate
2256 shape : iterable
2257 Shape of the variates to generate
2258 dim : int
2259 Dimension of the scale matrix
2260 df : int
2261 Degrees of freedom
2262 C : ndarray
2263 Cholesky factorization of the scale matrix, lower triangular.
2264 %(_doc_random_state)s
2266 Notes
2267 -----
2268 As this function does no argument checking, it should not be
2269 called directly; use 'rvs' instead.
2271 """
2272 random_state = self._get_random_state(random_state)
2273 # Calculate the matrices A, which are actually lower triangular
2274 # Cholesky factorizations of a matrix B such that B ~ W(df, I)
2275 A = self._standard_rvs(n, shape, dim, df, random_state)
2277 # Calculate SA = C A A' C', where SA ~ W(df, scale)
2278 # Note: this is the product of a (lower) (lower) (lower)' (lower)'
2279 # or, denoting B = AA', it is C B C' where C is the lower
2280 # triangular Cholesky factorization of the scale matrix.
2281 # this appears to conflict with the instructions in [1]_, which
2282 # suggest that it should be D' B D where D is the lower
2283 # triangular factorization of the scale matrix. However, it is
2284 # meant to refer to the Bartlett (1933) representation of a
2285 # Wishart random variate as L A A' L' where L is lower triangular
2286 # so it appears that understanding D' to be upper triangular
2287 # is either a typo in or misreading of [1]_.
2288 for index in np.ndindex(shape):
2289 CA = np.dot(C, A[index])
2290 A[index] = np.dot(CA, CA.T)
2292 return A
2294 def rvs(self, df, scale, size=1, random_state=None):
2295 """Draw random samples from a Wishart distribution.
2297 Parameters
2298 ----------
2299 %(_doc_default_callparams)s
2300 size : integer or iterable of integers, optional
2301 Number of samples to draw (default 1).
2302 %(_doc_random_state)s
2304 Returns
2305 -------
2306 rvs : ndarray
2307 Random variates of shape (`size`) + (`dim`, `dim), where `dim` is
2308 the dimension of the scale matrix.
2310 Notes
2311 -----
2312 %(_doc_callparams_note)s
2314 """
2315 n, shape = self._process_size(size)
2316 dim, df, scale = self._process_parameters(df, scale)
2318 # Cholesky decomposition of scale
2319 C = scipy.linalg.cholesky(scale, lower=True)
2321 out = self._rvs(n, shape, dim, df, C, random_state)
2323 return _squeeze_output(out)
2325 def _entropy(self, dim, df, log_det_scale):
2326 """Compute the differential entropy of the Wishart.
2328 Parameters
2329 ----------
2330 dim : int
2331 Dimension of the scale matrix
2332 df : int
2333 Degrees of freedom
2334 log_det_scale : float
2335 Logarithm of the determinant of the scale matrix
2337 Notes
2338 -----
2339 As this function does no argument checking, it should not be
2340 called directly; use 'entropy' instead.
2342 """
2343 return (
2344 0.5 * (dim+1) * log_det_scale +
2345 0.5 * dim * (dim+1) * _LOG_2 +
2346 multigammaln(0.5*df, dim) -
2347 0.5 * (df - dim - 1) * np.sum(
2348 [psi(0.5*(df + 1 - (i+1))) for i in range(dim)]
2349 ) +
2350 0.5 * df * dim
2351 )
2353 def entropy(self, df, scale):
2354 """Compute the differential entropy of the Wishart.
2356 Parameters
2357 ----------
2358 %(_doc_default_callparams)s
2360 Returns
2361 -------
2362 h : scalar
2363 Entropy of the Wishart distribution
2365 Notes
2366 -----
2367 %(_doc_callparams_note)s
2369 """
2370 dim, df, scale = self._process_parameters(df, scale)
2371 _, log_det_scale = self._cholesky_logdet(scale)
2372 return self._entropy(dim, df, log_det_scale)
2374 def _cholesky_logdet(self, scale):
2375 """Compute Cholesky decomposition and determine (log(det(scale)).
2377 Parameters
2378 ----------
2379 scale : ndarray
2380 Scale matrix.
2382 Returns
2383 -------
2384 c_decomp : ndarray
2385 The Cholesky decomposition of `scale`.
2386 logdet : scalar
2387 The log of the determinant of `scale`.
2389 Notes
2390 -----
2391 This computation of ``logdet`` is equivalent to
2392 ``np.linalg.slogdet(scale)``. It is ~2x faster though.
2394 """
2395 c_decomp = scipy.linalg.cholesky(scale, lower=True)
2396 logdet = 2 * np.sum(np.log(c_decomp.diagonal()))
2397 return c_decomp, logdet
2400wishart = wishart_gen()
2403class wishart_frozen(multi_rv_frozen):
2404 """Create a frozen Wishart distribution.
2406 Parameters
2407 ----------
2408 df : array_like
2409 Degrees of freedom of the distribution
2410 scale : array_like
2411 Scale matrix of the distribution
2412 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
2413 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
2414 singleton is used.
2415 If `seed` is an int, a new ``RandomState`` instance is used,
2416 seeded with `seed`.
2417 If `seed` is already a ``Generator`` or ``RandomState`` instance then
2418 that instance is used.
2420 """
2421 def __init__(self, df, scale, seed=None):
2422 self._dist = wishart_gen(seed)
2423 self.dim, self.df, self.scale = self._dist._process_parameters(
2424 df, scale)
2425 self.C, self.log_det_scale = self._dist._cholesky_logdet(self.scale)
2427 def logpdf(self, x):
2428 x = self._dist._process_quantiles(x, self.dim)
2430 out = self._dist._logpdf(x, self.dim, self.df, self.scale,
2431 self.log_det_scale, self.C)
2432 return _squeeze_output(out)
2434 def pdf(self, x):
2435 return np.exp(self.logpdf(x))
2437 def mean(self):
2438 out = self._dist._mean(self.dim, self.df, self.scale)
2439 return _squeeze_output(out)
2441 def mode(self):
2442 out = self._dist._mode(self.dim, self.df, self.scale)
2443 return _squeeze_output(out) if out is not None else out
2445 def var(self):
2446 out = self._dist._var(self.dim, self.df, self.scale)
2447 return _squeeze_output(out)
2449 def rvs(self, size=1, random_state=None):
2450 n, shape = self._dist._process_size(size)
2451 out = self._dist._rvs(n, shape, self.dim, self.df,
2452 self.C, random_state)
2453 return _squeeze_output(out)
2455 def entropy(self):
2456 return self._dist._entropy(self.dim, self.df, self.log_det_scale)
2459# Set frozen generator docstrings from corresponding docstrings in
2460# Wishart and fill in default strings in class docstrings
2461for name in ['logpdf', 'pdf', 'mean', 'mode', 'var', 'rvs', 'entropy']:
2462 method = wishart_gen.__dict__[name]
2463 method_frozen = wishart_frozen.__dict__[name]
2464 method_frozen.__doc__ = doccer.docformat(
2465 method.__doc__, wishart_docdict_noparams)
2466 method.__doc__ = doccer.docformat(method.__doc__, wishart_docdict_params)
2469def _cho_inv_batch(a, check_finite=True):
2470 """
2471 Invert the matrices a_i, using a Cholesky factorization of A, where
2472 a_i resides in the last two dimensions of a and the other indices describe
2473 the index i.
2475 Overwrites the data in a.
2477 Parameters
2478 ----------
2479 a : array
2480 Array of matrices to invert, where the matrices themselves are stored
2481 in the last two dimensions.
2482 check_finite : bool, optional
2483 Whether to check that the input matrices contain only finite numbers.
2484 Disabling may give a performance gain, but may result in problems
2485 (crashes, non-termination) if the inputs do contain infinities or NaNs.
2487 Returns
2488 -------
2489 x : array
2490 Array of inverses of the matrices ``a_i``.
2492 See Also
2493 --------
2494 scipy.linalg.cholesky : Cholesky factorization of a matrix
2496 """
2497 if check_finite:
2498 a1 = asarray_chkfinite(a)
2499 else:
2500 a1 = asarray(a)
2501 if len(a1.shape) < 2 or a1.shape[-2] != a1.shape[-1]:
2502 raise ValueError('expected square matrix in last two dimensions')
2504 potrf, potri = get_lapack_funcs(('potrf', 'potri'), (a1,))
2506 triu_rows, triu_cols = np.triu_indices(a.shape[-2], k=1)
2507 for index in np.ndindex(a1.shape[:-2]):
2509 # Cholesky decomposition
2510 a1[index], info = potrf(a1[index], lower=True, overwrite_a=False,
2511 clean=False)
2512 if info > 0:
2513 raise LinAlgError("%d-th leading minor not positive definite"
2514 % info)
2515 if info < 0:
2516 raise ValueError('illegal value in %d-th argument of internal'
2517 ' potrf' % -info)
2518 # Inversion
2519 a1[index], info = potri(a1[index], lower=True, overwrite_c=False)
2520 if info > 0:
2521 raise LinAlgError("the inverse could not be computed")
2522 if info < 0:
2523 raise ValueError('illegal value in %d-th argument of internal'
2524 ' potrf' % -info)
2526 # Make symmetric (dpotri only fills in the lower triangle)
2527 a1[index][triu_rows, triu_cols] = a1[index][triu_cols, triu_rows]
2529 return a1
2532class invwishart_gen(wishart_gen):
2533 r"""An inverse Wishart random variable.
2535 The `df` keyword specifies the degrees of freedom. The `scale` keyword
2536 specifies the scale matrix, which must be symmetric and positive definite.
2537 In this context, the scale matrix is often interpreted in terms of a
2538 multivariate normal covariance matrix.
2540 Methods
2541 -------
2542 pdf(x, df, scale)
2543 Probability density function.
2544 logpdf(x, df, scale)
2545 Log of the probability density function.
2546 rvs(df, scale, size=1, random_state=None)
2547 Draw random samples from an inverse Wishart distribution.
2548 entropy(df, scale)
2549 Differential entropy of the distribution.
2551 Parameters
2552 ----------
2553 %(_doc_default_callparams)s
2554 %(_doc_random_state)s
2556 Raises
2557 ------
2558 scipy.linalg.LinAlgError
2559 If the scale matrix `scale` is not positive definite.
2561 See Also
2562 --------
2563 wishart
2565 Notes
2566 -----
2567 %(_doc_callparams_note)s
2569 The scale matrix `scale` must be a symmetric positive definite
2570 matrix. Singular matrices, including the symmetric positive semi-definite
2571 case, are not supported. Symmetry is not checked; only the lower triangular
2572 portion is used.
2574 The inverse Wishart distribution is often denoted
2576 .. math::
2578 W_p^{-1}(\nu, \Psi)
2580 where :math:`\nu` is the degrees of freedom and :math:`\Psi` is the
2581 :math:`p \times p` scale matrix.
2583 The probability density function for `invwishart` has support over positive
2584 definite matrices :math:`S`; if :math:`S \sim W^{-1}_p(\nu, \Sigma)`,
2585 then its PDF is given by:
2587 .. math::
2589 f(S) = \frac{|\Sigma|^\frac{\nu}{2}}{2^{ \frac{\nu p}{2} }
2590 |S|^{\frac{\nu + p + 1}{2}} \Gamma_p \left(\frac{\nu}{2} \right)}
2591 \exp\left( -tr(\Sigma S^{-1}) / 2 \right)
2593 If :math:`S \sim W_p^{-1}(\nu, \Psi)` (inverse Wishart) then
2594 :math:`S^{-1} \sim W_p(\nu, \Psi^{-1})` (Wishart).
2596 If the scale matrix is 1-dimensional and equal to one, then the inverse
2597 Wishart distribution :math:`W_1(\nu, 1)` collapses to the
2598 inverse Gamma distribution with parameters shape = :math:`\frac{\nu}{2}`
2599 and scale = :math:`\frac{1}{2}`.
2601 .. versionadded:: 0.16.0
2603 References
2604 ----------
2605 .. [1] M.L. Eaton, "Multivariate Statistics: A Vector Space Approach",
2606 Wiley, 1983.
2607 .. [2] M.C. Jones, "Generating Inverse Wishart Matrices", Communications
2608 in Statistics - Simulation and Computation, vol. 14.2, pp.511-514,
2609 1985.
2610 .. [3] Gupta, M. and Srivastava, S. "Parametric Bayesian Estimation of
2611 Differential Entropy and Relative Entropy". Entropy 12, 818 - 843.
2612 2010.
2614 Examples
2615 --------
2616 >>> import numpy as np
2617 >>> import matplotlib.pyplot as plt
2618 >>> from scipy.stats import invwishart, invgamma
2619 >>> x = np.linspace(0.01, 1, 100)
2620 >>> iw = invwishart.pdf(x, df=6, scale=1)
2621 >>> iw[:3]
2622 array([ 1.20546865e-15, 5.42497807e-06, 4.45813929e-03])
2623 >>> ig = invgamma.pdf(x, 6/2., scale=1./2)
2624 >>> ig[:3]
2625 array([ 1.20546865e-15, 5.42497807e-06, 4.45813929e-03])
2626 >>> plt.plot(x, iw)
2627 >>> plt.show()
2629 The input quantiles can be any shape of array, as long as the last
2630 axis labels the components.
2632 Alternatively, the object may be called (as a function) to fix the degrees
2633 of freedom and scale parameters, returning a "frozen" inverse Wishart
2634 random variable:
2636 >>> rv = invwishart(df=1, scale=1)
2637 >>> # Frozen object with the same methods but holding the given
2638 >>> # degrees of freedom and scale fixed.
2640 """
2642 def __init__(self, seed=None):
2643 super().__init__(seed)
2644 self.__doc__ = doccer.docformat(self.__doc__, wishart_docdict_params)
2646 def __call__(self, df=None, scale=None, seed=None):
2647 """Create a frozen inverse Wishart distribution.
2649 See `invwishart_frozen` for more information.
2651 """
2652 return invwishart_frozen(df, scale, seed)
2654 def _logpdf(self, x, dim, df, scale, log_det_scale):
2655 """Log of the inverse Wishart probability density function.
2657 Parameters
2658 ----------
2659 x : ndarray
2660 Points at which to evaluate the log of the probability
2661 density function.
2662 dim : int
2663 Dimension of the scale matrix
2664 df : int
2665 Degrees of freedom
2666 scale : ndarray
2667 Scale matrix
2668 log_det_scale : float
2669 Logarithm of the determinant of the scale matrix
2671 Notes
2672 -----
2673 As this function does no argument checking, it should not be
2674 called directly; use 'logpdf' instead.
2676 """
2677 log_det_x = np.empty(x.shape[-1])
2678 x_inv = np.copy(x).T
2679 if dim > 1:
2680 _cho_inv_batch(x_inv) # works in-place
2681 else:
2682 x_inv = 1./x_inv
2683 tr_scale_x_inv = np.empty(x.shape[-1])
2685 for i in range(x.shape[-1]):
2686 C, lower = scipy.linalg.cho_factor(x[:, :, i], lower=True)
2688 log_det_x[i] = 2 * np.sum(np.log(C.diagonal()))
2690 tr_scale_x_inv[i] = np.dot(scale, x_inv[i]).trace()
2692 # Log PDF
2693 out = ((0.5 * df * log_det_scale - 0.5 * tr_scale_x_inv) -
2694 (0.5 * df * dim * _LOG_2 + 0.5 * (df + dim + 1) * log_det_x) -
2695 multigammaln(0.5*df, dim))
2697 return out
2699 def logpdf(self, x, df, scale):
2700 """Log of the inverse Wishart probability density function.
2702 Parameters
2703 ----------
2704 x : array_like
2705 Quantiles, with the last axis of `x` denoting the components.
2706 Each quantile must be a symmetric positive definite matrix.
2707 %(_doc_default_callparams)s
2709 Returns
2710 -------
2711 pdf : ndarray
2712 Log of the probability density function evaluated at `x`
2714 Notes
2715 -----
2716 %(_doc_callparams_note)s
2718 """
2719 dim, df, scale = self._process_parameters(df, scale)
2720 x = self._process_quantiles(x, dim)
2721 _, log_det_scale = self._cholesky_logdet(scale)
2722 out = self._logpdf(x, dim, df, scale, log_det_scale)
2723 return _squeeze_output(out)
2725 def pdf(self, x, df, scale):
2726 """Inverse Wishart probability density function.
2728 Parameters
2729 ----------
2730 x : array_like
2731 Quantiles, with the last axis of `x` denoting the components.
2732 Each quantile must be a symmetric positive definite matrix.
2733 %(_doc_default_callparams)s
2735 Returns
2736 -------
2737 pdf : ndarray
2738 Probability density function evaluated at `x`
2740 Notes
2741 -----
2742 %(_doc_callparams_note)s
2744 """
2745 return np.exp(self.logpdf(x, df, scale))
2747 def _mean(self, dim, df, scale):
2748 """Mean of the inverse Wishart distribution.
2750 Parameters
2751 ----------
2752 dim : int
2753 Dimension of the scale matrix
2754 %(_doc_default_callparams)s
2756 Notes
2757 -----
2758 As this function does no argument checking, it should not be
2759 called directly; use 'mean' instead.
2761 """
2762 if df > dim + 1:
2763 out = scale / (df - dim - 1)
2764 else:
2765 out = None
2766 return out
2768 def mean(self, df, scale):
2769 """Mean of the inverse Wishart distribution.
2771 Only valid if the degrees of freedom are greater than the dimension of
2772 the scale matrix plus one.
2774 Parameters
2775 ----------
2776 %(_doc_default_callparams)s
2778 Returns
2779 -------
2780 mean : float or None
2781 The mean of the distribution
2783 """
2784 dim, df, scale = self._process_parameters(df, scale)
2785 out = self._mean(dim, df, scale)
2786 return _squeeze_output(out) if out is not None else out
2788 def _mode(self, dim, df, scale):
2789 """Mode of the inverse Wishart distribution.
2791 Parameters
2792 ----------
2793 dim : int
2794 Dimension of the scale matrix
2795 %(_doc_default_callparams)s
2797 Notes
2798 -----
2799 As this function does no argument checking, it should not be
2800 called directly; use 'mode' instead.
2802 """
2803 return scale / (df + dim + 1)
2805 def mode(self, df, scale):
2806 """Mode of the inverse Wishart distribution.
2808 Parameters
2809 ----------
2810 %(_doc_default_callparams)s
2812 Returns
2813 -------
2814 mode : float
2815 The Mode of the distribution
2817 """
2818 dim, df, scale = self._process_parameters(df, scale)
2819 out = self._mode(dim, df, scale)
2820 return _squeeze_output(out)
2822 def _var(self, dim, df, scale):
2823 """Variance of the inverse Wishart distribution.
2825 Parameters
2826 ----------
2827 dim : int
2828 Dimension of the scale matrix
2829 %(_doc_default_callparams)s
2831 Notes
2832 -----
2833 As this function does no argument checking, it should not be
2834 called directly; use 'var' instead.
2836 """
2837 if df > dim + 3:
2838 var = (df - dim + 1) * scale**2
2839 diag = scale.diagonal() # 1 x dim array
2840 var += (df - dim - 1) * np.outer(diag, diag)
2841 var /= (df - dim) * (df - dim - 1)**2 * (df - dim - 3)
2842 else:
2843 var = None
2844 return var
2846 def var(self, df, scale):
2847 """Variance of the inverse Wishart distribution.
2849 Only valid if the degrees of freedom are greater than the dimension of
2850 the scale matrix plus three.
2852 Parameters
2853 ----------
2854 %(_doc_default_callparams)s
2856 Returns
2857 -------
2858 var : float
2859 The variance of the distribution
2860 """
2861 dim, df, scale = self._process_parameters(df, scale)
2862 out = self._var(dim, df, scale)
2863 return _squeeze_output(out) if out is not None else out
2865 def _rvs(self, n, shape, dim, df, C, random_state):
2866 """Draw random samples from an inverse Wishart distribution.
2868 Parameters
2869 ----------
2870 n : integer
2871 Number of variates to generate
2872 shape : iterable
2873 Shape of the variates to generate
2874 dim : int
2875 Dimension of the scale matrix
2876 df : int
2877 Degrees of freedom
2878 C : ndarray
2879 Cholesky factorization of the scale matrix, lower triagular.
2880 %(_doc_random_state)s
2882 Notes
2883 -----
2884 As this function does no argument checking, it should not be
2885 called directly; use 'rvs' instead.
2887 """
2888 random_state = self._get_random_state(random_state)
2889 # Get random draws A such that A ~ W(df, I)
2890 A = super()._standard_rvs(n, shape, dim, df, random_state)
2892 # Calculate SA = (CA)'^{-1} (CA)^{-1} ~ iW(df, scale)
2893 eye = np.eye(dim)
2894 trtrs = get_lapack_funcs(('trtrs'), (A,))
2896 for index in np.ndindex(A.shape[:-2]):
2897 # Calculate CA
2898 CA = np.dot(C, A[index])
2899 # Get (C A)^{-1} via triangular solver
2900 if dim > 1:
2901 CA, info = trtrs(CA, eye, lower=True)
2902 if info > 0:
2903 raise LinAlgError("Singular matrix.")
2904 if info < 0:
2905 raise ValueError('Illegal value in %d-th argument of'
2906 ' internal trtrs' % -info)
2907 else:
2908 CA = 1. / CA
2909 # Get SA
2910 A[index] = np.dot(CA.T, CA)
2912 return A
2914 def rvs(self, df, scale, size=1, random_state=None):
2915 """Draw random samples from an inverse Wishart distribution.
2917 Parameters
2918 ----------
2919 %(_doc_default_callparams)s
2920 size : integer or iterable of integers, optional
2921 Number of samples to draw (default 1).
2922 %(_doc_random_state)s
2924 Returns
2925 -------
2926 rvs : ndarray
2927 Random variates of shape (`size`) + (`dim`, `dim), where `dim` is
2928 the dimension of the scale matrix.
2930 Notes
2931 -----
2932 %(_doc_callparams_note)s
2934 """
2935 n, shape = self._process_size(size)
2936 dim, df, scale = self._process_parameters(df, scale)
2938 # Invert the scale
2939 eye = np.eye(dim)
2940 L, lower = scipy.linalg.cho_factor(scale, lower=True)
2941 inv_scale = scipy.linalg.cho_solve((L, lower), eye)
2942 # Cholesky decomposition of inverted scale
2943 C = scipy.linalg.cholesky(inv_scale, lower=True)
2945 out = self._rvs(n, shape, dim, df, C, random_state)
2947 return _squeeze_output(out)
2949 def _entropy(self, dim, df, log_det_scale):
2950 # reference: eq. (17) from ref. 3
2951 psi_eval_points = [0.5 * (df - dim + i) for i in range(1, dim + 1)]
2952 psi_eval_points = np.asarray(psi_eval_points)
2953 return multigammaln(0.5 * df, dim) + 0.5 * dim * df + \
2954 0.5 * (dim + 1) * (log_det_scale - _LOG_2) - \
2955 0.5 * (df + dim + 1) * \
2956 psi(psi_eval_points, out=psi_eval_points).sum()
2958 def entropy(self, df, scale):
2959 dim, df, scale = self._process_parameters(df, scale)
2960 _, log_det_scale = self._cholesky_logdet(scale)
2961 return self._entropy(dim, df, log_det_scale)
2964invwishart = invwishart_gen()
2967class invwishart_frozen(multi_rv_frozen):
2968 def __init__(self, df, scale, seed=None):
2969 """Create a frozen inverse Wishart distribution.
2971 Parameters
2972 ----------
2973 df : array_like
2974 Degrees of freedom of the distribution
2975 scale : array_like
2976 Scale matrix of the distribution
2977 seed : {None, int, `numpy.random.Generator`}, optional
2978 If `seed` is None the `numpy.random.Generator` singleton is used.
2979 If `seed` is an int, a new ``Generator`` instance is used,
2980 seeded with `seed`.
2981 If `seed` is already a ``Generator`` instance then that instance is
2982 used.
2984 """
2985 self._dist = invwishart_gen(seed)
2986 self.dim, self.df, self.scale = self._dist._process_parameters(
2987 df, scale
2988 )
2990 # Get the determinant via Cholesky factorization
2991 C, lower = scipy.linalg.cho_factor(self.scale, lower=True)
2992 self.log_det_scale = 2 * np.sum(np.log(C.diagonal()))
2994 # Get the inverse using the Cholesky factorization
2995 eye = np.eye(self.dim)
2996 self.inv_scale = scipy.linalg.cho_solve((C, lower), eye)
2998 # Get the Cholesky factorization of the inverse scale
2999 self.C = scipy.linalg.cholesky(self.inv_scale, lower=True)
3001 def logpdf(self, x):
3002 x = self._dist._process_quantiles(x, self.dim)
3003 out = self._dist._logpdf(x, self.dim, self.df, self.scale,
3004 self.log_det_scale)
3005 return _squeeze_output(out)
3007 def pdf(self, x):
3008 return np.exp(self.logpdf(x))
3010 def mean(self):
3011 out = self._dist._mean(self.dim, self.df, self.scale)
3012 return _squeeze_output(out) if out is not None else out
3014 def mode(self):
3015 out = self._dist._mode(self.dim, self.df, self.scale)
3016 return _squeeze_output(out)
3018 def var(self):
3019 out = self._dist._var(self.dim, self.df, self.scale)
3020 return _squeeze_output(out) if out is not None else out
3022 def rvs(self, size=1, random_state=None):
3023 n, shape = self._dist._process_size(size)
3025 out = self._dist._rvs(n, shape, self.dim, self.df,
3026 self.C, random_state)
3028 return _squeeze_output(out)
3030 def entropy(self):
3031 return self._dist._entropy(self.dim, self.df, self.log_det_scale)
3034# Set frozen generator docstrings from corresponding docstrings in
3035# inverse Wishart and fill in default strings in class docstrings
3036for name in ['logpdf', 'pdf', 'mean', 'mode', 'var', 'rvs']:
3037 method = invwishart_gen.__dict__[name]
3038 method_frozen = wishart_frozen.__dict__[name]
3039 method_frozen.__doc__ = doccer.docformat(
3040 method.__doc__, wishart_docdict_noparams)
3041 method.__doc__ = doccer.docformat(method.__doc__, wishart_docdict_params)
3043_multinomial_doc_default_callparams = """\
3044n : int
3045 Number of trials
3046p : array_like
3047 Probability of a trial falling into each category; should sum to 1
3048"""
3050_multinomial_doc_callparams_note = """\
3051`n` should be a nonnegative integer. Each element of `p` should be in the
3052interval :math:`[0,1]` and the elements should sum to 1. If they do not sum to
30531, the last element of the `p` array is not used and is replaced with the
3054remaining probability left over from the earlier elements.
3055"""
3057_multinomial_doc_frozen_callparams = ""
3059_multinomial_doc_frozen_callparams_note = """\
3060See class definition for a detailed description of parameters."""
3062multinomial_docdict_params = {
3063 '_doc_default_callparams': _multinomial_doc_default_callparams,
3064 '_doc_callparams_note': _multinomial_doc_callparams_note,
3065 '_doc_random_state': _doc_random_state
3066}
3068multinomial_docdict_noparams = {
3069 '_doc_default_callparams': _multinomial_doc_frozen_callparams,
3070 '_doc_callparams_note': _multinomial_doc_frozen_callparams_note,
3071 '_doc_random_state': _doc_random_state
3072}
3075class multinomial_gen(multi_rv_generic):
3076 r"""A multinomial random variable.
3078 Methods
3079 -------
3080 pmf(x, n, p)
3081 Probability mass function.
3082 logpmf(x, n, p)
3083 Log of the probability mass function.
3084 rvs(n, p, size=1, random_state=None)
3085 Draw random samples from a multinomial distribution.
3086 entropy(n, p)
3087 Compute the entropy of the multinomial distribution.
3088 cov(n, p)
3089 Compute the covariance matrix of the multinomial distribution.
3091 Parameters
3092 ----------
3093 %(_doc_default_callparams)s
3094 %(_doc_random_state)s
3096 Notes
3097 -----
3098 %(_doc_callparams_note)s
3100 The probability mass function for `multinomial` is
3102 .. math::
3104 f(x) = \frac{n!}{x_1! \cdots x_k!} p_1^{x_1} \cdots p_k^{x_k},
3106 supported on :math:`x=(x_1, \ldots, x_k)` where each :math:`x_i` is a
3107 nonnegative integer and their sum is :math:`n`.
3109 .. versionadded:: 0.19.0
3111 Examples
3112 --------
3114 >>> from scipy.stats import multinomial
3115 >>> rv = multinomial(8, [0.3, 0.2, 0.5])
3116 >>> rv.pmf([1, 3, 4])
3117 0.042000000000000072
3119 The multinomial distribution for :math:`k=2` is identical to the
3120 corresponding binomial distribution (tiny numerical differences
3121 notwithstanding):
3123 >>> from scipy.stats import binom
3124 >>> multinomial.pmf([3, 4], n=7, p=[0.4, 0.6])
3125 0.29030399999999973
3126 >>> binom.pmf(3, 7, 0.4)
3127 0.29030400000000012
3129 The functions ``pmf``, ``logpmf``, ``entropy``, and ``cov`` support
3130 broadcasting, under the convention that the vector parameters (``x`` and
3131 ``p``) are interpreted as if each row along the last axis is a single
3132 object. For instance:
3134 >>> multinomial.pmf([[3, 4], [3, 5]], n=[7, 8], p=[.3, .7])
3135 array([0.2268945, 0.25412184])
3137 Here, ``x.shape == (2, 2)``, ``n.shape == (2,)``, and ``p.shape == (2,)``,
3138 but following the rules mentioned above they behave as if the rows
3139 ``[3, 4]`` and ``[3, 5]`` in ``x`` and ``[.3, .7]`` in ``p`` were a single
3140 object, and as if we had ``x.shape = (2,)``, ``n.shape = (2,)``, and
3141 ``p.shape = ()``. To obtain the individual elements without broadcasting,
3142 we would do this:
3144 >>> multinomial.pmf([3, 4], n=7, p=[.3, .7])
3145 0.2268945
3146 >>> multinomial.pmf([3, 5], 8, p=[.3, .7])
3147 0.25412184
3149 This broadcasting also works for ``cov``, where the output objects are
3150 square matrices of size ``p.shape[-1]``. For example:
3152 >>> multinomial.cov([4, 5], [[.3, .7], [.4, .6]])
3153 array([[[ 0.84, -0.84],
3154 [-0.84, 0.84]],
3155 [[ 1.2 , -1.2 ],
3156 [-1.2 , 1.2 ]]])
3158 In this example, ``n.shape == (2,)`` and ``p.shape == (2, 2)``, and
3159 following the rules above, these broadcast as if ``p.shape == (2,)``.
3160 Thus the result should also be of shape ``(2,)``, but since each output is
3161 a :math:`2 \times 2` matrix, the result in fact has shape ``(2, 2, 2)``,
3162 where ``result[0]`` is equal to ``multinomial.cov(n=4, p=[.3, .7])`` and
3163 ``result[1]`` is equal to ``multinomial.cov(n=5, p=[.4, .6])``.
3165 Alternatively, the object may be called (as a function) to fix the `n` and
3166 `p` parameters, returning a "frozen" multinomial random variable:
3168 >>> rv = multinomial(n=7, p=[.3, .7])
3169 >>> # Frozen object with the same methods but holding the given
3170 >>> # degrees of freedom and scale fixed.
3172 See also
3173 --------
3174 scipy.stats.binom : The binomial distribution.
3175 numpy.random.Generator.multinomial : Sampling from the multinomial distribution.
3176 scipy.stats.multivariate_hypergeom :
3177 The multivariate hypergeometric distribution.
3178 """ # noqa: E501
3180 def __init__(self, seed=None):
3181 super().__init__(seed)
3182 self.__doc__ = \
3183 doccer.docformat(self.__doc__, multinomial_docdict_params)
3185 def __call__(self, n, p, seed=None):
3186 """Create a frozen multinomial distribution.
3188 See `multinomial_frozen` for more information.
3189 """
3190 return multinomial_frozen(n, p, seed)
3192 def _process_parameters(self, n, p, eps=1e-15):
3193 """Returns: n_, p_, npcond.
3195 n_ and p_ are arrays of the correct shape; npcond is a boolean array
3196 flagging values out of the domain.
3197 """
3198 p = np.array(p, dtype=np.float64, copy=True)
3199 p_adjusted = 1. - p[..., :-1].sum(axis=-1)
3200 i_adjusted = np.abs(p_adjusted) > eps
3201 p[i_adjusted, -1] = p_adjusted[i_adjusted]
3203 # true for bad p
3204 pcond = np.any(p < 0, axis=-1)
3205 pcond |= np.any(p > 1, axis=-1)
3207 n = np.array(n, dtype=np.int_, copy=True)
3209 # true for bad n
3210 ncond = n < 0
3212 return n, p, ncond | pcond
3214 def _process_quantiles(self, x, n, p):
3215 """Returns: x_, xcond.
3217 x_ is an int array; xcond is a boolean array flagging values out of the
3218 domain.
3219 """
3220 xx = np.asarray(x, dtype=np.int_)
3222 if xx.ndim == 0:
3223 raise ValueError("x must be an array.")
3225 if xx.size != 0 and not xx.shape[-1] == p.shape[-1]:
3226 raise ValueError("Size of each quantile should be size of p: "
3227 "received %d, but expected %d." %
3228 (xx.shape[-1], p.shape[-1]))
3230 # true for x out of the domain
3231 cond = np.any(xx != x, axis=-1)
3232 cond |= np.any(xx < 0, axis=-1)
3233 cond = cond | (np.sum(xx, axis=-1) != n)
3235 return xx, cond
3237 def _checkresult(self, result, cond, bad_value):
3238 result = np.asarray(result)
3240 if cond.ndim != 0:
3241 result[cond] = bad_value
3242 elif cond:
3243 if result.ndim == 0:
3244 return bad_value
3245 result[...] = bad_value
3246 return result
3248 def _logpmf(self, x, n, p):
3249 return gammaln(n+1) + np.sum(xlogy(x, p) - gammaln(x+1), axis=-1)
3251 def logpmf(self, x, n, p):
3252 """Log of the Multinomial probability mass function.
3254 Parameters
3255 ----------
3256 x : array_like
3257 Quantiles, with the last axis of `x` denoting the components.
3258 %(_doc_default_callparams)s
3260 Returns
3261 -------
3262 logpmf : ndarray or scalar
3263 Log of the probability mass function evaluated at `x`
3265 Notes
3266 -----
3267 %(_doc_callparams_note)s
3268 """
3269 n, p, npcond = self._process_parameters(n, p)
3270 x, xcond = self._process_quantiles(x, n, p)
3272 result = self._logpmf(x, n, p)
3274 # replace values for which x was out of the domain; broadcast
3275 # xcond to the right shape
3276 xcond_ = xcond | np.zeros(npcond.shape, dtype=np.bool_)
3277 result = self._checkresult(result, xcond_, -np.inf)
3279 # replace values bad for n or p; broadcast npcond to the right shape
3280 npcond_ = npcond | np.zeros(xcond.shape, dtype=np.bool_)
3281 return self._checkresult(result, npcond_, np.nan)
3283 def pmf(self, x, n, p):
3284 """Multinomial probability mass function.
3286 Parameters
3287 ----------
3288 x : array_like
3289 Quantiles, with the last axis of `x` denoting the components.
3290 %(_doc_default_callparams)s
3292 Returns
3293 -------
3294 pmf : ndarray or scalar
3295 Probability density function evaluated at `x`
3297 Notes
3298 -----
3299 %(_doc_callparams_note)s
3300 """
3301 return np.exp(self.logpmf(x, n, p))
3303 def mean(self, n, p):
3304 """Mean of the Multinomial distribution.
3306 Parameters
3307 ----------
3308 %(_doc_default_callparams)s
3310 Returns
3311 -------
3312 mean : float
3313 The mean of the distribution
3314 """
3315 n, p, npcond = self._process_parameters(n, p)
3316 result = n[..., np.newaxis]*p
3317 return self._checkresult(result, npcond, np.nan)
3319 def cov(self, n, p):
3320 """Covariance matrix of the multinomial distribution.
3322 Parameters
3323 ----------
3324 %(_doc_default_callparams)s
3326 Returns
3327 -------
3328 cov : ndarray
3329 The covariance matrix of the distribution
3330 """
3331 n, p, npcond = self._process_parameters(n, p)
3333 nn = n[..., np.newaxis, np.newaxis]
3334 result = nn * np.einsum('...j,...k->...jk', -p, p)
3336 # change the diagonal
3337 for i in range(p.shape[-1]):
3338 result[..., i, i] += n*p[..., i]
3340 return self._checkresult(result, npcond, np.nan)
3342 def entropy(self, n, p):
3343 r"""Compute the entropy of the multinomial distribution.
3345 The entropy is computed using this expression:
3347 .. math::
3349 f(x) = - \log n! - n\sum_{i=1}^k p_i \log p_i +
3350 \sum_{i=1}^k \sum_{x=0}^n \binom n x p_i^x(1-p_i)^{n-x} \log x!
3352 Parameters
3353 ----------
3354 %(_doc_default_callparams)s
3356 Returns
3357 -------
3358 h : scalar
3359 Entropy of the Multinomial distribution
3361 Notes
3362 -----
3363 %(_doc_callparams_note)s
3364 """
3365 n, p, npcond = self._process_parameters(n, p)
3367 x = np.r_[1:np.max(n)+1]
3369 term1 = n*np.sum(entr(p), axis=-1)
3370 term1 -= gammaln(n+1)
3372 n = n[..., np.newaxis]
3373 new_axes_needed = max(p.ndim, n.ndim) - x.ndim + 1
3374 x.shape += (1,)*new_axes_needed
3376 term2 = np.sum(binom.pmf(x, n, p)*gammaln(x+1),
3377 axis=(-1, -1-new_axes_needed))
3379 return self._checkresult(term1 + term2, npcond, np.nan)
3381 def rvs(self, n, p, size=None, random_state=None):
3382 """Draw random samples from a Multinomial distribution.
3384 Parameters
3385 ----------
3386 %(_doc_default_callparams)s
3387 size : integer or iterable of integers, optional
3388 Number of samples to draw (default 1).
3389 %(_doc_random_state)s
3391 Returns
3392 -------
3393 rvs : ndarray or scalar
3394 Random variates of shape (`size`, `len(p)`)
3396 Notes
3397 -----
3398 %(_doc_callparams_note)s
3399 """
3400 n, p, npcond = self._process_parameters(n, p)
3401 random_state = self._get_random_state(random_state)
3402 return random_state.multinomial(n, p, size)
3405multinomial = multinomial_gen()
3408class multinomial_frozen(multi_rv_frozen):
3409 r"""Create a frozen Multinomial distribution.
3411 Parameters
3412 ----------
3413 n : int
3414 number of trials
3415 p: array_like
3416 probability of a trial falling into each category; should sum to 1
3417 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
3418 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
3419 singleton is used.
3420 If `seed` is an int, a new ``RandomState`` instance is used,
3421 seeded with `seed`.
3422 If `seed` is already a ``Generator`` or ``RandomState`` instance then
3423 that instance is used.
3424 """
3425 def __init__(self, n, p, seed=None):
3426 self._dist = multinomial_gen(seed)
3427 self.n, self.p, self.npcond = self._dist._process_parameters(n, p)
3429 # monkey patch self._dist
3430 def _process_parameters(n, p):
3431 return self.n, self.p, self.npcond
3433 self._dist._process_parameters = _process_parameters
3435 def logpmf(self, x):
3436 return self._dist.logpmf(x, self.n, self.p)
3438 def pmf(self, x):
3439 return self._dist.pmf(x, self.n, self.p)
3441 def mean(self):
3442 return self._dist.mean(self.n, self.p)
3444 def cov(self):
3445 return self._dist.cov(self.n, self.p)
3447 def entropy(self):
3448 return self._dist.entropy(self.n, self.p)
3450 def rvs(self, size=1, random_state=None):
3451 return self._dist.rvs(self.n, self.p, size, random_state)
3454# Set frozen generator docstrings from corresponding docstrings in
3455# multinomial and fill in default strings in class docstrings
3456for name in ['logpmf', 'pmf', 'mean', 'cov', 'rvs']:
3457 method = multinomial_gen.__dict__[name]
3458 method_frozen = multinomial_frozen.__dict__[name]
3459 method_frozen.__doc__ = doccer.docformat(
3460 method.__doc__, multinomial_docdict_noparams)
3461 method.__doc__ = doccer.docformat(method.__doc__,
3462 multinomial_docdict_params)
3465class special_ortho_group_gen(multi_rv_generic):
3466 r"""A Special Orthogonal matrix (SO(N)) random variable.
3468 Return a random rotation matrix, drawn from the Haar distribution
3469 (the only uniform distribution on SO(N)) with a determinant of +1.
3471 The `dim` keyword specifies the dimension N.
3473 Methods
3474 -------
3475 rvs(dim=None, size=1, random_state=None)
3476 Draw random samples from SO(N).
3478 Parameters
3479 ----------
3480 dim : scalar
3481 Dimension of matrices
3482 seed : {None, int, np.random.RandomState, np.random.Generator}, optional
3483 Used for drawing random variates.
3484 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
3485 If `seed` is an int, a new ``RandomState`` instance is used, seeded
3486 with seed.
3487 If `seed` is already a ``RandomState`` or ``Generator`` instance,
3488 then that object is used.
3489 Default is `None`.
3491 Notes
3492 -----
3493 This class is wrapping the random_rot code from the MDP Toolkit,
3494 https://github.com/mdp-toolkit/mdp-toolkit
3496 Return a random rotation matrix, drawn from the Haar distribution
3497 (the only uniform distribution on SO(N)).
3498 The algorithm is described in the paper
3499 Stewart, G.W., "The efficient generation of random orthogonal
3500 matrices with an application to condition estimators", SIAM Journal
3501 on Numerical Analysis, 17(3), pp. 403-409, 1980.
3502 For more information see
3503 https://en.wikipedia.org/wiki/Orthogonal_matrix#Randomization
3505 See also the similar `ortho_group`. For a random rotation in three
3506 dimensions, see `scipy.spatial.transform.Rotation.random`.
3508 Examples
3509 --------
3510 >>> import numpy as np
3511 >>> from scipy.stats import special_ortho_group
3512 >>> x = special_ortho_group.rvs(3)
3514 >>> np.dot(x, x.T)
3515 array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16],
3516 [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16],
3517 [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]])
3519 >>> import scipy.linalg
3520 >>> scipy.linalg.det(x)
3521 1.0
3523 This generates one random matrix from SO(3). It is orthogonal and
3524 has a determinant of 1.
3526 Alternatively, the object may be called (as a function) to fix the `dim`
3527 parameter, returning a "frozen" special_ortho_group random variable:
3529 >>> rv = special_ortho_group(5)
3530 >>> # Frozen object with the same methods but holding the
3531 >>> # dimension parameter fixed.
3533 See Also
3534 --------
3535 ortho_group, scipy.spatial.transform.Rotation.random
3537 """
3539 def __init__(self, seed=None):
3540 super().__init__(seed)
3541 self.__doc__ = doccer.docformat(self.__doc__)
3543 def __call__(self, dim=None, seed=None):
3544 """Create a frozen SO(N) distribution.
3546 See `special_ortho_group_frozen` for more information.
3547 """
3548 return special_ortho_group_frozen(dim, seed=seed)
3550 def _process_parameters(self, dim):
3551 """Dimension N must be specified; it cannot be inferred."""
3552 if dim is None or not np.isscalar(dim) or dim <= 1 or dim != int(dim):
3553 raise ValueError("""Dimension of rotation must be specified,
3554 and must be a scalar greater than 1.""")
3556 return dim
3558 def rvs(self, dim, size=1, random_state=None):
3559 """Draw random samples from SO(N).
3561 Parameters
3562 ----------
3563 dim : integer
3564 Dimension of rotation space (N).
3565 size : integer, optional
3566 Number of samples to draw (default 1).
3568 Returns
3569 -------
3570 rvs : ndarray or scalar
3571 Random size N-dimensional matrices, dimension (size, dim, dim)
3573 """
3574 random_state = self._get_random_state(random_state)
3576 size = int(size)
3577 size = (size,) if size > 1 else ()
3579 dim = self._process_parameters(dim)
3581 # H represents a (dim, dim) matrix, while D represents the diagonal of
3582 # a (dim, dim) diagonal matrix. The algorithm that follows is
3583 # broadcasted on the leading shape in `size` to vectorize along
3584 # samples.
3585 H = np.empty(size + (dim, dim))
3586 H[..., :, :] = np.eye(dim)
3587 D = np.empty(size + (dim,))
3589 for n in range(dim-1):
3591 # x is a vector with length dim-n, xrow and xcol are views of it as
3592 # a row vector and column vector respectively. It's important they
3593 # are views and not copies because we are going to modify x
3594 # in-place.
3595 x = random_state.normal(size=size + (dim-n,))
3596 xrow = x[..., None, :]
3597 xcol = x[..., :, None]
3599 # This is the squared norm of x, without vectorization it would be
3600 # dot(x, x), to have proper broadcasting we use matmul and squeeze
3601 # out (convert to scalar) the resulting 1x1 matrix
3602 norm2 = np.matmul(xrow, xcol).squeeze((-2, -1))
3604 x0 = x[..., 0].copy()
3605 D[..., n] = np.where(x0 != 0, np.sign(x0), 1)
3606 x[..., 0] += D[..., n]*np.sqrt(norm2)
3608 # In renormalizing x we have to append an additional axis with
3609 # [..., None] to broadcast the scalar against the vector x
3610 x /= np.sqrt((norm2 - x0**2 + x[..., 0]**2) / 2.)[..., None]
3612 # Householder transformation, without vectorization the RHS can be
3613 # written as outer(H @ x, x) (apart from the slicing)
3614 H[..., :, n:] -= np.matmul(H[..., :, n:], xcol) * xrow
3616 D[..., -1] = (-1)**(dim-1)*D[..., :-1].prod(axis=-1)
3618 # Without vectorization this could be written as H = diag(D) @ H,
3619 # left-multiplication by a diagonal matrix amounts to multiplying each
3620 # row of H by an element of the diagonal, so we add a dummy axis for
3621 # the column index
3622 H *= D[..., :, None]
3623 return H
3626special_ortho_group = special_ortho_group_gen()
3629class special_ortho_group_frozen(multi_rv_frozen):
3630 def __init__(self, dim=None, seed=None):
3631 """Create a frozen SO(N) distribution.
3633 Parameters
3634 ----------
3635 dim : scalar
3636 Dimension of matrices
3637 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
3638 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
3639 singleton is used.
3640 If `seed` is an int, a new ``RandomState`` instance is used,
3641 seeded with `seed`.
3642 If `seed` is already a ``Generator`` or ``RandomState`` instance
3643 then that instance is used.
3645 Examples
3646 --------
3647 >>> from scipy.stats import special_ortho_group
3648 >>> g = special_ortho_group(5)
3649 >>> x = g.rvs()
3651 """
3652 self._dist = special_ortho_group_gen(seed)
3653 self.dim = self._dist._process_parameters(dim)
3655 def rvs(self, size=1, random_state=None):
3656 return self._dist.rvs(self.dim, size, random_state)
3659class ortho_group_gen(multi_rv_generic):
3660 r"""An Orthogonal matrix (O(N)) random variable.
3662 Return a random orthogonal matrix, drawn from the O(N) Haar
3663 distribution (the only uniform distribution on O(N)).
3665 The `dim` keyword specifies the dimension N.
3667 Methods
3668 -------
3669 rvs(dim=None, size=1, random_state=None)
3670 Draw random samples from O(N).
3672 Parameters
3673 ----------
3674 dim : scalar
3675 Dimension of matrices
3676 seed : {None, int, np.random.RandomState, np.random.Generator}, optional
3677 Used for drawing random variates.
3678 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
3679 If `seed` is an int, a new ``RandomState`` instance is used, seeded
3680 with seed.
3681 If `seed` is already a ``RandomState`` or ``Generator`` instance,
3682 then that object is used.
3683 Default is `None`.
3685 Notes
3686 -----
3687 This class is closely related to `special_ortho_group`.
3689 Some care is taken to avoid numerical error, as per the paper by Mezzadri.
3691 References
3692 ----------
3693 .. [1] F. Mezzadri, "How to generate random matrices from the classical
3694 compact groups", :arXiv:`math-ph/0609050v2`.
3696 Examples
3697 --------
3698 >>> import numpy as np
3699 >>> from scipy.stats import ortho_group
3700 >>> x = ortho_group.rvs(3)
3702 >>> np.dot(x, x.T)
3703 array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16],
3704 [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16],
3705 [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]])
3707 >>> import scipy.linalg
3708 >>> np.fabs(scipy.linalg.det(x))
3709 1.0
3711 This generates one random matrix from O(3). It is orthogonal and
3712 has a determinant of +1 or -1.
3714 Alternatively, the object may be called (as a function) to fix the `dim`
3715 parameter, returning a "frozen" ortho_group random variable:
3717 >>> rv = ortho_group(5)
3718 >>> # Frozen object with the same methods but holding the
3719 >>> # dimension parameter fixed.
3721 See Also
3722 --------
3723 special_ortho_group
3724 """
3726 def __init__(self, seed=None):
3727 super().__init__(seed)
3728 self.__doc__ = doccer.docformat(self.__doc__)
3730 def __call__(self, dim=None, seed=None):
3731 """Create a frozen O(N) distribution.
3733 See `ortho_group_frozen` for more information.
3734 """
3735 return ortho_group_frozen(dim, seed=seed)
3737 def _process_parameters(self, dim):
3738 """Dimension N must be specified; it cannot be inferred."""
3739 if dim is None or not np.isscalar(dim) or dim <= 1 or dim != int(dim):
3740 raise ValueError("Dimension of rotation must be specified,"
3741 "and must be a scalar greater than 1.")
3743 return dim
3745 def rvs(self, dim, size=1, random_state=None):
3746 """Draw random samples from O(N).
3748 Parameters
3749 ----------
3750 dim : integer
3751 Dimension of rotation space (N).
3752 size : integer, optional
3753 Number of samples to draw (default 1).
3755 Returns
3756 -------
3757 rvs : ndarray or scalar
3758 Random size N-dimensional matrices, dimension (size, dim, dim)
3760 """
3761 random_state = self._get_random_state(random_state)
3763 size = int(size)
3764 if size > 1 and NumpyVersion(np.__version__) < '1.22.0':
3765 return np.array([self.rvs(dim, size=1, random_state=random_state)
3766 for i in range(size)])
3768 dim = self._process_parameters(dim)
3770 size = (size,) if size > 1 else ()
3771 z = random_state.normal(size=size + (dim, dim))
3772 q, r = np.linalg.qr(z)
3773 # The last two dimensions are the rows and columns of R matrices.
3774 # Extract the diagonals. Note that this eliminates a dimension.
3775 d = r.diagonal(offset=0, axis1=-2, axis2=-1)
3776 # Add back a dimension for proper broadcasting: we're dividing
3777 # each row of each R matrix by the diagonal of the R matrix.
3778 q *= (d/abs(d))[..., np.newaxis, :] # to broadcast properly
3779 return q
3782ortho_group = ortho_group_gen()
3785class ortho_group_frozen(multi_rv_frozen):
3786 def __init__(self, dim=None, seed=None):
3787 """Create a frozen O(N) distribution.
3789 Parameters
3790 ----------
3791 dim : scalar
3792 Dimension of matrices
3793 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
3794 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
3795 singleton is used.
3796 If `seed` is an int, a new ``RandomState`` instance is used,
3797 seeded with `seed`.
3798 If `seed` is already a ``Generator`` or ``RandomState`` instance
3799 then that instance is used.
3801 Examples
3802 --------
3803 >>> from scipy.stats import ortho_group
3804 >>> g = ortho_group(5)
3805 >>> x = g.rvs()
3807 """
3808 self._dist = ortho_group_gen(seed)
3809 self.dim = self._dist._process_parameters(dim)
3811 def rvs(self, size=1, random_state=None):
3812 return self._dist.rvs(self.dim, size, random_state)
3815class random_correlation_gen(multi_rv_generic):
3816 r"""A random correlation matrix.
3818 Return a random correlation matrix, given a vector of eigenvalues.
3820 The `eigs` keyword specifies the eigenvalues of the correlation matrix,
3821 and implies the dimension.
3823 Methods
3824 -------
3825 rvs(eigs=None, random_state=None)
3826 Draw random correlation matrices, all with eigenvalues eigs.
3828 Parameters
3829 ----------
3830 eigs : 1d ndarray
3831 Eigenvalues of correlation matrix
3832 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
3833 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
3834 singleton is used.
3835 If `seed` is an int, a new ``RandomState`` instance is used,
3836 seeded with `seed`.
3837 If `seed` is already a ``Generator`` or ``RandomState`` instance
3838 then that instance is used.
3839 tol : float, optional
3840 Tolerance for input parameter checks
3841 diag_tol : float, optional
3842 Tolerance for deviation of the diagonal of the resulting
3843 matrix. Default: 1e-7
3845 Raises
3846 ------
3847 RuntimeError
3848 Floating point error prevented generating a valid correlation
3849 matrix.
3851 Returns
3852 -------
3853 rvs : ndarray or scalar
3854 Random size N-dimensional matrices, dimension (size, dim, dim),
3855 each having eigenvalues eigs.
3857 Notes
3858 -----
3860 Generates a random correlation matrix following a numerically stable
3861 algorithm spelled out by Davies & Higham. This algorithm uses a single O(N)
3862 similarity transformation to construct a symmetric positive semi-definite
3863 matrix, and applies a series of Givens rotations to scale it to have ones
3864 on the diagonal.
3866 References
3867 ----------
3869 .. [1] Davies, Philip I; Higham, Nicholas J; "Numerically stable generation
3870 of correlation matrices and their factors", BIT 2000, Vol. 40,
3871 No. 4, pp. 640 651
3873 Examples
3874 --------
3875 >>> import numpy as np
3876 >>> from scipy.stats import random_correlation
3877 >>> rng = np.random.default_rng()
3878 >>> x = random_correlation.rvs((.5, .8, 1.2, 1.5), random_state=rng)
3879 >>> x
3880 array([[ 1. , -0.02423399, 0.03130519, 0.4946965 ],
3881 [-0.02423399, 1. , 0.20334736, 0.04039817],
3882 [ 0.03130519, 0.20334736, 1. , 0.02694275],
3883 [ 0.4946965 , 0.04039817, 0.02694275, 1. ]])
3884 >>> import scipy.linalg
3885 >>> e, v = scipy.linalg.eigh(x)
3886 >>> e
3887 array([ 0.5, 0.8, 1.2, 1.5])
3889 """
3891 def __init__(self, seed=None):
3892 super().__init__(seed)
3893 self.__doc__ = doccer.docformat(self.__doc__)
3895 def __call__(self, eigs, seed=None, tol=1e-13, diag_tol=1e-7):
3896 """Create a frozen random correlation matrix.
3898 See `random_correlation_frozen` for more information.
3899 """
3900 return random_correlation_frozen(eigs, seed=seed, tol=tol,
3901 diag_tol=diag_tol)
3903 def _process_parameters(self, eigs, tol):
3904 eigs = np.asarray(eigs, dtype=float)
3905 dim = eigs.size
3907 if eigs.ndim != 1 or eigs.shape[0] != dim or dim <= 1:
3908 raise ValueError("Array 'eigs' must be a vector of length "
3909 "greater than 1.")
3911 if np.fabs(np.sum(eigs) - dim) > tol:
3912 raise ValueError("Sum of eigenvalues must equal dimensionality.")
3914 for x in eigs:
3915 if x < -tol:
3916 raise ValueError("All eigenvalues must be non-negative.")
3918 return dim, eigs
3920 def _givens_to_1(self, aii, ajj, aij):
3921 """Computes a 2x2 Givens matrix to put 1's on the diagonal.
3923 The input matrix is a 2x2 symmetric matrix M = [ aii aij ; aij ajj ].
3925 The output matrix g is a 2x2 anti-symmetric matrix of the form
3926 [ c s ; -s c ]; the elements c and s are returned.
3928 Applying the output matrix to the input matrix (as b=g.T M g)
3929 results in a matrix with bii=1, provided tr(M) - det(M) >= 1
3930 and floating point issues do not occur. Otherwise, some other
3931 valid rotation is returned. When tr(M)==2, also bjj=1.
3933 """
3934 aiid = aii - 1.
3935 ajjd = ajj - 1.
3937 if ajjd == 0:
3938 # ajj==1, so swap aii and ajj to avoid division by zero
3939 return 0., 1.
3941 dd = math.sqrt(max(aij**2 - aiid*ajjd, 0))
3943 # The choice of t should be chosen to avoid cancellation [1]
3944 t = (aij + math.copysign(dd, aij)) / ajjd
3945 c = 1. / math.sqrt(1. + t*t)
3946 if c == 0:
3947 # Underflow
3948 s = 1.0
3949 else:
3950 s = c*t
3951 return c, s
3953 def _to_corr(self, m):
3954 """
3955 Given a psd matrix m, rotate to put one's on the diagonal, turning it
3956 into a correlation matrix. This also requires the trace equal the
3957 dimensionality. Note: modifies input matrix
3958 """
3959 # Check requirements for in-place Givens
3960 if not (m.flags.c_contiguous and m.dtype == np.float64 and
3961 m.shape[0] == m.shape[1]):
3962 raise ValueError()
3964 d = m.shape[0]
3965 for i in range(d-1):
3966 if m[i, i] == 1:
3967 continue
3968 elif m[i, i] > 1:
3969 for j in range(i+1, d):
3970 if m[j, j] < 1:
3971 break
3972 else:
3973 for j in range(i+1, d):
3974 if m[j, j] > 1:
3975 break
3977 c, s = self._givens_to_1(m[i, i], m[j, j], m[i, j])
3979 # Use BLAS to apply Givens rotations in-place. Equivalent to:
3980 # g = np.eye(d)
3981 # g[i, i] = g[j,j] = c
3982 # g[j, i] = -s; g[i, j] = s
3983 # m = np.dot(g.T, np.dot(m, g))
3984 mv = m.ravel()
3985 drot(mv, mv, c, -s, n=d,
3986 offx=i*d, incx=1, offy=j*d, incy=1,
3987 overwrite_x=True, overwrite_y=True)
3988 drot(mv, mv, c, -s, n=d,
3989 offx=i, incx=d, offy=j, incy=d,
3990 overwrite_x=True, overwrite_y=True)
3992 return m
3994 def rvs(self, eigs, random_state=None, tol=1e-13, diag_tol=1e-7):
3995 """Draw random correlation matrices.
3997 Parameters
3998 ----------
3999 eigs : 1d ndarray
4000 Eigenvalues of correlation matrix
4001 tol : float, optional
4002 Tolerance for input parameter checks
4003 diag_tol : float, optional
4004 Tolerance for deviation of the diagonal of the resulting
4005 matrix. Default: 1e-7
4007 Raises
4008 ------
4009 RuntimeError
4010 Floating point error prevented generating a valid correlation
4011 matrix.
4013 Returns
4014 -------
4015 rvs : ndarray or scalar
4016 Random size N-dimensional matrices, dimension (size, dim, dim),
4017 each having eigenvalues eigs.
4019 """
4020 dim, eigs = self._process_parameters(eigs, tol=tol)
4022 random_state = self._get_random_state(random_state)
4024 m = ortho_group.rvs(dim, random_state=random_state)
4025 m = np.dot(np.dot(m, np.diag(eigs)), m.T) # Set the trace of m
4026 m = self._to_corr(m) # Carefully rotate to unit diagonal
4028 # Check diagonal
4029 if abs(m.diagonal() - 1).max() > diag_tol:
4030 raise RuntimeError("Failed to generate a valid correlation matrix")
4032 return m
4035random_correlation = random_correlation_gen()
4038class random_correlation_frozen(multi_rv_frozen):
4039 def __init__(self, eigs, seed=None, tol=1e-13, diag_tol=1e-7):
4040 """Create a frozen random correlation matrix distribution.
4042 Parameters
4043 ----------
4044 eigs : 1d ndarray
4045 Eigenvalues of correlation matrix
4046 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
4047 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
4048 singleton is used.
4049 If `seed` is an int, a new ``RandomState`` instance is used,
4050 seeded with `seed`.
4051 If `seed` is already a ``Generator`` or ``RandomState`` instance
4052 then that instance is used.
4053 tol : float, optional
4054 Tolerance for input parameter checks
4055 diag_tol : float, optional
4056 Tolerance for deviation of the diagonal of the resulting
4057 matrix. Default: 1e-7
4059 Raises
4060 ------
4061 RuntimeError
4062 Floating point error prevented generating a valid correlation
4063 matrix.
4065 Returns
4066 -------
4067 rvs : ndarray or scalar
4068 Random size N-dimensional matrices, dimension (size, dim, dim),
4069 each having eigenvalues eigs.
4070 """
4072 self._dist = random_correlation_gen(seed)
4073 self.tol = tol
4074 self.diag_tol = diag_tol
4075 _, self.eigs = self._dist._process_parameters(eigs, tol=self.tol)
4077 def rvs(self, random_state=None):
4078 return self._dist.rvs(self.eigs, random_state=random_state,
4079 tol=self.tol, diag_tol=self.diag_tol)
4082class unitary_group_gen(multi_rv_generic):
4083 r"""A matrix-valued U(N) random variable.
4085 Return a random unitary matrix.
4087 The `dim` keyword specifies the dimension N.
4089 Methods
4090 -------
4091 rvs(dim=None, size=1, random_state=None)
4092 Draw random samples from U(N).
4094 Parameters
4095 ----------
4096 dim : scalar
4097 Dimension of matrices
4098 seed : {None, int, np.random.RandomState, np.random.Generator}, optional
4099 Used for drawing random variates.
4100 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
4101 If `seed` is an int, a new ``RandomState`` instance is used, seeded
4102 with seed.
4103 If `seed` is already a ``RandomState`` or ``Generator`` instance,
4104 then that object is used.
4105 Default is `None`.
4107 Notes
4108 -----
4109 This class is similar to `ortho_group`.
4111 References
4112 ----------
4113 .. [1] F. Mezzadri, "How to generate random matrices from the classical
4114 compact groups", :arXiv:`math-ph/0609050v2`.
4116 Examples
4117 --------
4118 >>> import numpy as np
4119 >>> from scipy.stats import unitary_group
4120 >>> x = unitary_group.rvs(3)
4122 >>> np.dot(x, x.conj().T)
4123 array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16],
4124 [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16],
4125 [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]])
4127 This generates one random matrix from U(3). The dot product confirms that
4128 it is unitary up to machine precision.
4130 Alternatively, the object may be called (as a function) to fix the `dim`
4131 parameter, return a "frozen" unitary_group random variable:
4133 >>> rv = unitary_group(5)
4135 See Also
4136 --------
4137 ortho_group
4139 """
4141 def __init__(self, seed=None):
4142 super().__init__(seed)
4143 self.__doc__ = doccer.docformat(self.__doc__)
4145 def __call__(self, dim=None, seed=None):
4146 """Create a frozen (U(N)) n-dimensional unitary matrix distribution.
4148 See `unitary_group_frozen` for more information.
4149 """
4150 return unitary_group_frozen(dim, seed=seed)
4152 def _process_parameters(self, dim):
4153 """Dimension N must be specified; it cannot be inferred."""
4154 if dim is None or not np.isscalar(dim) or dim <= 1 or dim != int(dim):
4155 raise ValueError("Dimension of rotation must be specified,"
4156 "and must be a scalar greater than 1.")
4158 return dim
4160 def rvs(self, dim, size=1, random_state=None):
4161 """Draw random samples from U(N).
4163 Parameters
4164 ----------
4165 dim : integer
4166 Dimension of space (N).
4167 size : integer, optional
4168 Number of samples to draw (default 1).
4170 Returns
4171 -------
4172 rvs : ndarray or scalar
4173 Random size N-dimensional matrices, dimension (size, dim, dim)
4175 """
4176 random_state = self._get_random_state(random_state)
4178 size = int(size)
4179 if size > 1 and NumpyVersion(np.__version__) < '1.22.0':
4180 return np.array([self.rvs(dim, size=1, random_state=random_state)
4181 for i in range(size)])
4183 dim = self._process_parameters(dim)
4185 size = (size,) if size > 1 else ()
4186 z = 1/math.sqrt(2)*(random_state.normal(size=size + (dim, dim)) +
4187 1j*random_state.normal(size=size + (dim, dim)))
4188 q, r = np.linalg.qr(z)
4189 # The last two dimensions are the rows and columns of R matrices.
4190 # Extract the diagonals. Note that this eliminates a dimension.
4191 d = r.diagonal(offset=0, axis1=-2, axis2=-1)
4192 # Add back a dimension for proper broadcasting: we're dividing
4193 # each row of each R matrix by the diagonal of the R matrix.
4194 q *= (d/abs(d))[..., np.newaxis, :] # to broadcast properly
4195 return q
4198unitary_group = unitary_group_gen()
4201class unitary_group_frozen(multi_rv_frozen):
4202 def __init__(self, dim=None, seed=None):
4203 """Create a frozen (U(N)) n-dimensional unitary matrix distribution.
4205 Parameters
4206 ----------
4207 dim : scalar
4208 Dimension of matrices
4209 seed : {None, int, `numpy.random.Generator`, `numpy.random.RandomState`}, optional
4210 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
4211 singleton is used.
4212 If `seed` is an int, a new ``RandomState`` instance is used,
4213 seeded with `seed`.
4214 If `seed` is already a ``Generator`` or ``RandomState`` instance
4215 then that instance is used.
4217 Examples
4218 --------
4219 >>> from scipy.stats import unitary_group
4220 >>> x = unitary_group(3)
4221 >>> x.rvs()
4223 """
4224 self._dist = unitary_group_gen(seed)
4225 self.dim = self._dist._process_parameters(dim)
4227 def rvs(self, size=1, random_state=None):
4228 return self._dist.rvs(self.dim, size, random_state)
4231_mvt_doc_default_callparams = """\
4232loc : array_like, optional
4233 Location of the distribution. (default ``0``)
4234shape : array_like, optional
4235 Positive semidefinite matrix of the distribution. (default ``1``)
4236df : float, optional
4237 Degrees of freedom of the distribution; must be greater than zero.
4238 If ``np.inf`` then results are multivariate normal. The default is ``1``.
4239allow_singular : bool, optional
4240 Whether to allow a singular matrix. (default ``False``)
4241"""
4243_mvt_doc_callparams_note = """\
4244Setting the parameter `loc` to ``None`` is equivalent to having `loc`
4245be the zero-vector. The parameter `shape` can be a scalar, in which case
4246the shape matrix is the identity times that value, a vector of
4247diagonal entries for the shape matrix, or a two-dimensional array_like.
4248"""
4250_mvt_doc_frozen_callparams_note = """\
4251See class definition for a detailed description of parameters."""
4253mvt_docdict_params = {
4254 '_mvt_doc_default_callparams': _mvt_doc_default_callparams,
4255 '_mvt_doc_callparams_note': _mvt_doc_callparams_note,
4256 '_doc_random_state': _doc_random_state
4257}
4259mvt_docdict_noparams = {
4260 '_mvt_doc_default_callparams': "",
4261 '_mvt_doc_callparams_note': _mvt_doc_frozen_callparams_note,
4262 '_doc_random_state': _doc_random_state
4263}
4266class multivariate_t_gen(multi_rv_generic):
4267 r"""A multivariate t-distributed random variable.
4269 The `loc` parameter specifies the location. The `shape` parameter specifies
4270 the positive semidefinite shape matrix. The `df` parameter specifies the
4271 degrees of freedom.
4273 In addition to calling the methods below, the object itself may be called
4274 as a function to fix the location, shape matrix, and degrees of freedom
4275 parameters, returning a "frozen" multivariate t-distribution random.
4277 Methods
4278 -------
4279 pdf(x, loc=None, shape=1, df=1, allow_singular=False)
4280 Probability density function.
4281 logpdf(x, loc=None, shape=1, df=1, allow_singular=False)
4282 Log of the probability density function.
4283 cdf(x, loc=None, shape=1, df=1, allow_singular=False, *,
4284 maxpts=None, lower_limit=None, random_state=None)
4285 Cumulative distribution function.
4286 rvs(loc=None, shape=1, df=1, size=1, random_state=None)
4287 Draw random samples from a multivariate t-distribution.
4288 entropy(loc=None, shape=1, df=1)
4289 Differential entropy of a multivariate t-distribution.
4291 Parameters
4292 ----------
4293 %(_mvt_doc_default_callparams)s
4294 %(_doc_random_state)s
4296 Notes
4297 -----
4298 %(_mvt_doc_callparams_note)s
4299 The matrix `shape` must be a (symmetric) positive semidefinite matrix. The
4300 determinant and inverse of `shape` are computed as the pseudo-determinant
4301 and pseudo-inverse, respectively, so that `shape` does not need to have
4302 full rank.
4304 The probability density function for `multivariate_t` is
4306 .. math::
4308 f(x) = \frac{\Gamma((\nu + p)/2)}{\Gamma(\nu/2)\nu^{p/2}\pi^{p/2}|\Sigma|^{1/2}}
4309 \left[1 + \frac{1}{\nu} (\mathbf{x} - \boldsymbol{\mu})^{\top}
4310 \boldsymbol{\Sigma}^{-1}
4311 (\mathbf{x} - \boldsymbol{\mu}) \right]^{-(\nu + p)/2},
4313 where :math:`p` is the dimension of :math:`\mathbf{x}`,
4314 :math:`\boldsymbol{\mu}` is the :math:`p`-dimensional location,
4315 :math:`\boldsymbol{\Sigma}` the :math:`p \times p`-dimensional shape
4316 matrix, and :math:`\nu` is the degrees of freedom.
4318 .. versionadded:: 1.6.0
4320 References
4321 ----------
4322 [1] Arellano-Valle et al. "Shannon Entropy and Mutual Information for
4323 Multivariate Skew-Elliptical Distributions". Scandinavian Journal
4324 of Statistics. Vol. 40, issue 1.
4326 Examples
4327 --------
4328 The object may be called (as a function) to fix the `loc`, `shape`,
4329 `df`, and `allow_singular` parameters, returning a "frozen"
4330 multivariate_t random variable:
4332 >>> import numpy as np
4333 >>> from scipy.stats import multivariate_t
4334 >>> rv = multivariate_t([1.0, -0.5], [[2.1, 0.3], [0.3, 1.5]], df=2)
4335 >>> # Frozen object with the same methods but holding the given location,
4336 >>> # scale, and degrees of freedom fixed.
4338 Create a contour plot of the PDF.
4340 >>> import matplotlib.pyplot as plt
4341 >>> x, y = np.mgrid[-1:3:.01, -2:1.5:.01]
4342 >>> pos = np.dstack((x, y))
4343 >>> fig, ax = plt.subplots(1, 1)
4344 >>> ax.set_aspect('equal')
4345 >>> plt.contourf(x, y, rv.pdf(pos))
4347 """
4349 def __init__(self, seed=None):
4350 """Initialize a multivariate t-distributed random variable.
4352 Parameters
4353 ----------
4354 seed : Random state.
4356 """
4357 super().__init__(seed)
4358 self.__doc__ = doccer.docformat(self.__doc__, mvt_docdict_params)
4359 self._random_state = check_random_state(seed)
4361 def __call__(self, loc=None, shape=1, df=1, allow_singular=False,
4362 seed=None):
4363 """Create a frozen multivariate t-distribution.
4365 See `multivariate_t_frozen` for parameters.
4366 """
4367 if df == np.inf:
4368 return multivariate_normal_frozen(mean=loc, cov=shape,
4369 allow_singular=allow_singular,
4370 seed=seed)
4371 return multivariate_t_frozen(loc=loc, shape=shape, df=df,
4372 allow_singular=allow_singular, seed=seed)
4374 def pdf(self, x, loc=None, shape=1, df=1, allow_singular=False):
4375 """Multivariate t-distribution probability density function.
4377 Parameters
4378 ----------
4379 x : array_like
4380 Points at which to evaluate the probability density function.
4381 %(_mvt_doc_default_callparams)s
4383 Returns
4384 -------
4385 pdf : Probability density function evaluated at `x`.
4387 Examples
4388 --------
4389 >>> from scipy.stats import multivariate_t
4390 >>> x = [0.4, 5]
4391 >>> loc = [0, 1]
4392 >>> shape = [[1, 0.1], [0.1, 1]]
4393 >>> df = 7
4394 >>> multivariate_t.pdf(x, loc, shape, df)
4395 0.00075713
4397 """
4398 dim, loc, shape, df = self._process_parameters(loc, shape, df)
4399 x = self._process_quantiles(x, dim)
4400 shape_info = _PSD(shape, allow_singular=allow_singular)
4401 logpdf = self._logpdf(x, loc, shape_info.U, shape_info.log_pdet, df,
4402 dim, shape_info.rank)
4403 return np.exp(logpdf)
4405 def logpdf(self, x, loc=None, shape=1, df=1):
4406 """Log of the multivariate t-distribution probability density function.
4408 Parameters
4409 ----------
4410 x : array_like
4411 Points at which to evaluate the log of the probability density
4412 function.
4413 %(_mvt_doc_default_callparams)s
4415 Returns
4416 -------
4417 logpdf : Log of the probability density function evaluated at `x`.
4419 Examples
4420 --------
4421 >>> from scipy.stats import multivariate_t
4422 >>> x = [0.4, 5]
4423 >>> loc = [0, 1]
4424 >>> shape = [[1, 0.1], [0.1, 1]]
4425 >>> df = 7
4426 >>> multivariate_t.logpdf(x, loc, shape, df)
4427 -7.1859802
4429 See Also
4430 --------
4431 pdf : Probability density function.
4433 """
4434 dim, loc, shape, df = self._process_parameters(loc, shape, df)
4435 x = self._process_quantiles(x, dim)
4436 shape_info = _PSD(shape)
4437 return self._logpdf(x, loc, shape_info.U, shape_info.log_pdet, df, dim,
4438 shape_info.rank)
4440 def _logpdf(self, x, loc, prec_U, log_pdet, df, dim, rank):
4441 """Utility method `pdf`, `logpdf` for parameters.
4443 Parameters
4444 ----------
4445 x : ndarray
4446 Points at which to evaluate the log of the probability density
4447 function.
4448 loc : ndarray
4449 Location of the distribution.
4450 prec_U : ndarray
4451 A decomposition such that `np.dot(prec_U, prec_U.T)` is the inverse
4452 of the shape matrix.
4453 log_pdet : float
4454 Logarithm of the determinant of the shape matrix.
4455 df : float
4456 Degrees of freedom of the distribution.
4457 dim : int
4458 Dimension of the quantiles x.
4459 rank : int
4460 Rank of the shape matrix.
4462 Notes
4463 -----
4464 As this function does no argument checking, it should not be called
4465 directly; use 'logpdf' instead.
4467 """
4468 if df == np.inf:
4469 return multivariate_normal._logpdf(x, loc, prec_U, log_pdet, rank)
4471 dev = x - loc
4472 maha = np.square(np.dot(dev, prec_U)).sum(axis=-1)
4474 t = 0.5 * (df + dim)
4475 A = gammaln(t)
4476 B = gammaln(0.5 * df)
4477 C = dim/2. * np.log(df * np.pi)
4478 D = 0.5 * log_pdet
4479 E = -t * np.log(1 + (1./df) * maha)
4481 return _squeeze_output(A - B - C - D + E)
4483 def _cdf(self, x, loc, shape, df, dim, maxpts=None, lower_limit=None,
4484 random_state=None):
4486 # All of this - random state validation, maxpts, apply_along_axis,
4487 # etc. needs to go in this private method unless we want
4488 # frozen distribution's `cdf` method to duplicate it or call `cdf`,
4489 # which would require re-processing parameters
4490 if random_state is not None:
4491 rng = check_random_state(random_state)
4492 else:
4493 rng = self._random_state
4495 if not maxpts:
4496 maxpts = 1000 * dim
4498 x = self._process_quantiles(x, dim)
4499 lower_limit = (np.full(loc.shape, -np.inf)
4500 if lower_limit is None else lower_limit)
4502 # remove the mean
4503 x, lower_limit = x - loc, lower_limit - loc
4505 b, a = np.broadcast_arrays(x, lower_limit)
4506 i_swap = b < a
4507 signs = (-1)**(i_swap.sum(axis=-1)) # odd # of swaps -> negative
4508 a, b = a.copy(), b.copy()
4509 a[i_swap], b[i_swap] = b[i_swap], a[i_swap]
4510 n = x.shape[-1]
4511 limits = np.concatenate((a, b), axis=-1)
4513 def func1d(limits):
4514 a, b = limits[:n], limits[n:]
4515 return _qmvt(maxpts, df, shape, a, b, rng)[0]
4517 res = np.apply_along_axis(func1d, -1, limits) * signs
4518 # Fixing the output shape for existing distributions is a separate
4519 # issue. For now, let's keep this consistent with pdf.
4520 return _squeeze_output(res)
4522 def cdf(self, x, loc=None, shape=1, df=1, allow_singular=False, *,
4523 maxpts=None, lower_limit=None, random_state=None):
4524 """Multivariate t-distribution cumulative distribution function.
4526 Parameters
4527 ----------
4528 x : array_like
4529 Points at which to evaluate the cumulative distribution function.
4530 %(_mvt_doc_default_callparams)s
4531 maxpts : int, optional
4532 Maximum number of points to use for integration. The default is
4533 1000 times the number of dimensions.
4534 lower_limit : array_like, optional
4535 Lower limit of integration of the cumulative distribution function.
4536 Default is negative infinity. Must be broadcastable with `x`.
4537 %(_doc_random_state)s
4539 Returns
4540 -------
4541 cdf : ndarray or scalar
4542 Cumulative distribution function evaluated at `x`.
4544 Examples
4545 --------
4546 >>> from scipy.stats import multivariate_t
4547 >>> x = [0.4, 5]
4548 >>> loc = [0, 1]
4549 >>> shape = [[1, 0.1], [0.1, 1]]
4550 >>> df = 7
4551 >>> multivariate_t.cdf(x, loc, shape, df)
4552 0.64798491
4554 """
4555 dim, loc, shape, df = self._process_parameters(loc, shape, df)
4556 shape = _PSD(shape, allow_singular=allow_singular)._M
4558 return self._cdf(x, loc, shape, df, dim, maxpts,
4559 lower_limit, random_state)
4561 def _entropy(self, dim, df=1, shape=1):
4562 if df == np.inf:
4563 return multivariate_normal(None, cov=shape).entropy()
4565 shape_info = _PSD(shape)
4566 halfsum = 0.5 * (dim + df)
4567 half_df = 0.5 * df
4568 return (-gammaln(halfsum) + gammaln(half_df)
4569 + 0.5 * dim * np.log(df * np.pi) + halfsum
4570 * (psi(halfsum) - psi(half_df))
4571 + 0.5 * shape_info.log_pdet)
4573 def entropy(self, loc=None, shape=1, df=1):
4574 """Calculate the differential entropy of a multivariate
4575 t-distribution.
4577 Parameters
4578 ----------
4579 %(_mvt_doc_default_callparams)s
4581 Returns
4582 -------
4583 h : float
4584 Differential entropy
4586 """
4587 dim, loc, shape, df = self._process_parameters(None, shape, df)
4588 return self._entropy(dim, df, shape)
4590 def rvs(self, loc=None, shape=1, df=1, size=1, random_state=None):
4591 """Draw random samples from a multivariate t-distribution.
4593 Parameters
4594 ----------
4595 %(_mvt_doc_default_callparams)s
4596 size : integer, optional
4597 Number of samples to draw (default 1).
4598 %(_doc_random_state)s
4600 Returns
4601 -------
4602 rvs : ndarray or scalar
4603 Random variates of size (`size`, `P`), where `P` is the
4604 dimension of the random variable.
4606 Examples
4607 --------
4608 >>> from scipy.stats import multivariate_t
4609 >>> x = [0.4, 5]
4610 >>> loc = [0, 1]
4611 >>> shape = [[1, 0.1], [0.1, 1]]
4612 >>> df = 7
4613 >>> multivariate_t.rvs(loc, shape, df)
4614 array([[0.93477495, 3.00408716]])
4616 """
4617 # For implementation details, see equation (3):
4618 #
4619 # Hofert, "On Sampling from the Multivariatet Distribution", 2013
4620 # http://rjournal.github.io/archive/2013-2/hofert.pdf
4621 #
4622 dim, loc, shape, df = self._process_parameters(loc, shape, df)
4623 if random_state is not None:
4624 rng = check_random_state(random_state)
4625 else:
4626 rng = self._random_state
4628 if np.isinf(df):
4629 x = np.ones(size)
4630 else:
4631 x = rng.chisquare(df, size=size) / df
4633 z = rng.multivariate_normal(np.zeros(dim), shape, size=size)
4634 samples = loc + z / np.sqrt(x)[..., None]
4635 return _squeeze_output(samples)
4637 def _process_quantiles(self, x, dim):
4638 """
4639 Adjust quantiles array so that last axis labels the components of
4640 each data point.
4641 """
4642 x = np.asarray(x, dtype=float)
4643 if x.ndim == 0:
4644 x = x[np.newaxis]
4645 elif x.ndim == 1:
4646 if dim == 1:
4647 x = x[:, np.newaxis]
4648 else:
4649 x = x[np.newaxis, :]
4650 return x
4652 def _process_parameters(self, loc, shape, df):
4653 """
4654 Infer dimensionality from location array and shape matrix, handle
4655 defaults, and ensure compatible dimensions.
4656 """
4657 if loc is None and shape is None:
4658 loc = np.asarray(0, dtype=float)
4659 shape = np.asarray(1, dtype=float)
4660 dim = 1
4661 elif loc is None:
4662 shape = np.asarray(shape, dtype=float)
4663 if shape.ndim < 2:
4664 dim = 1
4665 else:
4666 dim = shape.shape[0]
4667 loc = np.zeros(dim)
4668 elif shape is None:
4669 loc = np.asarray(loc, dtype=float)
4670 dim = loc.size
4671 shape = np.eye(dim)
4672 else:
4673 shape = np.asarray(shape, dtype=float)
4674 loc = np.asarray(loc, dtype=float)
4675 dim = loc.size
4677 if dim == 1:
4678 loc = loc.reshape(1)
4679 shape = shape.reshape(1, 1)
4681 if loc.ndim != 1 or loc.shape[0] != dim:
4682 raise ValueError("Array 'loc' must be a vector of length %d." %
4683 dim)
4684 if shape.ndim == 0:
4685 shape = shape * np.eye(dim)
4686 elif shape.ndim == 1:
4687 shape = np.diag(shape)
4688 elif shape.ndim == 2 and shape.shape != (dim, dim):
4689 rows, cols = shape.shape
4690 if rows != cols:
4691 msg = ("Array 'cov' must be square if it is two dimensional,"
4692 " but cov.shape = %s." % str(shape.shape))
4693 else:
4694 msg = ("Dimension mismatch: array 'cov' is of shape %s,"
4695 " but 'loc' is a vector of length %d.")
4696 msg = msg % (str(shape.shape), len(loc))
4697 raise ValueError(msg)
4698 elif shape.ndim > 2:
4699 raise ValueError("Array 'cov' must be at most two-dimensional,"
4700 " but cov.ndim = %d" % shape.ndim)
4702 # Process degrees of freedom.
4703 if df is None:
4704 df = 1
4705 elif df <= 0:
4706 raise ValueError("'df' must be greater than zero.")
4707 elif np.isnan(df):
4708 raise ValueError("'df' is 'nan' but must be greater than zero or 'np.inf'.")
4710 return dim, loc, shape, df
4713class multivariate_t_frozen(multi_rv_frozen):
4715 def __init__(self, loc=None, shape=1, df=1, allow_singular=False,
4716 seed=None):
4717 """Create a frozen multivariate t distribution.
4719 Parameters
4720 ----------
4721 %(_mvt_doc_default_callparams)s
4723 Examples
4724 --------
4725 >>> import numpy as np
4726 >>> loc = np.zeros(3)
4727 >>> shape = np.eye(3)
4728 >>> df = 10
4729 >>> dist = multivariate_t(loc, shape, df)
4730 >>> dist.rvs()
4731 array([[ 0.81412036, -1.53612361, 0.42199647]])
4732 >>> dist.pdf([1, 1, 1])
4733 array([0.01237803])
4735 """
4736 self._dist = multivariate_t_gen(seed)
4737 dim, loc, shape, df = self._dist._process_parameters(loc, shape, df)
4738 self.dim, self.loc, self.shape, self.df = dim, loc, shape, df
4739 self.shape_info = _PSD(shape, allow_singular=allow_singular)
4741 def logpdf(self, x):
4742 x = self._dist._process_quantiles(x, self.dim)
4743 U = self.shape_info.U
4744 log_pdet = self.shape_info.log_pdet
4745 return self._dist._logpdf(x, self.loc, U, log_pdet, self.df, self.dim,
4746 self.shape_info.rank)
4748 def cdf(self, x, *, maxpts=None, lower_limit=None, random_state=None):
4749 x = self._dist._process_quantiles(x, self.dim)
4750 return self._dist._cdf(x, self.loc, self.shape, self.df, self.dim,
4751 maxpts, lower_limit, random_state)
4753 def pdf(self, x):
4754 return np.exp(self.logpdf(x))
4756 def rvs(self, size=1, random_state=None):
4757 return self._dist.rvs(loc=self.loc,
4758 shape=self.shape,
4759 df=self.df,
4760 size=size,
4761 random_state=random_state)
4763 def entropy(self):
4764 return self._dist._entropy(self.dim, self.df, self.shape)
4767multivariate_t = multivariate_t_gen()
4770# Set frozen generator docstrings from corresponding docstrings in
4771# multivariate_t_gen and fill in default strings in class docstrings
4772for name in ['logpdf', 'pdf', 'rvs', 'cdf', 'entropy']:
4773 method = multivariate_t_gen.__dict__[name]
4774 method_frozen = multivariate_t_frozen.__dict__[name]
4775 method_frozen.__doc__ = doccer.docformat(method.__doc__,
4776 mvt_docdict_noparams)
4777 method.__doc__ = doccer.docformat(method.__doc__, mvt_docdict_params)
4780_mhg_doc_default_callparams = """\
4781m : array_like
4782 The number of each type of object in the population.
4783 That is, :math:`m[i]` is the number of objects of
4784 type :math:`i`.
4785n : array_like
4786 The number of samples taken from the population.
4787"""
4789_mhg_doc_callparams_note = """\
4790`m` must be an array of positive integers. If the quantile
4791:math:`i` contains values out of the range :math:`[0, m_i]`
4792where :math:`m_i` is the number of objects of type :math:`i`
4793in the population or if the parameters are inconsistent with one
4794another (e.g. ``x.sum() != n``), methods return the appropriate
4795value (e.g. ``0`` for ``pmf``). If `m` or `n` contain negative
4796values, the result will contain ``nan`` there.
4797"""
4799_mhg_doc_frozen_callparams = ""
4801_mhg_doc_frozen_callparams_note = """\
4802See class definition for a detailed description of parameters."""
4804mhg_docdict_params = {
4805 '_doc_default_callparams': _mhg_doc_default_callparams,
4806 '_doc_callparams_note': _mhg_doc_callparams_note,
4807 '_doc_random_state': _doc_random_state
4808}
4810mhg_docdict_noparams = {
4811 '_doc_default_callparams': _mhg_doc_frozen_callparams,
4812 '_doc_callparams_note': _mhg_doc_frozen_callparams_note,
4813 '_doc_random_state': _doc_random_state
4814}
4817class multivariate_hypergeom_gen(multi_rv_generic):
4818 r"""A multivariate hypergeometric random variable.
4820 Methods
4821 -------
4822 pmf(x, m, n)
4823 Probability mass function.
4824 logpmf(x, m, n)
4825 Log of the probability mass function.
4826 rvs(m, n, size=1, random_state=None)
4827 Draw random samples from a multivariate hypergeometric
4828 distribution.
4829 mean(m, n)
4830 Mean of the multivariate hypergeometric distribution.
4831 var(m, n)
4832 Variance of the multivariate hypergeometric distribution.
4833 cov(m, n)
4834 Compute the covariance matrix of the multivariate
4835 hypergeometric distribution.
4837 Parameters
4838 ----------
4839 %(_doc_default_callparams)s
4840 %(_doc_random_state)s
4842 Notes
4843 -----
4844 %(_doc_callparams_note)s
4846 The probability mass function for `multivariate_hypergeom` is
4848 .. math::
4850 P(X_1 = x_1, X_2 = x_2, \ldots, X_k = x_k) = \frac{\binom{m_1}{x_1}
4851 \binom{m_2}{x_2} \cdots \binom{m_k}{x_k}}{\binom{M}{n}}, \\ \quad
4852 (x_1, x_2, \ldots, x_k) \in \mathbb{N}^k \text{ with }
4853 \sum_{i=1}^k x_i = n
4855 where :math:`m_i` are the number of objects of type :math:`i`, :math:`M`
4856 is the total number of objects in the population (sum of all the
4857 :math:`m_i`), and :math:`n` is the size of the sample to be taken
4858 from the population.
4860 .. versionadded:: 1.6.0
4862 Examples
4863 --------
4864 To evaluate the probability mass function of the multivariate
4865 hypergeometric distribution, with a dichotomous population of size
4866 :math:`10` and :math:`20`, at a sample of size :math:`12` with
4867 :math:`8` objects of the first type and :math:`4` objects of the
4868 second type, use:
4870 >>> from scipy.stats import multivariate_hypergeom
4871 >>> multivariate_hypergeom.pmf(x=[8, 4], m=[10, 20], n=12)
4872 0.0025207176631464523
4874 The `multivariate_hypergeom` distribution is identical to the
4875 corresponding `hypergeom` distribution (tiny numerical differences
4876 notwithstanding) when only two types (good and bad) of objects
4877 are present in the population as in the example above. Consider
4878 another example for a comparison with the hypergeometric distribution:
4880 >>> from scipy.stats import hypergeom
4881 >>> multivariate_hypergeom.pmf(x=[3, 1], m=[10, 5], n=4)
4882 0.4395604395604395
4883 >>> hypergeom.pmf(k=3, M=15, n=4, N=10)
4884 0.43956043956044005
4886 The functions ``pmf``, ``logpmf``, ``mean``, ``var``, ``cov``, and ``rvs``
4887 support broadcasting, under the convention that the vector parameters
4888 (``x``, ``m``, and ``n``) are interpreted as if each row along the last
4889 axis is a single object. For instance, we can combine the previous two
4890 calls to `multivariate_hypergeom` as
4892 >>> multivariate_hypergeom.pmf(x=[[8, 4], [3, 1]], m=[[10, 20], [10, 5]],
4893 ... n=[12, 4])
4894 array([0.00252072, 0.43956044])
4896 This broadcasting also works for ``cov``, where the output objects are
4897 square matrices of size ``m.shape[-1]``. For example:
4899 >>> multivariate_hypergeom.cov(m=[[7, 9], [10, 15]], n=[8, 12])
4900 array([[[ 1.05, -1.05],
4901 [-1.05, 1.05]],
4902 [[ 1.56, -1.56],
4903 [-1.56, 1.56]]])
4905 That is, ``result[0]`` is equal to
4906 ``multivariate_hypergeom.cov(m=[7, 9], n=8)`` and ``result[1]`` is equal
4907 to ``multivariate_hypergeom.cov(m=[10, 15], n=12)``.
4909 Alternatively, the object may be called (as a function) to fix the `m`
4910 and `n` parameters, returning a "frozen" multivariate hypergeometric
4911 random variable.
4913 >>> rv = multivariate_hypergeom(m=[10, 20], n=12)
4914 >>> rv.pmf(x=[8, 4])
4915 0.0025207176631464523
4917 See Also
4918 --------
4919 scipy.stats.hypergeom : The hypergeometric distribution.
4920 scipy.stats.multinomial : The multinomial distribution.
4922 References
4923 ----------
4924 .. [1] The Multivariate Hypergeometric Distribution,
4925 http://www.randomservices.org/random/urn/MultiHypergeometric.html
4926 .. [2] Thomas J. Sargent and John Stachurski, 2020,
4927 Multivariate Hypergeometric Distribution
4928 https://python.quantecon.org/_downloads/pdf/multi_hyper.pdf
4929 """
4930 def __init__(self, seed=None):
4931 super().__init__(seed)
4932 self.__doc__ = doccer.docformat(self.__doc__, mhg_docdict_params)
4934 def __call__(self, m, n, seed=None):
4935 """Create a frozen multivariate_hypergeom distribution.
4937 See `multivariate_hypergeom_frozen` for more information.
4938 """
4939 return multivariate_hypergeom_frozen(m, n, seed=seed)
4941 def _process_parameters(self, m, n):
4942 m = np.asarray(m)
4943 n = np.asarray(n)
4944 if m.size == 0:
4945 m = m.astype(int)
4946 if n.size == 0:
4947 n = n.astype(int)
4948 if not np.issubdtype(m.dtype, np.integer):
4949 raise TypeError("'m' must an array of integers.")
4950 if not np.issubdtype(n.dtype, np.integer):
4951 raise TypeError("'n' must an array of integers.")
4952 if m.ndim == 0:
4953 raise ValueError("'m' must be an array with"
4954 " at least one dimension.")
4956 # check for empty arrays
4957 if m.size != 0:
4958 n = n[..., np.newaxis]
4960 m, n = np.broadcast_arrays(m, n)
4962 # check for empty arrays
4963 if m.size != 0:
4964 n = n[..., 0]
4966 mcond = m < 0
4968 M = m.sum(axis=-1)
4970 ncond = (n < 0) | (n > M)
4971 return M, m, n, mcond, ncond, np.any(mcond, axis=-1) | ncond
4973 def _process_quantiles(self, x, M, m, n):
4974 x = np.asarray(x)
4975 if not np.issubdtype(x.dtype, np.integer):
4976 raise TypeError("'x' must an array of integers.")
4977 if x.ndim == 0:
4978 raise ValueError("'x' must be an array with"
4979 " at least one dimension.")
4980 if not x.shape[-1] == m.shape[-1]:
4981 raise ValueError(f"Size of each quantile must be size of 'm': "
4982 f"received {x.shape[-1]}, "
4983 f"but expected {m.shape[-1]}.")
4985 # check for empty arrays
4986 if m.size != 0:
4987 n = n[..., np.newaxis]
4988 M = M[..., np.newaxis]
4990 x, m, n, M = np.broadcast_arrays(x, m, n, M)
4992 # check for empty arrays
4993 if m.size != 0:
4994 n, M = n[..., 0], M[..., 0]
4996 xcond = (x < 0) | (x > m)
4997 return (x, M, m, n, xcond,
4998 np.any(xcond, axis=-1) | (x.sum(axis=-1) != n))
5000 def _checkresult(self, result, cond, bad_value):
5001 result = np.asarray(result)
5002 if cond.ndim != 0:
5003 result[cond] = bad_value
5004 elif cond:
5005 return bad_value
5006 if result.ndim == 0:
5007 return result[()]
5008 return result
5010 def _logpmf(self, x, M, m, n, mxcond, ncond):
5011 # This equation of the pmf comes from the relation,
5012 # n combine r = beta(n+1, 1) / beta(r+1, n-r+1)
5013 num = np.zeros_like(m, dtype=np.float_)
5014 den = np.zeros_like(n, dtype=np.float_)
5015 m, x = m[~mxcond], x[~mxcond]
5016 M, n = M[~ncond], n[~ncond]
5017 num[~mxcond] = (betaln(m+1, 1) - betaln(x+1, m-x+1))
5018 den[~ncond] = (betaln(M+1, 1) - betaln(n+1, M-n+1))
5019 num[mxcond] = np.nan
5020 den[ncond] = np.nan
5021 num = num.sum(axis=-1)
5022 return num - den
5024 def logpmf(self, x, m, n):
5025 """Log of the multivariate hypergeometric probability mass function.
5027 Parameters
5028 ----------
5029 x : array_like
5030 Quantiles, with the last axis of `x` denoting the components.
5031 %(_doc_default_callparams)s
5033 Returns
5034 -------
5035 logpmf : ndarray or scalar
5036 Log of the probability mass function evaluated at `x`
5038 Notes
5039 -----
5040 %(_doc_callparams_note)s
5041 """
5042 M, m, n, mcond, ncond, mncond = self._process_parameters(m, n)
5043 (x, M, m, n, xcond,
5044 xcond_reduced) = self._process_quantiles(x, M, m, n)
5045 mxcond = mcond | xcond
5046 ncond = ncond | np.zeros(n.shape, dtype=np.bool_)
5048 result = self._logpmf(x, M, m, n, mxcond, ncond)
5050 # replace values for which x was out of the domain; broadcast
5051 # xcond to the right shape
5052 xcond_ = xcond_reduced | np.zeros(mncond.shape, dtype=np.bool_)
5053 result = self._checkresult(result, xcond_, -np.inf)
5055 # replace values bad for n or m; broadcast
5056 # mncond to the right shape
5057 mncond_ = mncond | np.zeros(xcond_reduced.shape, dtype=np.bool_)
5058 return self._checkresult(result, mncond_, np.nan)
5060 def pmf(self, x, m, n):
5061 """Multivariate hypergeometric probability mass function.
5063 Parameters
5064 ----------
5065 x : array_like
5066 Quantiles, with the last axis of `x` denoting the components.
5067 %(_doc_default_callparams)s
5069 Returns
5070 -------
5071 pmf : ndarray or scalar
5072 Probability density function evaluated at `x`
5074 Notes
5075 -----
5076 %(_doc_callparams_note)s
5077 """
5078 out = np.exp(self.logpmf(x, m, n))
5079 return out
5081 def mean(self, m, n):
5082 """Mean of the multivariate hypergeometric distribution.
5084 Parameters
5085 ----------
5086 %(_doc_default_callparams)s
5088 Returns
5089 -------
5090 mean : array_like or scalar
5091 The mean of the distribution
5092 """
5093 M, m, n, _, _, mncond = self._process_parameters(m, n)
5094 # check for empty arrays
5095 if m.size != 0:
5096 M, n = M[..., np.newaxis], n[..., np.newaxis]
5097 cond = (M == 0)
5098 M = np.ma.masked_array(M, mask=cond)
5099 mu = n*(m/M)
5100 if m.size != 0:
5101 mncond = (mncond[..., np.newaxis] |
5102 np.zeros(mu.shape, dtype=np.bool_))
5103 return self._checkresult(mu, mncond, np.nan)
5105 def var(self, m, n):
5106 """Variance of the multivariate hypergeometric distribution.
5108 Parameters
5109 ----------
5110 %(_doc_default_callparams)s
5112 Returns
5113 -------
5114 array_like
5115 The variances of the components of the distribution. This is
5116 the diagonal of the covariance matrix of the distribution
5117 """
5118 M, m, n, _, _, mncond = self._process_parameters(m, n)
5119 # check for empty arrays
5120 if m.size != 0:
5121 M, n = M[..., np.newaxis], n[..., np.newaxis]
5122 cond = (M == 0) & (M-1 == 0)
5123 M = np.ma.masked_array(M, mask=cond)
5124 output = n * m/M * (M-m)/M * (M-n)/(M-1)
5125 if m.size != 0:
5126 mncond = (mncond[..., np.newaxis] |
5127 np.zeros(output.shape, dtype=np.bool_))
5128 return self._checkresult(output, mncond, np.nan)
5130 def cov(self, m, n):
5131 """Covariance matrix of the multivariate hypergeometric distribution.
5133 Parameters
5134 ----------
5135 %(_doc_default_callparams)s
5137 Returns
5138 -------
5139 cov : array_like
5140 The covariance matrix of the distribution
5141 """
5142 # see [1]_ for the formula and [2]_ for implementation
5143 # cov( x_i,x_j ) = -n * (M-n)/(M-1) * (K_i*K_j) / (M**2)
5144 M, m, n, _, _, mncond = self._process_parameters(m, n)
5145 # check for empty arrays
5146 if m.size != 0:
5147 M = M[..., np.newaxis, np.newaxis]
5148 n = n[..., np.newaxis, np.newaxis]
5149 cond = (M == 0) & (M-1 == 0)
5150 M = np.ma.masked_array(M, mask=cond)
5151 output = (-n * (M-n)/(M-1) *
5152 np.einsum("...i,...j->...ij", m, m) / (M**2))
5153 # check for empty arrays
5154 if m.size != 0:
5155 M, n = M[..., 0, 0], n[..., 0, 0]
5156 cond = cond[..., 0, 0]
5157 dim = m.shape[-1]
5158 # diagonal entries need to be computed differently
5159 for i in range(dim):
5160 output[..., i, i] = (n * (M-n) * m[..., i]*(M-m[..., i]))
5161 output[..., i, i] = output[..., i, i] / (M-1)
5162 output[..., i, i] = output[..., i, i] / (M**2)
5163 if m.size != 0:
5164 mncond = (mncond[..., np.newaxis, np.newaxis] |
5165 np.zeros(output.shape, dtype=np.bool_))
5166 return self._checkresult(output, mncond, np.nan)
5168 def rvs(self, m, n, size=None, random_state=None):
5169 """Draw random samples from a multivariate hypergeometric distribution.
5171 Parameters
5172 ----------
5173 %(_doc_default_callparams)s
5174 size : integer or iterable of integers, optional
5175 Number of samples to draw. Default is ``None``, in which case a
5176 single variate is returned as an array with shape ``m.shape``.
5177 %(_doc_random_state)s
5179 Returns
5180 -------
5181 rvs : array_like
5182 Random variates of shape ``size`` or ``m.shape``
5183 (if ``size=None``).
5185 Notes
5186 -----
5187 %(_doc_callparams_note)s
5189 Also note that NumPy's `multivariate_hypergeometric` sampler is not
5190 used as it doesn't support broadcasting.
5191 """
5192 M, m, n, _, _, _ = self._process_parameters(m, n)
5194 random_state = self._get_random_state(random_state)
5196 if size is not None and isinstance(size, int):
5197 size = (size, )
5199 if size is None:
5200 rvs = np.empty(m.shape, dtype=m.dtype)
5201 else:
5202 rvs = np.empty(size + (m.shape[-1], ), dtype=m.dtype)
5203 rem = M
5205 # This sampler has been taken from numpy gh-13794
5206 # https://github.com/numpy/numpy/pull/13794
5207 for c in range(m.shape[-1] - 1):
5208 rem = rem - m[..., c]
5209 n0mask = n == 0
5210 rvs[..., c] = (~n0mask *
5211 random_state.hypergeometric(m[..., c],
5212 rem + n0mask,
5213 n + n0mask,
5214 size=size))
5215 n = n - rvs[..., c]
5216 rvs[..., m.shape[-1] - 1] = n
5218 return rvs
5221multivariate_hypergeom = multivariate_hypergeom_gen()
5224class multivariate_hypergeom_frozen(multi_rv_frozen):
5225 def __init__(self, m, n, seed=None):
5226 self._dist = multivariate_hypergeom_gen(seed)
5227 (self.M, self.m, self.n,
5228 self.mcond, self.ncond,
5229 self.mncond) = self._dist._process_parameters(m, n)
5231 # monkey patch self._dist
5232 def _process_parameters(m, n):
5233 return (self.M, self.m, self.n,
5234 self.mcond, self.ncond,
5235 self.mncond)
5236 self._dist._process_parameters = _process_parameters
5238 def logpmf(self, x):
5239 return self._dist.logpmf(x, self.m, self.n)
5241 def pmf(self, x):
5242 return self._dist.pmf(x, self.m, self.n)
5244 def mean(self):
5245 return self._dist.mean(self.m, self.n)
5247 def var(self):
5248 return self._dist.var(self.m, self.n)
5250 def cov(self):
5251 return self._dist.cov(self.m, self.n)
5253 def rvs(self, size=1, random_state=None):
5254 return self._dist.rvs(self.m, self.n,
5255 size=size,
5256 random_state=random_state)
5259# Set frozen generator docstrings from corresponding docstrings in
5260# multivariate_hypergeom and fill in default strings in class docstrings
5261for name in ['logpmf', 'pmf', 'mean', 'var', 'cov', 'rvs']:
5262 method = multivariate_hypergeom_gen.__dict__[name]
5263 method_frozen = multivariate_hypergeom_frozen.__dict__[name]
5264 method_frozen.__doc__ = doccer.docformat(
5265 method.__doc__, mhg_docdict_noparams)
5266 method.__doc__ = doccer.docformat(method.__doc__,
5267 mhg_docdict_params)
5270class random_table_gen(multi_rv_generic):
5271 r"""Contingency tables from independent samples with fixed marginal sums.
5273 This is the distribution of random tables with given row and column vector
5274 sums. This distribution represents the set of random tables under the null
5275 hypothesis that rows and columns are independent. It is used in hypothesis
5276 tests of independence.
5278 Because of assumed independence, the expected frequency of each table
5279 element can be computed from the row and column sums, so that the
5280 distribution is completely determined by these two vectors.
5282 Methods
5283 -------
5284 logpmf(x)
5285 Log-probability of table `x` to occur in the distribution.
5286 pmf(x)
5287 Probability of table `x` to occur in the distribution.
5288 mean(row, col)
5289 Mean table.
5290 rvs(row, col, size=None, method=None, random_state=None)
5291 Draw random tables with given row and column vector sums.
5293 Parameters
5294 ----------
5295 %(_doc_row_col)s
5296 %(_doc_random_state)s
5298 Notes
5299 -----
5300 %(_doc_row_col_note)s
5302 Random elements from the distribution are generated either with Boyett's
5303 [1]_ or Patefield's algorithm [2]_. Boyett's algorithm has
5304 O(N) time and space complexity, where N is the total sum of entries in the
5305 table. Patefield's algorithm has O(K x log(N)) time complexity, where K is
5306 the number of cells in the table and requires only a small constant work
5307 space. By default, the `rvs` method selects the fastest algorithm based on
5308 the input, but you can specify the algorithm with the keyword `method`.
5309 Allowed values are "boyett" and "patefield".
5311 .. versionadded:: 1.10.0
5313 Examples
5314 --------
5315 >>> from scipy.stats import random_table
5317 >>> row = [1, 5]
5318 >>> col = [2, 3, 1]
5319 >>> random_table.mean(row, col)
5320 array([[0.33333333, 0.5 , 0.16666667],
5321 [1.66666667, 2.5 , 0.83333333]])
5323 Alternatively, the object may be called (as a function) to fix the row
5324 and column vector sums, returning a "frozen" distribution.
5326 >>> dist = random_table(row, col)
5327 >>> dist.rvs(random_state=123)
5328 array([[1., 0., 0.],
5329 [1., 3., 1.]])
5331 References
5332 ----------
5333 .. [1] J. Boyett, AS 144 Appl. Statist. 28 (1979) 329-332
5334 .. [2] W.M. Patefield, AS 159 Appl. Statist. 30 (1981) 91-97
5335 """
5337 def __init__(self, seed=None):
5338 super().__init__(seed)
5340 def __call__(self, row, col, *, seed=None):
5341 """Create a frozen distribution of tables with given marginals.
5343 See `random_table_frozen` for more information.
5344 """
5345 return random_table_frozen(row, col, seed=seed)
5347 def logpmf(self, x, row, col):
5348 """Log-probability of table to occur in the distribution.
5350 Parameters
5351 ----------
5352 %(_doc_x)s
5353 %(_doc_row_col)s
5355 Returns
5356 -------
5357 logpmf : ndarray or scalar
5358 Log of the probability mass function evaluated at `x`.
5360 Notes
5361 -----
5362 %(_doc_row_col_note)s
5364 If row and column marginals of `x` do not match `row` and `col`,
5365 negative infinity is returned.
5367 Examples
5368 --------
5369 >>> from scipy.stats import random_table
5370 >>> import numpy as np
5372 >>> x = [[1, 5, 1], [2, 3, 1]]
5373 >>> row = np.sum(x, axis=1)
5374 >>> col = np.sum(x, axis=0)
5375 >>> random_table.logpmf(x, row, col)
5376 -1.6306401200847027
5378 Alternatively, the object may be called (as a function) to fix the row
5379 and column vector sums, returning a "frozen" distribution.
5381 >>> d = random_table(row, col)
5382 >>> d.logpmf(x)
5383 -1.6306401200847027
5384 """
5385 r, c, n = self._process_parameters(row, col)
5386 x = np.asarray(x)
5388 if x.ndim < 2:
5389 raise ValueError("`x` must be at least two-dimensional")
5391 dtype_is_int = np.issubdtype(x.dtype, np.integer)
5392 with np.errstate(invalid='ignore'):
5393 if not dtype_is_int and not np.all(x.astype(int) == x):
5394 raise ValueError("`x` must contain only integral values")
5396 # x does not contain NaN if we arrive here
5397 if np.any(x < 0):
5398 raise ValueError("`x` must contain only non-negative values")
5400 r2 = np.sum(x, axis=-1)
5401 c2 = np.sum(x, axis=-2)
5403 if r2.shape[-1] != len(r):
5404 raise ValueError("shape of `x` must agree with `row`")
5406 if c2.shape[-1] != len(c):
5407 raise ValueError("shape of `x` must agree with `col`")
5409 res = np.empty(x.shape[:-2])
5411 mask = np.all(r2 == r, axis=-1) & np.all(c2 == c, axis=-1)
5413 def lnfac(x):
5414 return gammaln(x + 1)
5416 res[mask] = (np.sum(lnfac(r), axis=-1) + np.sum(lnfac(c), axis=-1)
5417 - lnfac(n) - np.sum(lnfac(x[mask]), axis=(-1, -2)))
5418 res[~mask] = -np.inf
5420 return res[()]
5422 def pmf(self, x, row, col):
5423 """Probability of table to occur in the distribution.
5425 Parameters
5426 ----------
5427 %(_doc_x)s
5428 %(_doc_row_col)s
5430 Returns
5431 -------
5432 pmf : ndarray or scalar
5433 Probability mass function evaluated at `x`.
5435 Notes
5436 -----
5437 %(_doc_row_col_note)s
5439 If row and column marginals of `x` do not match `row` and `col`,
5440 zero is returned.
5442 Examples
5443 --------
5444 >>> from scipy.stats import random_table
5445 >>> import numpy as np
5447 >>> x = [[1, 5, 1], [2, 3, 1]]
5448 >>> row = np.sum(x, axis=1)
5449 >>> col = np.sum(x, axis=0)
5450 >>> random_table.pmf(x, row, col)
5451 0.19580419580419592
5453 Alternatively, the object may be called (as a function) to fix the row
5454 and column vector sums, returning a "frozen" distribution.
5456 >>> d = random_table(row, col)
5457 >>> d.pmf(x)
5458 0.19580419580419592
5459 """
5460 return np.exp(self.logpmf(x, row, col))
5462 def mean(self, row, col):
5463 """Mean of distribution of conditional tables.
5464 %(_doc_mean_params)s
5466 Returns
5467 -------
5468 mean: ndarray
5469 Mean of the distribution.
5471 Notes
5472 -----
5473 %(_doc_row_col_note)s
5475 Examples
5476 --------
5477 >>> from scipy.stats import random_table
5479 >>> row = [1, 5]
5480 >>> col = [2, 3, 1]
5481 >>> random_table.mean(row, col)
5482 array([[0.33333333, 0.5 , 0.16666667],
5483 [1.66666667, 2.5 , 0.83333333]])
5485 Alternatively, the object may be called (as a function) to fix the row
5486 and column vector sums, returning a "frozen" distribution.
5488 >>> d = random_table(row, col)
5489 >>> d.mean()
5490 array([[0.33333333, 0.5 , 0.16666667],
5491 [1.66666667, 2.5 , 0.83333333]])
5492 """
5493 r, c, n = self._process_parameters(row, col)
5494 return np.outer(r, c) / n
5496 def rvs(self, row, col, *, size=None, method=None, random_state=None):
5497 """Draw random tables with fixed column and row marginals.
5499 Parameters
5500 ----------
5501 %(_doc_row_col)s
5502 size : integer, optional
5503 Number of samples to draw (default 1).
5504 method : str, optional
5505 Which method to use, "boyett" or "patefield". If None (default),
5506 selects the fastest method for this input.
5507 %(_doc_random_state)s
5509 Returns
5510 -------
5511 rvs : ndarray
5512 Random 2D tables of shape (`size`, `len(row)`, `len(col)`).
5514 Notes
5515 -----
5516 %(_doc_row_col_note)s
5518 Examples
5519 --------
5520 >>> from scipy.stats import random_table
5522 >>> row = [1, 5]
5523 >>> col = [2, 3, 1]
5524 >>> random_table.rvs(row, col, random_state=123)
5525 array([[1., 0., 0.],
5526 [1., 3., 1.]])
5528 Alternatively, the object may be called (as a function) to fix the row
5529 and column vector sums, returning a "frozen" distribution.
5531 >>> d = random_table(row, col)
5532 >>> d.rvs(random_state=123)
5533 array([[1., 0., 0.],
5534 [1., 3., 1.]])
5535 """
5536 r, c, n = self._process_parameters(row, col)
5537 size, shape = self._process_size_shape(size, r, c)
5539 random_state = self._get_random_state(random_state)
5540 meth = self._process_rvs_method(method, r, c, n)
5542 return meth(r, c, n, size, random_state).reshape(shape)
5544 @staticmethod
5545 def _process_parameters(row, col):
5546 """
5547 Check that row and column vectors are one-dimensional, that they do
5548 not contain negative or non-integer entries, and that the sums over
5549 both vectors are equal.
5550 """
5551 r = np.array(row, dtype=np.int64, copy=True)
5552 c = np.array(col, dtype=np.int64, copy=True)
5554 if np.ndim(r) != 1:
5555 raise ValueError("`row` must be one-dimensional")
5556 if np.ndim(c) != 1:
5557 raise ValueError("`col` must be one-dimensional")
5559 if np.any(r < 0):
5560 raise ValueError("each element of `row` must be non-negative")
5561 if np.any(c < 0):
5562 raise ValueError("each element of `col` must be non-negative")
5564 n = np.sum(r)
5565 if n != np.sum(c):
5566 raise ValueError("sums over `row` and `col` must be equal")
5568 if not np.all(r == np.asarray(row)):
5569 raise ValueError("each element of `row` must be an integer")
5570 if not np.all(c == np.asarray(col)):
5571 raise ValueError("each element of `col` must be an integer")
5573 return r, c, n
5575 @staticmethod
5576 def _process_size_shape(size, r, c):
5577 """
5578 Compute the number of samples to be drawn and the shape of the output
5579 """
5580 shape = (len(r), len(c))
5582 if size is None:
5583 return 1, shape
5585 size = np.atleast_1d(size)
5586 if not np.issubdtype(size.dtype, np.integer) or np.any(size < 0):
5587 raise ValueError("`size` must be a non-negative integer or `None`")
5589 return np.prod(size), tuple(size) + shape
5591 @classmethod
5592 def _process_rvs_method(cls, method, r, c, n):
5593 known_methods = {
5594 None: cls._rvs_select(r, c, n),
5595 "boyett": cls._rvs_boyett,
5596 "patefield": cls._rvs_patefield,
5597 }
5598 try:
5599 return known_methods[method]
5600 except KeyError:
5601 raise ValueError(f"'{method}' not recognized, "
5602 f"must be one of {set(known_methods)}")
5604 @classmethod
5605 def _rvs_select(cls, r, c, n):
5606 fac = 1.0 # benchmarks show that this value is about 1
5607 k = len(r) * len(c) # number of cells
5608 # n + 1 guards against failure if n == 0
5609 if n > fac * np.log(n + 1) * k:
5610 return cls._rvs_patefield
5611 return cls._rvs_boyett
5613 @staticmethod
5614 def _rvs_boyett(row, col, ntot, size, random_state):
5615 return _rcont.rvs_rcont1(row, col, ntot, size, random_state)
5617 @staticmethod
5618 def _rvs_patefield(row, col, ntot, size, random_state):
5619 return _rcont.rvs_rcont2(row, col, ntot, size, random_state)
5622random_table = random_table_gen()
5625class random_table_frozen(multi_rv_frozen):
5626 def __init__(self, row, col, *, seed=None):
5627 self._dist = random_table_gen(seed)
5628 self._params = self._dist._process_parameters(row, col)
5630 # monkey patch self._dist
5631 def _process_parameters(r, c):
5632 return self._params
5633 self._dist._process_parameters = _process_parameters
5635 def logpmf(self, x):
5636 return self._dist.logpmf(x, None, None)
5638 def pmf(self, x):
5639 return self._dist.pmf(x, None, None)
5641 def mean(self):
5642 return self._dist.mean(None, None)
5644 def rvs(self, size=None, method=None, random_state=None):
5645 # optimisations are possible here
5646 return self._dist.rvs(None, None, size=size, method=method,
5647 random_state=random_state)
5650_ctab_doc_row_col = """\
5651row : array_like
5652 Sum of table entries in each row.
5653col : array_like
5654 Sum of table entries in each column."""
5656_ctab_doc_x = """\
5657x : array-like
5658 Two-dimensional table of non-negative integers, or a
5659 multi-dimensional array with the last two dimensions
5660 corresponding with the tables."""
5662_ctab_doc_row_col_note = """\
5663The row and column vectors must be one-dimensional, not empty,
5664and each sum up to the same value. They cannot contain negative
5665or noninteger entries."""
5667_ctab_doc_mean_params = f"""
5668Parameters
5669----------
5670{_ctab_doc_row_col}"""
5672_ctab_doc_row_col_note_frozen = """\
5673See class definition for a detailed description of parameters."""
5675_ctab_docdict = {
5676 "_doc_random_state": _doc_random_state,
5677 "_doc_row_col": _ctab_doc_row_col,
5678 "_doc_x": _ctab_doc_x,
5679 "_doc_mean_params": _ctab_doc_mean_params,
5680 "_doc_row_col_note": _ctab_doc_row_col_note,
5681}
5683_ctab_docdict_frozen = _ctab_docdict.copy()
5684_ctab_docdict_frozen.update({
5685 "_doc_row_col": "",
5686 "_doc_mean_params": "",
5687 "_doc_row_col_note": _ctab_doc_row_col_note_frozen,
5688})
5691def _docfill(obj, docdict, template=None):
5692 obj.__doc__ = doccer.docformat(template or obj.__doc__, docdict)
5695# Set frozen generator docstrings from corresponding docstrings in
5696# random_table and fill in default strings in class docstrings
5697_docfill(random_table_gen, _ctab_docdict)
5698for name in ['logpmf', 'pmf', 'mean', 'rvs']:
5699 method = random_table_gen.__dict__[name]
5700 method_frozen = random_table_frozen.__dict__[name]
5701 _docfill(method_frozen, _ctab_docdict_frozen, method.__doc__)
5702 _docfill(method, _ctab_docdict)
5705class uniform_direction_gen(multi_rv_generic):
5706 r"""A vector-valued uniform direction.
5708 Return a random direction (unit vector). The `dim` keyword specifies
5709 the dimensionality of the space.
5711 Methods
5712 -------
5713 rvs(dim=None, size=1, random_state=None)
5714 Draw random directions.
5716 Parameters
5717 ----------
5718 dim : scalar
5719 Dimension of directions.
5720 seed : {None, int, `numpy.random.Generator`,
5721 `numpy.random.RandomState`}, optional
5723 Used for drawing random variates.
5724 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
5725 If `seed` is an int, a new ``RandomState`` instance is used, seeded
5726 with seed.
5727 If `seed` is already a ``RandomState`` or ``Generator`` instance,
5728 then that object is used.
5729 Default is `None`.
5731 Notes
5732 -----
5733 This distribution generates unit vectors uniformly distributed on
5734 the surface of a hypersphere. These can be interpreted as random
5735 directions.
5736 For example, if `dim` is 3, 3D vectors from the surface of :math:`S^2`
5737 will be sampled.
5739 References
5740 ----------
5741 .. [1] Marsaglia, G. (1972). "Choosing a Point from the Surface of a
5742 Sphere". Annals of Mathematical Statistics. 43 (2): 645-646.
5744 Examples
5745 --------
5746 >>> import numpy as np
5747 >>> from scipy.stats import uniform_direction
5748 >>> x = uniform_direction.rvs(3)
5749 >>> np.linalg.norm(x)
5750 1.
5752 This generates one random direction, a vector on the surface of
5753 :math:`S^2`.
5755 Alternatively, the object may be called (as a function) to return a frozen
5756 distribution with fixed `dim` parameter. Here,
5757 we create a `uniform_direction` with ``dim=3`` and draw 5 observations.
5758 The samples are then arranged in an array of shape 5x3.
5760 >>> rng = np.random.default_rng()
5761 >>> uniform_sphere_dist = uniform_direction(3)
5762 >>> unit_vectors = uniform_sphere_dist.rvs(5, random_state=rng)
5763 >>> unit_vectors
5764 array([[ 0.56688642, -0.1332634 , -0.81294566],
5765 [-0.427126 , -0.74779278, 0.50830044],
5766 [ 0.3793989 , 0.92346629, 0.05715323],
5767 [ 0.36428383, -0.92449076, -0.11231259],
5768 [-0.27733285, 0.94410968, -0.17816678]])
5769 """
5771 def __init__(self, seed=None):
5772 super().__init__(seed)
5773 self.__doc__ = doccer.docformat(self.__doc__)
5775 def __call__(self, dim=None, seed=None):
5776 """Create a frozen n-dimensional uniform direction distribution.
5778 See `uniform_direction` for more information.
5779 """
5780 return uniform_direction_frozen(dim, seed=seed)
5782 def _process_parameters(self, dim):
5783 """Dimension N must be specified; it cannot be inferred."""
5784 if dim is None or not np.isscalar(dim) or dim < 1 or dim != int(dim):
5785 raise ValueError("Dimension of vector must be specified, "
5786 "and must be an integer greater than 0.")
5788 return int(dim)
5790 def rvs(self, dim, size=None, random_state=None):
5791 """Draw random samples from S(N-1).
5793 Parameters
5794 ----------
5795 dim : integer
5796 Dimension of space (N).
5797 size : int or tuple of ints, optional
5798 Given a shape of, for example, (m,n,k), m*n*k samples are
5799 generated, and packed in an m-by-n-by-k arrangement.
5800 Because each sample is N-dimensional, the output shape
5801 is (m,n,k,N). If no shape is specified, a single (N-D)
5802 sample is returned.
5803 random_state : {None, int, `numpy.random.Generator`,
5804 `numpy.random.RandomState`}, optional
5806 Pseudorandom number generator state used to generate resamples.
5808 If `random_state` is ``None`` (or `np.random`), the
5809 `numpy.random.RandomState` singleton is used.
5810 If `random_state` is an int, a new ``RandomState`` instance is
5811 used, seeded with `random_state`.
5812 If `random_state` is already a ``Generator`` or ``RandomState``
5813 instance then that instance is used.
5815 Returns
5816 -------
5817 rvs : ndarray
5818 Random direction vectors
5820 """
5821 random_state = self._get_random_state(random_state)
5822 if size is None:
5823 size = np.array([], dtype=int)
5824 size = np.atleast_1d(size)
5826 dim = self._process_parameters(dim)
5828 samples = _sample_uniform_direction(dim, size, random_state)
5829 return samples
5832uniform_direction = uniform_direction_gen()
5835class uniform_direction_frozen(multi_rv_frozen):
5836 def __init__(self, dim=None, seed=None):
5837 """Create a frozen n-dimensional uniform direction distribution.
5839 Parameters
5840 ----------
5841 dim : int
5842 Dimension of matrices
5843 seed : {None, int, `numpy.random.Generator`,
5844 `numpy.random.RandomState`}, optional
5846 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
5847 singleton is used.
5848 If `seed` is an int, a new ``RandomState`` instance is used,
5849 seeded with `seed`.
5850 If `seed` is already a ``Generator`` or ``RandomState`` instance
5851 then that instance is used.
5853 Examples
5854 --------
5855 >>> from scipy.stats import uniform_direction
5856 >>> x = uniform_direction(3)
5857 >>> x.rvs()
5859 """
5860 self._dist = uniform_direction_gen(seed)
5861 self.dim = self._dist._process_parameters(dim)
5863 def rvs(self, size=None, random_state=None):
5864 return self._dist.rvs(self.dim, size, random_state)
5867def _sample_uniform_direction(dim, size, random_state):
5868 """
5869 Private method to generate uniform directions
5870 Reference: Marsaglia, G. (1972). "Choosing a Point from the Surface of a
5871 Sphere". Annals of Mathematical Statistics. 43 (2): 645-646.
5872 """
5873 samples_shape = np.append(size, dim)
5874 samples = random_state.standard_normal(samples_shape)
5875 samples /= np.linalg.norm(samples, axis=-1, keepdims=True)
5876 return samples
5879_dirichlet_mn_doc_default_callparams = """\
5880alpha : array_like
5881 The concentration parameters. The number of entries along the last axis
5882 determines the dimensionality of the distribution. Each entry must be
5883 strictly positive.
5884n : int or array_like
5885 The number of trials. Each element must be a strictly positive integer.
5886"""
5888_dirichlet_mn_doc_frozen_callparams = ""
5890_dirichlet_mn_doc_frozen_callparams_note = """\
5891See class definition for a detailed description of parameters."""
5893dirichlet_mn_docdict_params = {
5894 '_dirichlet_mn_doc_default_callparams': _dirichlet_mn_doc_default_callparams, # noqa
5895 '_doc_random_state': _doc_random_state
5896}
5898dirichlet_mn_docdict_noparams = {
5899 '_dirichlet_mn_doc_default_callparams': _dirichlet_mn_doc_frozen_callparams, # noqa
5900 '_doc_random_state': _doc_random_state
5901}
5904def _dirichlet_multinomial_check_parameters(alpha, n, x=None):
5906 alpha = np.asarray(alpha)
5907 n = np.asarray(n)
5909 if x is not None:
5910 # Ensure that `x` and `alpha` are arrays. If the shapes are
5911 # incompatible, NumPy will raise an appropriate error.
5912 try:
5913 x, alpha = np.broadcast_arrays(x, alpha)
5914 except ValueError as e:
5915 msg = "`x` and `alpha` must be broadcastable."
5916 raise ValueError(msg) from e
5918 x_int = np.floor(x)
5919 if np.any(x < 0) or np.any(x != x_int):
5920 raise ValueError("`x` must contain only non-negative integers.")
5921 x = x_int
5923 if np.any(alpha <= 0):
5924 raise ValueError("`alpha` must contain only positive values.")
5926 n_int = np.floor(n)
5927 if np.any(n <= 0) or np.any(n != n_int):
5928 raise ValueError("`n` must be a positive integer.")
5929 n = n_int
5931 sum_alpha = np.sum(alpha, axis=-1)
5932 sum_alpha, n = np.broadcast_arrays(sum_alpha, n)
5934 return (alpha, sum_alpha, n) if x is None else (alpha, sum_alpha, n, x)
5937class dirichlet_multinomial_gen(multi_rv_generic):
5938 r"""A Dirichlet multinomial random variable.
5940 The Dirichlet multinomial distribution is a compound probability
5941 distribution: it is the multinomial distribution with number of trials
5942 `n` and class probabilities ``p`` randomly sampled from a Dirichlet
5943 distribution with concentration parameters ``alpha``.
5945 Methods
5946 -------
5947 logpmf(x, alpha, n):
5948 Log of the probability mass function.
5949 pmf(x, alpha, n):
5950 Probability mass function.
5951 mean(alpha, n):
5952 Mean of the Dirichlet multinomial distribution.
5953 var(alpha, n):
5954 Variance of the Dirichlet multinomial distribution.
5955 cov(alpha, n):
5956 The covariance of the Dirichlet multinomial distribution.
5958 Parameters
5959 ----------
5960 %(_dirichlet_mn_doc_default_callparams)s
5961 %(_doc_random_state)s
5963 See Also
5964 --------
5965 scipy.stats.dirichlet : The dirichlet distribution.
5966 scipy.stats.multinomial : The multinomial distribution.
5968 References
5969 ----------
5970 .. [1] Dirichlet-multinomial distribution, Wikipedia,
5971 https://www.wikipedia.org/wiki/Dirichlet-multinomial_distribution
5973 Examples
5974 --------
5975 >>> from scipy.stats import dirichlet_multinomial
5977 Get the PMF
5979 >>> n = 6 # number of trials
5980 >>> alpha = [3, 4, 5] # concentration parameters
5981 >>> x = [1, 2, 3] # counts
5982 >>> dirichlet_multinomial.pmf(x, alpha, n)
5983 0.08484162895927604
5985 If the sum of category counts does not equal the number of trials,
5986 the probability mass is zero.
5988 >>> dirichlet_multinomial.pmf(x, alpha, n=7)
5989 0.0
5991 Get the log of the PMF
5993 >>> dirichlet_multinomial.logpmf(x, alpha, n)
5994 -2.4669689491013327
5996 Get the mean
5998 >>> dirichlet_multinomial.mean(alpha, n)
5999 array([1.5, 2. , 2.5])
6001 Get the variance
6003 >>> dirichlet_multinomial.var(alpha, n)
6004 array([1.55769231, 1.84615385, 2.01923077])
6006 Get the covariance
6008 >>> dirichlet_multinomial.cov(alpha, n)
6009 array([[ 1.55769231, -0.69230769, -0.86538462],
6010 [-0.69230769, 1.84615385, -1.15384615],
6011 [-0.86538462, -1.15384615, 2.01923077]])
6013 Alternatively, the object may be called (as a function) to fix the
6014 `alpha` and `n` parameters, returning a "frozen" Dirichlet multinomial
6015 random variable.
6017 >>> dm = dirichlet_multinomial(alpha, n)
6018 >>> dm.pmf(x)
6019 0.08484162895927579
6021 All methods are fully vectorized. Each element of `x` and `alpha` is
6022 a vector (along the last axis), each element of `n` is an
6023 integer (scalar), and the result is computed element-wise.
6025 >>> x = [[1, 2, 3], [4, 5, 6]]
6026 >>> alpha = [[1, 2, 3], [4, 5, 6]]
6027 >>> n = [6, 15]
6028 >>> dirichlet_multinomial.pmf(x, alpha, n)
6029 array([0.06493506, 0.02626937])
6031 >>> dirichlet_multinomial.cov(alpha, n).shape # both covariance matrices
6032 (2, 3, 3)
6034 Broadcasting according to standard NumPy conventions is supported. Here,
6035 we have four sets of concentration parameters (each a two element vector)
6036 for each of three numbers of trials (each a scalar).
6038 >>> alpha = [[3, 4], [4, 5], [5, 6], [6, 7]]
6039 >>> n = [[6], [7], [8]]
6040 >>> dirichlet_multinomial.mean(alpha, n).shape
6041 (3, 4, 2)
6043 """
6044 def __init__(self, seed=None):
6045 super().__init__(seed)
6046 self.__doc__ = doccer.docformat(self.__doc__,
6047 dirichlet_mn_docdict_params)
6049 def __call__(self, alpha, n, seed=None):
6050 return dirichlet_multinomial_frozen(alpha, n, seed=seed)
6052 def logpmf(self, x, alpha, n):
6053 """The log of the probability mass function.
6055 Parameters
6056 ----------
6057 x: ndarray
6058 Category counts (non-negative integers). Must be broadcastable
6059 with shape parameter ``alpha``. If multidimensional, the last axis
6060 must correspond with the categories.
6061 %(_dirichlet_mn_doc_default_callparams)s
6063 Returns
6064 -------
6065 out: ndarray or scalar
6066 Log of the probability mass function.
6068 """
6070 a, Sa, n, x = _dirichlet_multinomial_check_parameters(alpha, n, x)
6072 out = np.asarray(loggamma(Sa) + loggamma(n + 1) - loggamma(n + Sa))
6073 out += (loggamma(x + a) - (loggamma(a) + loggamma(x + 1))).sum(axis=-1)
6074 np.place(out, n != x.sum(axis=-1), -np.inf)
6075 return out[()]
6077 def pmf(self, x, alpha, n):
6078 """Probability mass function for a Dirichlet multinomial distribution.
6080 Parameters
6081 ----------
6082 x: ndarray
6083 Category counts (non-negative integers). Must be broadcastable
6084 with shape parameter ``alpha``. If multidimensional, the last axis
6085 must correspond with the categories.
6086 %(_dirichlet_mn_doc_default_callparams)s
6088 Returns
6089 -------
6090 out: ndarray or scalar
6091 Probability mass function.
6093 """
6094 return np.exp(self.logpmf(x, alpha, n))
6096 def mean(self, alpha, n):
6097 """Mean of a Dirichlet multinomial distribution.
6099 Parameters
6100 ----------
6101 %(_dirichlet_mn_doc_default_callparams)s
6103 Returns
6104 -------
6105 out: ndarray
6106 Mean of a Dirichlet multinomial distribution.
6108 """
6109 a, Sa, n = _dirichlet_multinomial_check_parameters(alpha, n)
6110 n, Sa = n[..., np.newaxis], Sa[..., np.newaxis]
6111 return n * a / Sa
6113 def var(self, alpha, n):
6114 """The variance of the Dirichlet multinomial distribution.
6116 Parameters
6117 ----------
6118 %(_dirichlet_mn_doc_default_callparams)s
6120 Returns
6121 -------
6122 out: array_like
6123 The variances of the components of the distribution. This is
6124 the diagonal of the covariance matrix of the distribution.
6126 """
6127 a, Sa, n = _dirichlet_multinomial_check_parameters(alpha, n)
6128 n, Sa = n[..., np.newaxis], Sa[..., np.newaxis]
6129 return n * a / Sa * (1 - a/Sa) * (n + Sa) / (1 + Sa)
6131 def cov(self, alpha, n):
6132 """Covariance matrix of a Dirichlet multinomial distribution.
6134 Parameters
6135 ----------
6136 %(_dirichlet_mn_doc_default_callparams)s
6138 Returns
6139 -------
6140 out : array_like
6141 The covariance matrix of the distribution.
6143 """
6144 a, Sa, n = _dirichlet_multinomial_check_parameters(alpha, n)
6145 var = dirichlet_multinomial.var(a, n)
6147 n, Sa = n[..., np.newaxis, np.newaxis], Sa[..., np.newaxis, np.newaxis]
6148 aiaj = a[..., :, np.newaxis] * a[..., np.newaxis, :]
6149 cov = -n * aiaj / Sa ** 2 * (n + Sa) / (1 + Sa)
6151 ii = np.arange(cov.shape[-1])
6152 cov[..., ii, ii] = var
6153 return cov
6156dirichlet_multinomial = dirichlet_multinomial_gen()
6159class dirichlet_multinomial_frozen(multi_rv_frozen):
6160 def __init__(self, alpha, n, seed=None):
6161 alpha, Sa, n = _dirichlet_multinomial_check_parameters(alpha, n)
6162 self.alpha = alpha
6163 self.n = n
6164 self._dist = dirichlet_multinomial_gen(seed)
6166 def logpmf(self, x):
6167 return self._dist.logpmf(x, self.alpha, self.n)
6169 def pmf(self, x):
6170 return self._dist.pmf(x, self.alpha, self.n)
6172 def mean(self):
6173 return self._dist.mean(self.alpha, self.n)
6175 def var(self):
6176 return self._dist.var(self.alpha, self.n)
6178 def cov(self):
6179 return self._dist.cov(self.alpha, self.n)
6182# Set frozen generator docstrings from corresponding docstrings in
6183# dirichlet_multinomial and fill in default strings in class docstrings.
6184for name in ['logpmf', 'pmf', 'mean', 'var', 'cov']:
6185 method = dirichlet_multinomial_gen.__dict__[name]
6186 method_frozen = dirichlet_multinomial_frozen.__dict__[name]
6187 method_frozen.__doc__ = doccer.docformat(
6188 method.__doc__, dirichlet_mn_docdict_noparams)
6189 method.__doc__ = doccer.docformat(method.__doc__,
6190 dirichlet_mn_docdict_params)
6193class vonmises_fisher_gen(multi_rv_generic):
6194 r"""A von Mises-Fisher variable.
6196 The `mu` keyword specifies the mean direction vector. The `kappa` keyword
6197 specifies the concentration parameter.
6199 Methods
6200 -------
6201 pdf(x, mu=None, kappa=1)
6202 Probability density function.
6203 logpdf(x, mu=None, kappa=1)
6204 Log of the probability density function.
6205 rvs(mu=None, kappa=1, size=1, random_state=None)
6206 Draw random samples from a von Mises-Fisher distribution.
6207 entropy(mu=None, kappa=1)
6208 Compute the differential entropy of the von Mises-Fisher distribution.
6209 fit(data)
6210 Fit a von Mises-Fisher distribution to data.
6212 Parameters
6213 ----------
6214 mu : array_like
6215 Mean direction of the distribution. Must be a one-dimensional unit
6216 vector of norm 1.
6217 kappa : float
6218 Concentration parameter. Must be positive.
6219 seed : {None, int, np.random.RandomState, np.random.Generator}, optional
6220 Used for drawing random variates.
6221 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
6222 If `seed` is an int, a new ``RandomState`` instance is used, seeded
6223 with seed.
6224 If `seed` is already a ``RandomState`` or ``Generator`` instance,
6225 then that object is used.
6226 Default is `None`.
6228 See Also
6229 --------
6230 scipy.stats.vonmises : Von-Mises Fisher distribution in 2D on a circle
6231 uniform_direction : uniform distribution on the surface of a hypersphere
6233 Notes
6234 -----
6235 The von Mises-Fisher distribution is a directional distribution on the
6236 surface of the unit hypersphere. The probability density
6237 function of a unit vector :math:`\mathbf{x}` is
6239 .. math::
6241 f(\mathbf{x}) = \frac{\kappa^{d/2-1}}{(2\pi)^{d/2}I_{d/2-1}(\kappa)}
6242 \exp\left(\kappa \mathbf{\mu}^T\mathbf{x}\right),
6244 where :math:`\mathbf{\mu}` is the mean direction, :math:`\kappa` the
6245 concentration parameter, :math:`d` the dimension and :math:`I` the
6246 modified Bessel function of the first kind. As :math:`\mu` represents
6247 a direction, it must be a unit vector or in other words, a point
6248 on the hypersphere: :math:`\mathbf{\mu}\in S^{d-1}`. :math:`\kappa` is a
6249 concentration parameter, which means that it must be positive
6250 (:math:`\kappa>0`) and that the distribution becomes more narrow with
6251 increasing :math:`\kappa`. In that sense, the reciprocal value
6252 :math:`1/\kappa` resembles the variance parameter of the normal
6253 distribution.
6255 The von Mises-Fisher distribution often serves as an analogue of the
6256 normal distribution on the sphere. Intuitively, for unit vectors, a
6257 useful distance measure is given by the angle :math:`\alpha` between
6258 them. This is exactly what the scalar product
6259 :math:`\mathbf{\mu}^T\mathbf{x}=\cos(\alpha)` in the
6260 von Mises-Fisher probability density function describes: the angle
6261 between the mean direction :math:`\mathbf{\mu}` and the vector
6262 :math:`\mathbf{x}`. The larger the angle between them, the smaller the
6263 probability to observe :math:`\mathbf{x}` for this particular mean
6264 direction :math:`\mathbf{\mu}`.
6266 In dimensions 2 and 3, specialized algorithms are used for fast sampling
6267 [2]_, [3]_. For dimenions of 4 or higher the rejection sampling algorithm
6268 described in [4]_ is utilized. This implementation is partially based on
6269 the geomstats package [5]_, [6]_.
6271 .. versionadded:: 1.11
6273 References
6274 ----------
6275 .. [1] Von Mises-Fisher distribution, Wikipedia,
6276 https://en.wikipedia.org/wiki/Von_Mises%E2%80%93Fisher_distribution
6277 .. [2] Mardia, K., and Jupp, P. Directional statistics. Wiley, 2000.
6278 .. [3] J. Wenzel. Numerically stable sampling of the von Mises Fisher
6279 distribution on S2.
6280 https://www.mitsuba-renderer.org/~wenzel/files/vmf.pdf
6281 .. [4] Wood, A. Simulation of the von mises fisher distribution.
6282 Communications in statistics-simulation and computation 23,
6283 1 (1994), 157-164. https://doi.org/10.1080/03610919408813161
6284 .. [5] geomstats, Github. MIT License. Accessed: 06.01.2023.
6285 https://github.com/geomstats/geomstats
6286 .. [6] Miolane, N. et al. Geomstats: A Python Package for Riemannian
6287 Geometry in Machine Learning. Journal of Machine Learning Research
6288 21 (2020). http://jmlr.org/papers/v21/19-027.html
6290 Examples
6291 --------
6292 **Visualization of the probability density**
6294 Plot the probability density in three dimensions for increasing
6295 concentration parameter. The density is calculated by the ``pdf``
6296 method.
6298 >>> import numpy as np
6299 >>> import matplotlib.pyplot as plt
6300 >>> from scipy.stats import vonmises_fisher
6301 >>> from matplotlib.colors import Normalize
6302 >>> n_grid = 100
6303 >>> u = np.linspace(0, np.pi, n_grid)
6304 >>> v = np.linspace(0, 2 * np.pi, n_grid)
6305 >>> u_grid, v_grid = np.meshgrid(u, v)
6306 >>> vertices = np.stack([np.cos(v_grid) * np.sin(u_grid),
6307 ... np.sin(v_grid) * np.sin(u_grid),
6308 ... np.cos(u_grid)],
6309 ... axis=2)
6310 >>> x = np.outer(np.cos(v), np.sin(u))
6311 >>> y = np.outer(np.sin(v), np.sin(u))
6312 >>> z = np.outer(np.ones_like(u), np.cos(u))
6313 >>> def plot_vmf_density(ax, x, y, z, vertices, mu, kappa):
6314 ... vmf = vonmises_fisher(mu, kappa)
6315 ... pdf_values = vmf.pdf(vertices)
6316 ... pdfnorm = Normalize(vmin=pdf_values.min(), vmax=pdf_values.max())
6317 ... ax.plot_surface(x, y, z, rstride=1, cstride=1,
6318 ... facecolors=plt.cm.viridis(pdfnorm(pdf_values)),
6319 ... linewidth=0)
6320 ... ax.set_aspect('equal')
6321 ... ax.view_init(azim=-130, elev=0)
6322 ... ax.axis('off')
6323 ... ax.set_title(rf"$\kappa={kappa}$")
6324 >>> fig, axes = plt.subplots(nrows=1, ncols=3, figsize=(9, 4),
6325 ... subplot_kw={"projection": "3d"})
6326 >>> left, middle, right = axes
6327 >>> mu = np.array([-np.sqrt(0.5), -np.sqrt(0.5), 0])
6328 >>> plot_vmf_density(left, x, y, z, vertices, mu, 5)
6329 >>> plot_vmf_density(middle, x, y, z, vertices, mu, 20)
6330 >>> plot_vmf_density(right, x, y, z, vertices, mu, 100)
6331 >>> plt.subplots_adjust(top=1, bottom=0.0, left=0.0, right=1.0, wspace=0.)
6332 >>> plt.show()
6334 As we increase the concentration parameter, the points are getting more
6335 clustered together around the mean direction.
6337 **Sampling**
6339 Draw 5 samples from the distribution using the ``rvs`` method resulting
6340 in a 5x3 array.
6342 >>> rng = np.random.default_rng()
6343 >>> mu = np.array([0, 0, 1])
6344 >>> samples = vonmises_fisher(mu, 20).rvs(5, random_state=rng)
6345 >>> samples
6346 array([[ 0.3884594 , -0.32482588, 0.86231516],
6347 [ 0.00611366, -0.09878289, 0.99509023],
6348 [-0.04154772, -0.01637135, 0.99900239],
6349 [-0.14613735, 0.12553507, 0.98126695],
6350 [-0.04429884, -0.23474054, 0.97104814]])
6352 These samples are unit vectors on the sphere :math:`S^2`. To verify,
6353 let us calculate their euclidean norms:
6355 >>> np.linalg.norm(samples, axis=1)
6356 array([1., 1., 1., 1., 1.])
6358 Plot 20 observations drawn from the von Mises-Fisher distribution for
6359 increasing concentration parameter :math:`\kappa`. The red dot highlights
6360 the mean direction :math:`\mu`.
6362 >>> def plot_vmf_samples(ax, x, y, z, mu, kappa):
6363 ... vmf = vonmises_fisher(mu, kappa)
6364 ... samples = vmf.rvs(20)
6365 ... ax.plot_surface(x, y, z, rstride=1, cstride=1, linewidth=0,
6366 ... alpha=0.2)
6367 ... ax.scatter(samples[:, 0], samples[:, 1], samples[:, 2], c='k', s=5)
6368 ... ax.scatter(mu[0], mu[1], mu[2], c='r', s=30)
6369 ... ax.set_aspect('equal')
6370 ... ax.view_init(azim=-130, elev=0)
6371 ... ax.axis('off')
6372 ... ax.set_title(rf"$\kappa={kappa}$")
6373 >>> mu = np.array([-np.sqrt(0.5), -np.sqrt(0.5), 0])
6374 >>> fig, axes = plt.subplots(nrows=1, ncols=3,
6375 ... subplot_kw={"projection": "3d"},
6376 ... figsize=(9, 4))
6377 >>> left, middle, right = axes
6378 >>> plot_vmf_samples(left, x, y, z, mu, 5)
6379 >>> plot_vmf_samples(middle, x, y, z, mu, 20)
6380 >>> plot_vmf_samples(right, x, y, z, mu, 100)
6381 >>> plt.subplots_adjust(top=1, bottom=0.0, left=0.0,
6382 ... right=1.0, wspace=0.)
6383 >>> plt.show()
6385 The plots show that with increasing concentration :math:`\kappa` the
6386 resulting samples are centered more closely around the mean direction.
6388 **Fitting the distribution parameters**
6390 The distribution can be fitted to data using the ``fit`` method returning
6391 the estimated parameters. As a toy example let's fit the distribution to
6392 samples drawn from a known von Mises-Fisher distribution.
6394 >>> mu, kappa = np.array([0, 0, 1]), 20
6395 >>> samples = vonmises_fisher(mu, kappa).rvs(1000, random_state=rng)
6396 >>> mu_fit, kappa_fit = vonmises_fisher.fit(samples)
6397 >>> mu_fit, kappa_fit
6398 (array([0.01126519, 0.01044501, 0.99988199]), 19.306398751730995)
6400 We see that the estimated parameters `mu_fit` and `kappa_fit` are
6401 very close to the ground truth parameters.
6403 """
6404 def __init__(self, seed=None):
6405 super().__init__(seed)
6407 def __call__(self, mu=None, kappa=1, seed=None):
6408 """Create a frozen von Mises-Fisher distribution.
6410 See `vonmises_fisher_frozen` for more information.
6411 """
6412 return vonmises_fisher_frozen(mu, kappa, seed=seed)
6414 def _process_parameters(self, mu, kappa):
6415 """
6416 Infer dimensionality from mu and ensure that mu is a one-dimensional
6417 unit vector and kappa positive.
6418 """
6419 mu = np.asarray(mu)
6420 if mu.ndim > 1:
6421 raise ValueError("'mu' must have one-dimensional shape.")
6422 if not np.allclose(np.linalg.norm(mu), 1.):
6423 raise ValueError("'mu' must be a unit vector of norm 1.")
6424 if not mu.size > 1:
6425 raise ValueError("'mu' must have at least two entries.")
6426 kappa_error_msg = "'kappa' must be a positive scalar."
6427 if not np.isscalar(kappa) or kappa < 0:
6428 raise ValueError(kappa_error_msg)
6429 if float(kappa) == 0.:
6430 raise ValueError("For 'kappa=0' the von Mises-Fisher distribution "
6431 "becomes the uniform distribution on the sphere "
6432 "surface. Consider using "
6433 "'scipy.stats.uniform_direction' instead.")
6434 dim = mu.size
6436 return dim, mu, kappa
6438 def _check_data_vs_dist(self, x, dim):
6439 if x.shape[-1] != dim:
6440 raise ValueError("The dimensionality of the last axis of 'x' must "
6441 "match the dimensionality of the "
6442 "von Mises Fisher distribution.")
6443 if not np.allclose(np.linalg.norm(x, axis=-1), 1.):
6444 msg = "'x' must be unit vectors of norm 1 along last dimension."
6445 raise ValueError(msg)
6447 def _log_norm_factor(self, dim, kappa):
6448 # normalization factor is given by
6449 # c = kappa**(dim/2-1)/((2*pi)**(dim/2)*I[dim/2-1](kappa))
6450 # = kappa**(dim/2-1)*exp(-kappa) /
6451 # ((2*pi)**(dim/2)*I[dim/2-1](kappa)*exp(-kappa)
6452 # = kappa**(dim/2-1)*exp(-kappa) /
6453 # ((2*pi)**(dim/2)*ive[dim/2-1](kappa)
6454 # Then the log is given by
6455 # log c = 1/2*(dim -1)*log(kappa) - kappa - -1/2*dim*ln(2*pi) -
6456 # ive[dim/2-1](kappa)
6457 halfdim = 0.5 * dim
6458 return (0.5 * (dim - 2)*np.log(kappa) - halfdim * _LOG_2PI -
6459 np.log(ive(halfdim - 1, kappa)) - kappa)
6461 def _logpdf(self, x, dim, mu, kappa):
6462 """Log of the von Mises-Fisher probability density function.
6464 As this function does no argument checking, it should not be
6465 called directly; use 'logpdf' instead.
6467 """
6468 x = np.asarray(x)
6469 self._check_data_vs_dist(x, dim)
6470 dotproducts = np.einsum('i,...i->...', mu, x)
6471 return self._log_norm_factor(dim, kappa) + kappa * dotproducts
6473 def logpdf(self, x, mu=None, kappa=1):
6474 """Log of the von Mises-Fisher probability density function.
6476 Parameters
6477 ----------
6478 x : array_like
6479 Points at which to evaluate the log of the probability
6480 density function. The last axis of `x` must correspond
6481 to unit vectors of the same dimensionality as the distribution.
6482 mu : array_like, default: None
6483 Mean direction of the distribution. Must be a one-dimensional unit
6484 vector of norm 1.
6485 kappa : float, default: 1
6486 Concentration parameter. Must be positive.
6488 Returns
6489 -------
6490 logpdf : ndarray or scalar
6491 Log of the probability density function evaluated at `x`.
6493 """
6494 dim, mu, kappa = self._process_parameters(mu, kappa)
6495 return self._logpdf(x, dim, mu, kappa)
6497 def pdf(self, x, mu=None, kappa=1):
6498 """Von Mises-Fisher probability density function.
6500 Parameters
6501 ----------
6502 x : array_like
6503 Points at which to evaluate the probability
6504 density function. The last axis of `x` must correspond
6505 to unit vectors of the same dimensionality as the distribution.
6506 mu : array_like
6507 Mean direction of the distribution. Must be a one-dimensional unit
6508 vector of norm 1.
6509 kappa : float
6510 Concentration parameter. Must be positive.
6512 Returns
6513 -------
6514 pdf : ndarray or scalar
6515 Probability density function evaluated at `x`.
6517 """
6518 dim, mu, kappa = self._process_parameters(mu, kappa)
6519 return np.exp(self._logpdf(x, dim, mu, kappa))
6521 def _rvs_2d(self, mu, kappa, size, random_state):
6522 """
6523 In 2D, the von Mises-Fisher distribution reduces to the
6524 von Mises distribution which can be efficiently sampled by numpy.
6525 This method is much faster than the general rejection
6526 sampling based algorithm.
6528 """
6529 mean_angle = np.arctan2(mu[1], mu[0])
6530 angle_samples = random_state.vonmises(mean_angle, kappa, size=size)
6531 samples = np.stack([np.cos(angle_samples), np.sin(angle_samples)],
6532 axis=-1)
6533 return samples
6535 def _rvs_3d(self, kappa, size, random_state):
6536 """
6537 Generate samples from a von Mises-Fisher distribution
6538 with mu = [1, 0, 0] and kappa. Samples then have to be
6539 rotated towards the desired mean direction mu.
6540 This method is much faster than the general rejection
6541 sampling based algorithm.
6542 Reference: https://www.mitsuba-renderer.org/~wenzel/files/vmf.pdf
6544 """
6545 if size is None:
6546 sample_size = 1
6547 else:
6548 sample_size = size
6550 # compute x coordinate acc. to equation from section 3.1
6551 x = random_state.random(sample_size)
6552 x = 1. + np.log(x + (1. - x) * np.exp(-2 * kappa))/kappa
6554 # (y, z) are random 2D vectors that only have to be
6555 # normalized accordingly. Then (x, y z) follow a VMF distribution
6556 temp = np.sqrt(1. - np.square(x))
6557 uniformcircle = _sample_uniform_direction(2, sample_size, random_state)
6558 samples = np.stack([x, temp * uniformcircle[..., 0],
6559 temp * uniformcircle[..., 1]],
6560 axis=-1)
6561 if size is None:
6562 samples = np.squeeze(samples)
6563 return samples
6565 def _rejection_sampling(self, dim, kappa, size, random_state):
6566 """
6567 Generate samples from a n-dimensional von Mises-Fisher distribution
6568 with mu = [1, 0, ..., 0] and kappa via rejection sampling.
6569 Samples then have to be rotated towards the desired mean direction mu.
6570 Reference: https://doi.org/10.1080/03610919408813161
6571 """
6572 dim_minus_one = dim - 1
6573 # calculate number of requested samples
6574 if size is not None:
6575 if not np.iterable(size):
6576 size = (size, )
6577 n_samples = math.prod(size)
6578 else:
6579 n_samples = 1
6580 # calculate envelope for rejection sampler (eq. 4)
6581 sqrt = np.sqrt(4 * kappa ** 2. + dim_minus_one ** 2)
6582 envelop_param = (-2 * kappa + sqrt) / dim_minus_one
6583 if envelop_param == 0:
6584 # the regular formula suffers from loss of precision for high
6585 # kappa. This can only be detected by checking for 0 here.
6586 # Workaround: expansion for sqrt variable
6587 # https://www.wolframalpha.com/input?i=sqrt%284*x%5E2%2Bd%5E2%29
6588 # e = (-2 * k + sqrt(k**2 + d**2)) / d
6589 # ~ (-2 * k + 2 * k + d**2/(4 * k) - d**4/(64 * k**3)) / d
6590 # = d/(4 * k) - d**3/(64 * k**3)
6591 envelop_param = (dim_minus_one/4 * kappa**-1.
6592 - dim_minus_one**3/64 * kappa**-3.)
6593 # reference step 0
6594 node = (1. - envelop_param) / (1. + envelop_param)
6595 # t = ln(1 - ((1-x)/(1+x))**2)
6596 # = ln(4 * x / (1+x)**2)
6597 # = ln(4) + ln(x) - 2*log1p(x)
6598 correction = (kappa * node + dim_minus_one
6599 * (np.log(4) + np.log(envelop_param)
6600 - 2 * np.log1p(envelop_param)))
6601 n_accepted = 0
6602 x = np.zeros((n_samples, ))
6603 halfdim = 0.5 * dim_minus_one
6604 # main loop
6605 while n_accepted < n_samples:
6606 # generate candidates acc. to reference step 1
6607 sym_beta = random_state.beta(halfdim, halfdim,
6608 size=n_samples - n_accepted)
6609 coord_x = (1 - (1 + envelop_param) * sym_beta) / (
6610 1 - (1 - envelop_param) * sym_beta)
6611 # accept or reject: reference step 2
6612 # reformulation for numerical stability:
6613 # t = ln(1 - (1-x)/(1+x) * y)
6614 # = ln((1 + x - y +x*y)/(1 +x))
6615 accept_tol = random_state.random(n_samples - n_accepted)
6616 criterion = (
6617 kappa * coord_x
6618 + dim_minus_one * (np.log((1 + envelop_param - coord_x
6619 + coord_x * envelop_param) / (1 + envelop_param)))
6620 - correction) > np.log(accept_tol)
6621 accepted_iter = np.sum(criterion)
6622 x[n_accepted:n_accepted + accepted_iter] = coord_x[criterion]
6623 n_accepted += accepted_iter
6624 # concatenate x and remaining coordinates: step 3
6625 coord_rest = _sample_uniform_direction(dim_minus_one, n_accepted,
6626 random_state)
6627 coord_rest = np.einsum(
6628 '...,...i->...i', np.sqrt(1 - x ** 2), coord_rest)
6629 samples = np.concatenate([x[..., None], coord_rest], axis=1)
6630 # reshape output to (size, dim)
6631 if size is not None:
6632 samples = samples.reshape(size + (dim, ))
6633 else:
6634 samples = np.squeeze(samples)
6635 return samples
6637 def _rotate_samples(self, samples, mu, dim):
6638 """A QR decomposition is used to find the rotation that maps the
6639 north pole (1, 0,...,0) to the vector mu. This rotation is then
6640 applied to all samples.
6642 Parameters
6643 ----------
6644 samples: array_like, shape = [..., n]
6645 mu : array-like, shape=[n, ]
6646 Point to parametrise the rotation.
6648 Returns
6649 -------
6650 samples : rotated samples
6652 """
6653 base_point = np.zeros((dim, ))
6654 base_point[0] = 1.
6655 embedded = np.concatenate([mu[None, :], np.zeros((dim - 1, dim))])
6656 rotmatrix, _ = np.linalg.qr(np.transpose(embedded))
6657 if np.allclose(np.matmul(rotmatrix, base_point[:, None])[:, 0], mu):
6658 rotsign = 1
6659 else:
6660 rotsign = -1
6662 # apply rotation
6663 samples = np.einsum('ij,...j->...i', rotmatrix, samples) * rotsign
6664 return samples
6666 def _rvs(self, dim, mu, kappa, size, random_state):
6667 if dim == 2:
6668 samples = self._rvs_2d(mu, kappa, size, random_state)
6669 elif dim == 3:
6670 samples = self._rvs_3d(kappa, size, random_state)
6671 else:
6672 samples = self._rejection_sampling(dim, kappa, size,
6673 random_state)
6675 if dim != 2:
6676 samples = self._rotate_samples(samples, mu, dim)
6677 return samples
6679 def rvs(self, mu=None, kappa=1, size=1, random_state=None):
6680 """Draw random samples from a von Mises-Fisher distribution.
6682 Parameters
6683 ----------
6684 mu : array_like
6685 Mean direction of the distribution. Must be a one-dimensional unit
6686 vector of norm 1.
6687 kappa : float
6688 Concentration parameter. Must be positive.
6689 size : int or tuple of ints, optional
6690 Given a shape of, for example, (m,n,k), m*n*k samples are
6691 generated, and packed in an m-by-n-by-k arrangement.
6692 Because each sample is N-dimensional, the output shape
6693 is (m,n,k,N). If no shape is specified, a single (N-D)
6694 sample is returned.
6695 random_state : {None, int, np.random.RandomState, np.random.Generator},
6696 optional
6697 Used for drawing random variates.
6698 If `seed` is `None`, the `~np.random.RandomState` singleton is used.
6699 If `seed` is an int, a new ``RandomState`` instance is used, seeded
6700 with seed.
6701 If `seed` is already a ``RandomState`` or ``Generator`` instance,
6702 then that object is used.
6703 Default is `None`.
6705 Returns
6706 -------
6707 rvs : ndarray
6708 Random variates of shape (`size`, `N`), where `N` is the
6709 dimension of the distribution.
6711 """
6712 dim, mu, kappa = self._process_parameters(mu, kappa)
6713 random_state = self._get_random_state(random_state)
6714 samples = self._rvs(dim, mu, kappa, size, random_state)
6715 return samples
6717 def _entropy(self, dim, kappa):
6718 halfdim = 0.5 * dim
6719 return (-self._log_norm_factor(dim, kappa) - kappa *
6720 ive(halfdim, kappa) / ive(halfdim - 1, kappa))
6722 def entropy(self, mu=None, kappa=1):
6723 """Compute the differential entropy of the von Mises-Fisher
6724 distribution.
6726 Parameters
6727 ----------
6728 mu : array_like, default: None
6729 Mean direction of the distribution. Must be a one-dimensional unit
6730 vector of norm 1.
6731 kappa : float, default: 1
6732 Concentration parameter. Must be positive.
6734 Returns
6735 -------
6736 h : scalar
6737 Entropy of the von Mises-Fisher distribution.
6739 """
6740 dim, _, kappa = self._process_parameters(mu, kappa)
6741 return self._entropy(dim, kappa)
6743 def fit(self, x):
6744 """Fit the von Mises-Fisher distribution to data.
6746 Parameters
6747 ----------
6748 x : array-like
6749 Data the distribution is fitted to. Must be two dimensional.
6750 The second axis of `x` must be unit vectors of norm 1 and
6751 determine the dimensionality of the fitted
6752 von Mises-Fisher distribution.
6754 Returns
6755 -------
6756 mu : ndarray
6757 Estimated mean direction.
6758 kappa : float
6759 Estimated concentration parameter.
6761 """
6762 # validate input data
6763 x = np.asarray(x)
6764 if x.ndim != 2:
6765 raise ValueError("'x' must be two dimensional.")
6766 if not np.allclose(np.linalg.norm(x, axis=-1), 1.):
6767 msg = "'x' must be unit vectors of norm 1 along last dimension."
6768 raise ValueError(msg)
6769 dim = x.shape[-1]
6771 # mu is simply the directional mean
6772 dirstats = directional_stats(x)
6773 mu = dirstats.mean_direction
6774 r = dirstats.mean_resultant_length
6776 # kappa is the solution to the equation:
6777 # r = I[dim/2](kappa) / I[dim/2 -1](kappa)
6778 # = I[dim/2](kappa) * exp(-kappa) / I[dim/2 -1](kappa) * exp(-kappa)
6779 # = ive(dim/2, kappa) / ive(dim/2 -1, kappa)
6781 halfdim = 0.5 * dim
6783 def solve_for_kappa(kappa):
6784 bessel_vals = ive([halfdim, halfdim - 1], kappa)
6785 return bessel_vals[0]/bessel_vals[1] - r
6787 root_res = root_scalar(solve_for_kappa, method="brentq",
6788 bracket=(1e-8, 1e9))
6789 kappa = root_res.root
6790 return mu, kappa
6793vonmises_fisher = vonmises_fisher_gen()
6796class vonmises_fisher_frozen(multi_rv_frozen):
6797 def __init__(self, mu=None, kappa=1, seed=None):
6798 """Create a frozen von Mises-Fisher distribution.
6800 Parameters
6801 ----------
6802 mu : array_like, default: None
6803 Mean direction of the distribution.
6804 kappa : float, default: 1
6805 Concentration parameter. Must be positive.
6806 seed : {None, int, `numpy.random.Generator`,
6807 `numpy.random.RandomState`}, optional
6808 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
6809 singleton is used.
6810 If `seed` is an int, a new ``RandomState`` instance is used,
6811 seeded with `seed`.
6812 If `seed` is already a ``Generator`` or ``RandomState`` instance
6813 then that instance is used.
6815 """
6816 self._dist = vonmises_fisher_gen(seed)
6817 self.dim, self.mu, self.kappa = (
6818 self._dist._process_parameters(mu, kappa)
6819 )
6821 def logpdf(self, x):
6822 """
6823 Parameters
6824 ----------
6825 x : array_like
6826 Points at which to evaluate the log of the probability
6827 density function. The last axis of `x` must correspond
6828 to unit vectors of the same dimensionality as the distribution.
6830 Returns
6831 -------
6832 logpdf : ndarray or scalar
6833 Log of probability density function evaluated at `x`.
6835 """
6836 return self._dist._logpdf(x, self.dim, self.mu, self.kappa)
6838 def pdf(self, x):
6839 """
6840 Parameters
6841 ----------
6842 x : array_like
6843 Points at which to evaluate the log of the probability
6844 density function. The last axis of `x` must correspond
6845 to unit vectors of the same dimensionality as the distribution.
6847 Returns
6848 -------
6849 pdf : ndarray or scalar
6850 Probability density function evaluated at `x`.
6852 """
6853 return np.exp(self.logpdf(x))
6855 def rvs(self, size=1, random_state=None):
6856 """Draw random variates from the Von Mises-Fisher distribution.
6858 Parameters
6859 ----------
6860 size : int or tuple of ints, optional
6861 Given a shape of, for example, (m,n,k), m*n*k samples are
6862 generated, and packed in an m-by-n-by-k arrangement.
6863 Because each sample is N-dimensional, the output shape
6864 is (m,n,k,N). If no shape is specified, a single (N-D)
6865 sample is returned.
6866 random_state : {None, int, `numpy.random.Generator`,
6867 `numpy.random.RandomState`}, optional
6868 If `seed` is None (or `np.random`), the `numpy.random.RandomState`
6869 singleton is used.
6870 If `seed` is an int, a new ``RandomState`` instance is used,
6871 seeded with `seed`.
6872 If `seed` is already a ``Generator`` or ``RandomState`` instance
6873 then that instance is used.
6875 Returns
6876 -------
6877 rvs : ndarray or scalar
6878 Random variates of size (`size`, `N`), where `N` is the
6879 dimension of the distribution.
6881 """
6882 random_state = self._dist._get_random_state(random_state)
6883 return self._dist._rvs(self.dim, self.mu, self.kappa, size,
6884 random_state)
6886 def entropy(self):
6887 """
6888 Calculate the differential entropy of the von Mises-Fisher
6889 distribution.
6891 Returns
6892 -------
6893 h: float
6894 Entropy of the Von Mises-Fisher distribution.
6896 """
6897 return self._dist._entropy(self.dim, self.kappa)