&]\c @`sddlmZmZmZddlZddlZddlmZmZddl Z ddl m Z ddl mZmZmZmZmZddlmZddlmZddlmZdd lmZd d lmZd d lmZd dddddddddg Zej dej!Z"ej dZ#ej ej!Z$dZ%dZ&e'e'dZ(ddZ)de*fdYZ+de*fd YZ,d!e*fd"YZ-d#Z.d$Z/d%Z0d&Z1ie.d'6e/d(6e%d)6Z2ie0d'6e1d(6e%d)6Z3d*e,fd+YZ4e4Z5d,e-fd-YZ6xgd.d/d0d1d2gD]PZ7e4j8e7Z9e6j8e7Z:e j;e9j<e3e:_<e j;e9j<e2e9_<qNWd3Z=d4Z>d%Z?d&Z@ie=d56e>d66e%d)6ZAie?d56e@d66e%d)6ZBd7e,fd8YZCeCZDd9e-fd:YZExad.d/d2gD]PZ7eCj8e7Z9eEj8e7Z:e j;e9j<eBe:_<e j;e9j<eAe9_<q5Wd;ZFd%ZGd&ZHieFd<6e%d)6ZIieGd<6e%d)6ZJd=ZKd>ZLd?ZMd@e,fdAYZNeNZOdBe-fdCYZPxjd.d/d2dDdEdFgD]PZ7eNj8e7Z9ePj8e7Z:e j;e9j<eJe:_<e j;e9j<eIe9_<q,WdGZQd%ZRd%ZSd&ZTieQdH6eRdI6e%d)6ZUieSdH6eTdI6e%d)6ZVdJe,fdKYZWeWZXdLe-fdMYZYxmd.d/dDdNdEd2dFgD]PZ7eWj8e7Z9eYj8e7Z:e j;e9j<eVe:_<e j;e9j<eUe9_<qWeZdOZ[dPeWfdQYZ\e\Z]dRe-fdSYZ^xjd.d/dDdNdEd2gD]PZ7e\j8e7Z9eYj8e7Z:e j;e9j<eVe:_<e j;e9j<eUe9_<qWdTZ_dUZ`d%Zad&Zbie_dH6e`dI6e%d)6ZcieadH6ebdI6e%d)6ZddVe,fdWYZeeeZfdXe-fdYYZgxgdZd[dDd\d2gD]PZ7eej8e7Z9egj8e7Z:e j;e9j<ede:_<e j;e9j<ece9_<qWd]e,fd^YZhehZid_e-fd`YZjdae,fdbYZkekZldce,fddYZmemZndee,fdfYZoeoZpdS(gi(tdivisiontprint_functiontabsolute_importN(tasarray_chkfinitetasarray(tdoccer(tgammalntpsit multigammalntxlogytentr(tcheck_random_state(tdrot(t LinAlgError(tget_lapack_funcsi(tbinom(tmvntmultivariate_normalt matrix_normalt dirichlettwishartt invwishartt multinomialtspecial_ortho_groupt ortho_grouptrandom_correlationt unitary_groupisrandom_state : None or int or np.random.RandomState instance, optional If int or RandomState, use it for drawing the random variates. If None (or np.random), the global np.random state is used. Default is None. cC`s,|j}|jdkr(|d}n|S(s` Remove single-dimensional entries from array and convert to scalar, if necessary. i((tsqueezetndim(tout((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_squeeze_output,s  cC`s|dk r|}n|dkrd|jjj}idd6dd6}||tj|j}n|tjt|}|S(s Determine which eigenvalues are "small" given the spectrum. This is for compatibility across various linear algebra functions that should agree about whether or not a Hermitian matrix is numerically singular and what is its numerical matrix rank. This is designed to be compatible with scipy.linalg.pinvh. Parameters ---------- spectrum : 1d ndarray Array of eigenvalues of a Hermitian matrix. cond, rcond : float, optional Cutoff for small eigenvalues. Singular values smaller than rcond * largest_eigenvalue are considered zero. If None or -1, suitable machine precision is used. Returns ------- eps : float Magnitude cutoff for numerical negligibility. ig@@tfg.AtdN(Ni( tNonetdtypetchartlowertnptfinfotepstmaxtabs(tspectrumtcondtrcondtttfactorR'((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_eigvalsh_to_eps8s   gh㈵>cC`sBtjg|D](}t||kr+dnd|^q dtS(s A helper function for computing the pseudoinverse. Parameters ---------- v : iterable of numbers This may be thought of as a vector of eigenvalues or singular values. eps : float Values with magnitude no greater than eps are considered negligible. Returns ------- v_pinv : 1d float ndarray A vector of pseudo-inverted numbers. iiR"(R%tarrayR)tfloat(tvR'tx((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_pinv_1d[st_PSDcB`s5eZdZddeeedZedZRS(sN Compute coordinated functions of a symmetric positive semidefinite matrix. This class addresses two issues. Firstly it allows the pseudoinverse, the logarithm of the pseudo-determinant, and the rank of the matrix to be computed using one call to eigh instead of three. Secondly it allows these functions to be computed in a way that gives mutually compatible results. All of the functions are computed with a common understanding as to which of the eigenvalues are to be considered negligibly small. The functions are designed to coordinate with scipy.linalg.pinvh() but not necessarily with np.linalg.det() or with np.linalg.matrix_rank(). Parameters ---------- M : array_like Symmetric positive semidefinite matrix (2-D). cond, rcond : float, optional Cutoff for small eigenvalues. Singular values smaller than rcond * largest_eigenvalue are considered zero. If None or -1, suitable machine precision is used. lower : bool, optional Whether the pertinent array data is taken from the lower or upper triangle of M. (Default: lower) check_finite : bool, optional Whether to check that the input matrices contain only finite numbers. Disabling may give a performance gain, but may result in problems (crashes, non-termination) if the inputs do contain infinities or NaNs. allow_singular : bool, optional Whether to allow a singular matrix. (Default: True) Notes ----- The arguments are similar to those of scipy.linalg.pinvh(). c C`s tjj|d|d|\}}t|||} tj|| kr[tdn||| k} t| t|kr| rtjjdnt || } tj |tj | } t| |_ | |_ tjtj| |_d|_dS(NR$t check_finites.the input matrix must be positive semidefinitessingular matrix(tscipytlinalgteighR/R%tmint ValueErrortlenR R4tmultiplytsqrttranktUtsumtlogtlog_pdetR!t_pinv( tselftMR+R,R$R6tallow_singulartstuR'R ts_pinvR@((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt__init__s$ cC`s7|jdkr0tj|j|jj|_n|jS(N(RDR!R%tdotR@tT(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytpinvs!N(t__name__t __module__t__doc__R!tTrueRKtpropertyRN(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR5os& tmulti_rv_genericcB`sDeZdZddZedZejdZdZRS(sd Class which encapsulates common functionality between all multivariate distributions. cC`s&tt|jt||_dS(N(tsuperRTRKR t _random_state(REtseed((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKscC`s|jS(sY Get or set the RandomState object for generating random variates. This can be either None or an existing RandomState object. If None (or np.random), use the RandomState singleton used by np.random. If already a RandomState instance, use it. If an int, use a new RandomState instance seeded with seed. (RV(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt random_states cC`st||_dS(N(R RV(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRXscC`s!|dk rt|S|jSdS(N(R!R RV(RERX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_get_random_states  N( RORPRQR!RKRSRXtsetterRY(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRTs   tmulti_rv_frozencB`s/eZdZedZejdZRS(sj Class which encapsulates common functionality between all frozen multivariate distributions. cC`s |jjS(N(t_distRV(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRXscC`st||j_dS(N(R R\RV(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRXs(RORPRQRSRXRZ(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR[ssmean : array_like, optional Mean of the distribution (default zero) cov : array_like, optional Covariance matrix of the distribution (default one) allow_singular : bool, optional Whether to allow a singular covariance matrix. (Default: False) s2Setting the parameter `mean` to `None` is equivalent to having `mean` be the zero-vector. The parameter `cov` can be a scalar, in which case the covariance matrix is the identity times that value, a vector of diagonal entries for the covariance matrix, or a two-dimensional array_like. ts>See class definition for a detailed description of parameters.t_mvn_doc_default_callparamst_mvn_doc_callparams_notet_doc_random_statetmultivariate_normal_gencB`seZdZddZddeddZdZdZdZ ddedZ ddedZ d Z ddedd d d Z ddedd d d Zddddd ZdddZRS(s A multivariate normal random variable. The `mean` keyword specifies the mean. The `cov` keyword specifies the covariance matrix. Methods ------- ``pdf(x, mean=None, cov=1, allow_singular=False)`` Probability density function. ``logpdf(x, mean=None, cov=1, allow_singular=False)`` Log of the probability density function. ``cdf(x, mean=None, cov=1, allow_singular=False, maxpts=1000000*dim, abseps=1e-5, releps=1e-5)`` Cumulative distribution function. ``logcdf(x, mean=None, cov=1, allow_singular=False, maxpts=1000000*dim, abseps=1e-5, releps=1e-5)`` Log of the cumulative distribution function. ``rvs(mean=None, cov=1, size=1, random_state=None)`` Draw random samples from a multivariate normal distribution. ``entropy()`` Compute the differential entropy of the multivariate normal. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_mvn_doc_default_callparams)s %(_doc_random_state)s Alternatively, the object may be called (as a function) to fix the mean and covariance parameters, returning a "frozen" multivariate normal random variable: rv = multivariate_normal(mean=None, cov=1, allow_singular=False) - Frozen object with the same methods but holding the given mean and covariance fixed. Notes ----- %(_mvn_doc_callparams_note)s The covariance matrix `cov` must be a (symmetric) positive semi-definite matrix. The determinant and inverse of `cov` are computed as the pseudo-determinant and pseudo-inverse, respectively, so that `cov` does not need to have full rank. The probability density function for `multivariate_normal` is .. math:: f(x) = \frac{1}{\sqrt{(2 \pi)^k \det \Sigma}} \exp\left( -\frac{1}{2} (x - \mu)^T \Sigma^{-1} (x - \mu) \right), where :math:`\mu` is the mean, :math:`\Sigma` the covariance matrix, and :math:`k` is the dimension of the space where :math:`x` takes values. .. versionadded:: 0.14.0 Examples -------- >>> import matplotlib.pyplot as plt >>> from scipy.stats import multivariate_normal >>> x = np.linspace(0, 5, 10, endpoint=False) >>> y = multivariate_normal.pdf(x, mean=2.5, cov=0.5); y array([ 0.00108914, 0.01033349, 0.05946514, 0.20755375, 0.43939129, 0.56418958, 0.43939129, 0.20755375, 0.05946514, 0.01033349]) >>> fig1 = plt.figure() >>> ax = fig1.add_subplot(111) >>> ax.plot(x, y) The input quantiles can be any shape of array, as long as the last axis labels the components. This allows us for instance to display the frozen pdf for a non-isotropic random variable in 2D as follows: >>> x, y = np.mgrid[-1:1:.01, -1:1:.01] >>> pos = np.dstack((x, y)) >>> rv = multivariate_normal([0.5, -0.2], [[2.0, 0.3], [0.3, 0.5]]) >>> fig2 = plt.figure() >>> ax2 = fig2.add_subplot(111) >>> ax2.contourf(x, y, rv.pdf(pos)) cC`s2tt|j|tj|jt|_dS(N(RURaRKRt docformatRQtmvn_docdict_params(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK^sicC`st||d|d|S(s Create a frozen multivariate normal distribution. See `multivariate_normal_frozen` for more information. RGRW(tmultivariate_normal_frozen(REtmeantcovRGRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt__call__bs cC`s]|d kr|d krj|d kr-d}qtj|dt}|jdkrZd}q|jd}qtj|dt}|j}ntj|stdn|d krtj |}ntj|dt}|d krd}ntj|dt}|dkr'd |_d |_n|jdksI|jd|kr\td|n|jdkr|tj |}n|jdkrtj |}n|jdkr+|j||fkr+|j\}}||krdt |j}n%d }|t |jt |f}t|n%|jdkrPtd |jn|||fS(s Infer dimensionality from mean or covariance matrix, ensure that mean and covariance are full vector resp. matrix. iR"iis.Dimension of random variable must be a scalar.g?s+Array 'mean' must be a vector of length %d.sHArray 'cov' must be square if it is two dimensional, but cov.shape = %s.sTDimension mismatch: array 'cov' is of shape %s, but 'mean' is a vector of length %d.s>Array 'cov' must be at most two-dimensional, but cov.ndim = %dN(i(ii(R!R%RR1RtshapetsizetisscalarR;tzerosteyetdiagtstrR<(REtdimReRftrowstcolstmsg((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_process_parametersmsP            " $ cC`stj|dt}|jdkr4|tj}nS|jdkr|dkrk|ddtjf}q|tjddf}n|S(sm Adjust quantiles array so that last axis labels the components of each data point. R"iiN(R%RR1Rtnewaxis(RER3Ro((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_process_quantiless cC`sH||}tjtjtj||dd}d|t||S(s Parameters ---------- x : ndarray Points at which to evaluate the log of the probability density function mean : ndarray Mean of the distribution prec_U : ndarray A decomposition such that np.dot(prec_U, prec_U.T) is the precision matrix, i.e. inverse of the covariance matrix. log_det_cov : float Logarithm of the determinant of the covariance matrix rank : int Rank of the covariance matrix. Notes ----- As this function does no argument checking, it should not be called directly; use 'logpdf' instead. taxisig(R%RAtsquareRLt_LOG_2PI(RER3Retprec_Ut log_det_covR?tdevtmaha((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_logpdfs *cC`sp|jd||\}}}|j||}t|d|}|j|||j|j|j}t|S(s Log of the multivariate normal probability density function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_mvn_doc_default_callparams)s Returns ------- pdf : ndarray or scalar Log of the probability density function evaluated at `x` Notes ----- %(_mvn_doc_callparams_note)s RGN( RsR!RuR5R}R@RCR?R(RER3ReRfRGRotpsdR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytlogpdfs $cC`sy|jd||\}}}|j||}t|d|}tj|j|||j|j|j }t |S(s Multivariate normal probability density function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_mvn_doc_default_callparams)s Returns ------- pdf : ndarray or scalar Probability density function evaluated at `x` Notes ----- %(_mvn_doc_callparams_note)s RGN( RsR!RuR5R%texpR}R@RCR?R(RER3ReRfRGRoR~R((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytpdfs -c `sVtjjtj fd}tj|d|}t|S(s Parameters ---------- x : ndarray Points at which to evaluate the cumulative distribution function. mean : ndarray Mean of the distribution cov : array_like Covariance matrix of the distribution maxpts: integer The maximum number of points to use for integration abseps: float Absolute error tolerance releps: float Relative error tolerance Notes ----- As this function does no argument checking, it should not be called directly; use 'cdf' instead. .. versionadded:: 1.0.0 c`s#tj|dS(Ni(Rtmvnun(tx_slice(tabsepsRfR$tmaxptsRetreleps(s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt(si(R%tfullRhtinftapply_along_axisR( RER3ReRfRRRtfunc1dR((RRfR$RReRs8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_cdf sgh㈵>c C`s~|jd||\}}}|j||}t|d||sSd|}ntj|j||||||} | S(s Log of the multivariate normal cumulative distribution function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_mvn_doc_default_callparams)s maxpts: integer, optional The maximum number of points to use for integration (default `1000000*dim`) abseps: float, optional Absolute error tolerance (default 1e-5) releps: float, optional Relative error tolerance (default 1e-5) Returns ------- cdf : ndarray or scalar Log of the cumulative distribution function evaluated at `x` Notes ----- %(_mvn_doc_callparams_note)s .. versionadded:: 1.0.0 RGi@BN(RsR!RuR5R%RBR( RER3ReRfRGRRRRoR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytlogcdf-s 'c C`su|jd||\}}}|j||}t|d||sSd|}n|j||||||} | S(s Multivariate normal cumulative distribution function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_mvn_doc_default_callparams)s maxpts: integer, optional The maximum number of points to use for integration (default `1000000*dim`) abseps: float, optional Absolute error tolerance (default 1e-5) releps: float, optional Relative error tolerance (default 1e-5) Returns ------- cdf : ndarray or scalar Cumulative distribution function evaluated at `x` Notes ----- %(_mvn_doc_callparams_note)s .. versionadded:: 1.0.0 RGi@BN(RsR!RuR5R( RER3ReRfRGRRRRoR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytcdfTs cC`sL|jd||\}}}|j|}|j|||}t|S(s Draw random samples from a multivariate normal distribution. Parameters ---------- %(_mvn_doc_default_callparams)s size : integer, optional Number of samples to draw (default 1). %(_doc_random_state)s Returns ------- rvs : ndarray or scalar Random variates of size (`size`, `N`), where `N` is the dimension of the random variable. Notes ----- %(_mvn_doc_callparams_note)s N(RsR!RYRR(REReRfRiRXRoR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytrvs{scC`sP|jd||\}}}tjjdtjtj|\}}d|S(sP Compute the differential entropy of the multivariate normal. Parameters ---------- %(_mvn_doc_default_callparams)s Returns ------- h : scalar Entropy of the multivariate normal distribution Notes ----- %(_mvn_doc_callparams_note)s ig?N(RsR!R%R8tslogdettpite(REReRfRot_tlogdet((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytentropys*N(RORPRQR!RKtFalseRgRsRuR}RRRRRRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRa sS  ?   & &RdcB`sbeZd ded d dddZdZdZdZdZdd dZ dZ RS( igh㈵>cC`st||_|jjd||\|_|_|_t|jd||_|sgd|j}n||_ ||_ ||_ dS(s& Create a frozen multivariate normal distribution. Parameters ---------- mean : array_like, optional Mean of the distribution (default zero) cov : array_like, optional Covariance matrix of the distribution (default one) allow_singular : bool, optional If this flag is True then tolerate a singular covariance matrix (default False). seed : None or int or np.random.RandomState instance, optional This parameter defines the RandomState object to use for drawing random variates. If None (or np.random), the global np.random state is used. If integer, it is used to seed the local RandomState instance Default is None. maxpts: integer, optional The maximum number of points to use for integration of the cumulative distribution function (default `1000000*dim`) abseps: float, optional Absolute error tolerance for the cumulative distribution function (default 1e-5) releps: float, optional Relative error tolerance for the cumulative distribution function (default 1e-5) Examples -------- When called with the default parameters, this will create a 1D random variable with mean 0 and covariance 1: >>> from scipy.stats import multivariate_normal >>> r = multivariate_normal() >>> r.mean array([ 0.]) >>> r.cov array([[1.]]) RGi@BN( RaR\RsR!RoReRfR5tcov_infoRRR(REReRfRGRWRRR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKs+ !  cC`sU|jj||j}|jj||j|jj|jj|jj}t |S(N( R\RuRoR}ReRR@RCR?R(RER3R((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRscC`stj|j|S(N(R%RR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRscC`stj|j|S(N(R%RBR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRscC`sR|jj||j}|jj||j|j|j|j|j}t |S(N( R\RuRoRReRfRRRR(RER3R((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs$ cC`s|jj|j|j||S(N(R\RReRf(RERiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRscC`s,|jj}|jj}d|td|S(s Computes the differential entropy of the multivariate normal. Returns ------- h : scalar Entropy of the multivariate normal distribution g?i(RRCR?Rx(RERCR?((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs  N( RORPR!RRKRRRRRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRds 4    RRRRRsmean : array_like, optional Mean of the distribution (default: `None`) rowcov : array_like, optional Among-row covariance matrix of the distribution (default: `1`) colcov : array_like, optional Among-column covariance matrix of the distribution (default: `1`) sIf `mean` is set to `None` then a matrix of zeros is used for the mean. The dimensions of this matrix are inferred from the shape of `rowcov` and `colcov`, if these are provided, or set to `1` if ambiguous. `rowcov` and `colcov` can be two-dimensional array_likes specifying the covariance matrices directly. Alternatively, a one-dimensional array will be be interpreted as the entries of a diagonal matrix, and a scalar or zero-dimensional array will be interpreted as this value times the identity matrix. t_matnorm_doc_default_callparamst_matnorm_doc_callparams_notetmatrix_normal_gencB`seZdZd dZd ddd dZdZdZdZd dddZ d dddZ d dddd d Z RS( sd A matrix normal random variable. The `mean` keyword specifies the mean. The `rowcov` keyword specifies the among-row covariance matrix. The 'colcov' keyword specifies the among-column covariance matrix. Methods ------- ``pdf(X, mean=None, rowcov=1, colcov=1)`` Probability density function. ``logpdf(X, mean=None, rowcov=1, colcov=1)`` Log of the probability density function. ``rvs(mean=None, rowcov=1, colcov=1, size=1, random_state=None)`` Draw random samples. Parameters ---------- X : array_like Quantiles, with the last two axes of `X` denoting the components. %(_matnorm_doc_default_callparams)s %(_doc_random_state)s Alternatively, the object may be called (as a function) to fix the mean and covariance parameters, returning a "frozen" matrix normal random variable: rv = matrix_normal(mean=None, rowcov=1, colcov=1) - Frozen object with the same methods but holding the given mean and covariance fixed. Notes ----- %(_matnorm_doc_callparams_note)s The covariance matrices specified by `rowcov` and `colcov` must be (symmetric) positive definite. If the samples in `X` are :math:`m \times n`, then `rowcov` must be :math:`m \times m` and `colcov` must be :math:`n \times n`. `mean` must be the same shape as `X`. The probability density function for `matrix_normal` is .. math:: f(X) = (2 \pi)^{-\frac{mn}{2}}|U|^{-\frac{n}{2}} |V|^{-\frac{m}{2}} \exp\left( -\frac{1}{2} \mathrm{Tr}\left[ U^{-1} (X-M) V^{-1} (X-M)^T \right] \right), where :math:`M` is the mean, :math:`U` the among-row covariance matrix, :math:`V` the among-column covariance matrix. The `allow_singular` behaviour of the `multivariate_normal` distribution is not currently supported. Covariance matrices must be full rank. The `matrix_normal` distribution is closely related to the `multivariate_normal` distribution. Specifically, :math:`\mathrm{Vec}(X)` (the vector formed by concatenating the columns of :math:`X`) has a multivariate normal distribution with mean :math:`\mathrm{Vec}(M)` and covariance :math:`V \otimes U` (where :math:`\otimes` is the Kronecker product). Sampling and pdf evaluation are :math:`\mathcal{O}(m^3 + n^3 + m^2 n + m n^2)` for the matrix normal, but :math:`\mathcal{O}(m^3 n^3)` for the equivalent multivariate normal, making this equivalent form algorithmically inefficient. .. versionadded:: 0.17.0 Examples -------- >>> from scipy.stats import matrix_normal >>> M = np.arange(6).reshape(3,2); M array([[0, 1], [2, 3], [4, 5]]) >>> U = np.diag([1,2,3]); U array([[1, 0, 0], [0, 2, 0], [0, 0, 3]]) >>> V = 0.3*np.identity(2); V array([[ 0.3, 0. ], [ 0. , 0.3]]) >>> X = M + 0.1; X array([[ 0.1, 1.1], [ 2.1, 3.1], [ 4.1, 5.1]]) >>> matrix_normal.pdf(X, mean=M, rowcov=U, colcov=V) 0.023410202050005054 >>> # Equivalent multivariate normal >>> from scipy.stats import multivariate_normal >>> vectorised_X = X.T.flatten() >>> equiv_mean = M.T.flatten() >>> equiv_cov = np.kron(V,U) >>> multivariate_normal.pdf(vectorised_X, mean=equiv_mean, cov=equiv_cov) 0.023410202050005054 cC`s2tt|j|tj|jt|_dS(N(RURRKRRbRQtmatnorm_docdict_params(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKsicC`st|||d|S(sx Create a frozen matrix normal distribution. See `matrix_normal_frozen` for more information. RW(tmatrix_normal_frozen(RERetrowcovtcolcovRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRgsc C`s|dk rrtj|dt}|j}t|dkrKtdntj|dkrrtdqrntj|dt}|jdkr|dk r|tj |d}q|tj d}n!|jdkrtj |}n|j}t|dkrtdn|d|dkr@tdn|ddkr_td n|d}tj|dt}|jdkr|dk r|tj |d}q|tj d}n!|jdkrtj |}n|j}t|dkrtd n|d|dkr7td n|ddkrVtd n|d}|dk r|d|krtd n|d|krtdqntj ||f}||f} | |||fS(s Infer dimensionality from mean or covariance matrices. Handle defaults. Ensure compatible dimensions. R"is%Array `mean` must be two dimensional.isArray `mean` has invalid shape.is(`rowcov` must be a scalar or a 2D array.sArray `rowcov` must be square.s!Array `rowcov` has invalid shape.s(`colcov` must be a scalar or a 2D array.sArray `colcov` must be square.s!Array `colcov` has invalid shape.s=Arrays `mean` and `rowcov` must have the same number of rows.s@Arrays `mean` and `colcov` must have the same number of columns.N( R!R%RR1RhR<R;tanyRtidentityRmRk( REReRRt meanshapetrowshapetnumrowstcolshapetnumcolstdims((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRssZ          cC`sftj|dt}|jdkr@|tjddf}n|jd|krbtdn|S(sq Adjust quantiles array so that last two axes labels the components of each data point. R"iNisJThe shape of array `X` is not compatible with the distribution parameters.(R%RR1RRtRhR;(REtXR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRus c C`s|\}} tj||dddd} tj|jtj| |d} tjtjtj| dddd} d|| t| |||| S(s Parameters ---------- dims : tuple Dimensions of the matrix variates X : ndarray Points at which to evaluate the log of the probability density function mean : ndarray Mean of the distribution row_prec_rt : ndarray A decomposition such that np.dot(row_prec_rt, row_prec_rt.T) is the inverse of the among-row covariance matrix log_det_rowcov : float Logarithm of the determinant of the among-row covariance matrix col_prec_rt : ndarray A decomposition such that np.dot(col_prec_rt, col_prec_rt.T) is the inverse of the among-column covariance matrix log_det_colcov : float Logarithm of the determinant of the among-column covariance matrix Notes ----- As this function does no argument checking, it should not be called directly; use 'logpdf' instead. Rvitstartiig(R%trollaxist tensordotRMRLRARwRx( RERRRet row_prec_rttlog_det_rowcovt col_prec_rttlog_det_colcovRRtroll_devt scale_devR|((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR}s  -c C`s|j|||\}}}}|j||}t|dt}t|dt}|j||||j|j|j|j}t|S(s Log of the matrix normal probability density function. Parameters ---------- X : array_like Quantiles, with the last two axes of `X` denoting the components. %(_matnorm_doc_default_callparams)s Returns ------- logpdf : ndarray Log of the probability density function evaluated at `X` Notes ----- %(_matnorm_doc_callparams_note)s RG(RsRuR5RR}R@RCR( RERReRRRtrowpsdtcolpsdR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR"s ! cC`stj|j||||S(s Matrix normal probability density function. Parameters ---------- X : array_like Quantiles, with the last two axes of `X` denoting the components. %(_matnorm_doc_default_callparams)s Returns ------- pdf : ndarray Probability density function evaluated at `X` Notes ----- %(_matnorm_doc_callparams_note)s (R%RR(RERReRR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR?sc C`st|}|j|||\}}}}tjj|dt}tjj|dt}|j|}|jd|d||df} tj |tj | |j d} tj | j dddd|tj ddddf} |dkr| j|j} n| S(s  Draw random samples from a matrix normal distribution. Parameters ---------- %(_matnorm_doc_default_callparams)s size : integer, optional Number of samples to draw (default 1). %(_doc_random_state)s Returns ------- rvs : ndarray or scalar Random variates of size (`size`, `dims`), where `dims` is the dimension of the random matrices. Notes ----- %(_matnorm_doc_callparams_note)s R$RiiiRvRN(tintRsR7R8tcholeskyRRRYtstandard_normalR%RRLRMRRttreshapeRh( REReRRRiRXRtrowcholtcolcholtstd_normtroll_rvsR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRUs  #$> N( RORPRQR!RKRgRsRuR}RRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR;sb  C  %RcB`s>eZdddddZdZdZdddZRS(icC`sst||_|jj|||\|_|_|_|_t|jdt|_ t|jdt|_ dS(sM Create a frozen matrix normal distribution. Parameters ---------- %(_matnorm_doc_default_callparams)s seed : None or int or np.random.RandomState instance, optional If int or RandomState, use it for drawing the random variates. If None (or np.random), the global np.random state is used. Default is None. Examples -------- >>> from scipy.stats import matrix_normal >>> distn = matrix_normal(mean=np.zeros((3,3))) >>> X = distn.rvs(); X array([[-0.02976962, 0.93339138, -0.09663178], [ 0.67405524, 0.28250467, -0.93308929], [-0.31144782, 0.74535536, 1.30412916]]) >>> distn.pdf(X) 2.5160642368346784e-05 >>> distn.logpdf(X) -10.590229595124615 RGN( RR\RsRReRRR5RRR(REReRRRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK}s0cC`sd|jj||j}|jj|j||j|jj|jj|jj|jj}t |S(N( R\RuRR}ReRR@RCRR(RERR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs !cC`stj|j|S(N(R%RR(RER((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRscC`s%|jj|j|j|j||S(N(R\RReRR(RERiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRsN(RORPR!RKRRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR|s  salpha : array_like The concentration parameters. The number of entries determines the dimensionality of the distribution. t!_dirichlet_doc_default_callparamscC`s_tj|}tj|dkr3tdn(|jdkr[td|jfn|S(Nis%All parameters must be greater than 0is?Parameter vector 'a' must be one dimensional, but a.shape = %s.(R%RR:R;RRh(talpha((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_dirichlet_check_parametersscC`stj|}|jdd|jdkrf|jd|jdkrftd|j|jfn|jd|jdkrtjdtj|dg}|jdkrtj||}q|jdkrtj||f}qtdntj |dkr tdntj |dkrDtdn|dk}|dk}|j|jkrtj ||jdd dj |j}ntj ||}tj|rtd ntjtj|dd d kjrtd tj|dn|S(NiisVector 'x' must have either the same number of entries as, or one entry fewer than, parameter vector 'a', but alpha.shape = %s and x.shape = %s.isUThe input must be one dimensional or a two dimensional matrix containing the entries.s8Each entry in 'x' must be greater than or equal to zero.s/Each entry in 'x' must be smaller or equal one.iRvsJEach entry in 'x' must be greater than zero if its alpha is less than one.g?g& .>sOThe input vector 'x' must lie within the normal simplex. but np.sum(x, 0) = %s.(R%RRhR;R0RARtappendtvstackR:R(trepeatRt logical_andR)R(RR3txktxeq0talphalt1tchk((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_dirichlet_check_inputs48"  .+cC`s&tjt|ttj|S(sz Internal helper function to compute the log of the useful quotient .. math:: B(\alpha) = \frac{\prod_{i=1}{K}\Gamma(\alpha_i)} {\Gamma\left(\sum_{i=1}^{K} \alpha_i \right)} Parameters ---------- %(_dirichlet_doc_default_callparams)s Returns ------- B : scalar Helper quotient, internal use only (R%RAR(R((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_lnBst dirichlet_gencB`skeZdZd dZd dZdZdZdZdZ dZ dZ d d d Z RS( sX A Dirichlet random variable. The `alpha` keyword specifies the concentration parameters of the distribution. .. versionadded:: 0.15.0 Methods ------- ``pdf(x, alpha)`` Probability density function. ``logpdf(x, alpha)`` Log of the probability density function. ``rvs(alpha, size=1, random_state=None)`` Draw random samples from a Dirichlet distribution. ``mean(alpha)`` The mean of the Dirichlet distribution ``var(alpha)`` The variance of the Dirichlet distribution ``entropy(alpha)`` Compute the differential entropy of the Dirichlet distribution. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_dirichlet_doc_default_callparams)s %(_doc_random_state)s Alternatively, the object may be called (as a function) to fix concentration parameters, returning a "frozen" Dirichlet random variable: rv = dirichlet(alpha) - Frozen object with the same methods but holding the given concentration parameters fixed. Notes ----- Each :math:`\alpha` entry must be positive. The distribution has only support on the simplex defined by .. math:: \sum_{i=1}^{K} x_i \le 1 The probability density function for `dirichlet` is .. math:: f(x) = \frac{1}{\mathrm{B}(\boldsymbol\alpha)} \prod_{i=1}^K x_i^{\alpha_i - 1} where .. math:: \mathrm{B}(\boldsymbol\alpha) = \frac{\prod_{i=1}^K \Gamma(\alpha_i)} {\Gamma\bigl(\sum_{i=1}^K \alpha_i\bigr)} and :math:`\boldsymbol\alpha=(\alpha_1,\ldots,\alpha_K)`, the concentration parameters and :math:`K` is the dimension of the space where :math:`x` takes values. Note that the dirichlet interface is somewhat inconsistent. The array returned by the rvs function is transposed with respect to the format expected by the pdf and logpdf. Examples -------- >>> from scipy.stats import dirichlet Generate a dirichlet random variable >>> quantiles = np.array([0.2, 0.2, 0.6]) # specify quantiles >>> alpha = np.array([0.4, 5, 15]) # specify concentration parameters >>> dirichlet.pdf(quantiles, alpha) 0.2843831684937255 The same PDF but following a log scale >>> dirichlet.logpdf(quantiles, alpha) -1.2574327653159187 Once we specify the dirichlet distribution we can then calculate quantities of interest >>> dirichlet.mean(alpha) # get the mean of the distribution array([0.01960784, 0.24509804, 0.73529412]) >>> dirichlet.var(alpha) # get variance array([0.00089829, 0.00864603, 0.00909517]) >>> dirichlet.entropy(alpha) # calculate the differential entropy -4.3280162474082715 We can also return random samples from the distribution >>> dirichlet.rvs(alpha, size=1, random_state=1) array([[0.00766178, 0.24670518, 0.74563305]]) >>> dirichlet.rvs(alpha, size=2, random_state=2) array([[0.01639427, 0.1292273 , 0.85437844], [0.00156917, 0.19033695, 0.80809388]]) cC`s2tt|j|tj|jt|_dS(N(RURRKRRbRQtdirichlet_docdict_params(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKscC`st|d|S(NRW(tdirichlet_frozen(RERRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRgscC`s4t|}| tjt|d|jjdS(sc Parameters ---------- x : ndarray Points at which to evaluate the log of the probability density function %(_dirichlet_doc_default_callparams)s Notes ----- As this function does no argument checking, it should not be called directly; use 'logpdf' instead. ii(RR%RAR RM(RER3RtlnB((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR}s cC`s7t|}t||}|j||}t|S(s Log of the Dirichlet probability density function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_dirichlet_doc_default_callparams)s Returns ------- pdf : ndarray or scalar Log of the probability density function evaluated at `x`. (RRR}R(RER3RR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs cC`s@t|}t||}tj|j||}t|S(st The Dirichlet probability density function. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_dirichlet_doc_default_callparams)s Returns ------- pdf : ndarray or scalar The probability density function evaluated at `x`. (RRR%RR}R(RER3RR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs cC`s)t|}|tj|}t|S(s Compute the mean of the dirichlet distribution. Parameters ---------- %(_dirichlet_doc_default_callparams)s Returns ------- mu : ndarray or scalar Mean of the Dirichlet distribution. (RR%RAR(RERR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRes cC`sCt|}tj|}||||||d}t|S(s Compute the variance of the dirichlet distribution. Parameters ---------- %(_dirichlet_doc_default_callparams)s Returns ------- v : ndarray or scalar Variance of the Dirichlet distribution. i(RR%RAR(RERtalpha0R((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytvars cC`s}t|}tj|}t|}|jd}|||tjj|tj|dtjj|}t|S(s  Compute the differential entropy of the dirichlet distribution. Parameters ---------- %(_dirichlet_doc_default_callparams)s Returns ------- h : scalar Entropy of the Dirichlet distribution ii( RR%RARRhR7tspecialRR(RERRRtKR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs   !icC`s.t|}|j|}|j|d|S(s Draw random samples from a Dirichlet distribution. Parameters ---------- %(_dirichlet_doc_default_callparams)s size : int, optional Number of samples to draw (default 1). %(_doc_random_state)s Returns ------- rvs : ndarray or scalar Random variates of size (`size`, `N`), where `N` is the dimension of the random variable. Ri(RRYR(RERRiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs N( RORPRQR!RKRgR}RRReRRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRsg        RcB`sPeZddZdZdZdZdZdZdddZ RS( cC`s"t||_t||_dS(N(RRRR\(RERRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`s|jj||jS(N(R\RR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR$scC`s|jj||jS(N(R\RR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR'scC`s|jj|jS(N(R\ReR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRe*scC`s|jj|jS(N(R\RR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR-scC`s|jj|jS(N(R\RR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR0sicC`s|jj|j||S(N(R\RR(RERiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR3sN( RORPR!RKRRReRRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs      ReRRsdf : int Degrees of freedom, must be greater than or equal to dimension of the scale matrix scale : array_like Symmetric positive definite scale matrix of the distribution t_doc_default_callparamst_doc_callparams_notet wishart_gencB`seZdZddZddddZdZdZdZdZ dZ dZ d Z d Z d Zd Zd ZdZdZdZdddZdZdZdZRS(sx A Wishart random variable. The `df` keyword specifies the degrees of freedom. The `scale` keyword specifies the scale matrix, which must be symmetric and positive definite. In this context, the scale matrix is often interpreted in terms of a multivariate normal precision matrix (the inverse of the covariance matrix). Methods ------- ``pdf(x, df, scale)`` Probability density function. ``logpdf(x, df, scale)`` Log of the probability density function. ``rvs(df, scale, size=1, random_state=None)`` Draw random samples from a Wishart distribution. ``entropy()`` Compute the differential entropy of the Wishart distribution. Parameters ---------- x : array_like Quantiles, with the last axis of `x` denoting the components. %(_doc_default_callparams)s %(_doc_random_state)s Alternatively, the object may be called (as a function) to fix the degrees of freedom and scale parameters, returning a "frozen" Wishart random variable: rv = wishart(df=1, scale=1) - Frozen object with the same methods but holding the given degrees of freedom and scale fixed. See Also -------- invwishart, chi2 Notes ----- %(_doc_callparams_note)s The scale matrix `scale` must be a symmetric positive definite matrix. Singular matrices, including the symmetric positive semi-definite case, are not supported. The Wishart distribution is often denoted .. math:: W_p(\nu, \Sigma) where :math:`\nu` is the degrees of freedom and :math:`\Sigma` is the :math:`p \times p` scale matrix. The probability density function for `wishart` has support over positive definite matrices :math:`S`; if :math:`S \sim W_p(\nu, \Sigma)`, then its PDF is given by: .. math:: f(S) = \frac{|S|^{\frac{\nu - p - 1}{2}}}{2^{ \frac{\nu p}{2} } |\Sigma|^\frac{\nu}{2} \Gamma_p \left ( \frac{\nu}{2} \right )} \exp\left( -tr(\Sigma^{-1} S) / 2 \right) If :math:`S \sim W_p(\nu, \Sigma)` (Wishart) then :math:`S^{-1} \sim W_p^{-1}(\nu, \Sigma^{-1})` (inverse Wishart). If the scale matrix is 1-dimensional and equal to one, then the Wishart distribution :math:`W_1(\nu, 1)` collapses to the :math:`\chi^2(\nu)` distribution. .. versionadded:: 0.16.0 References ---------- .. [1] M.L. Eaton, "Multivariate Statistics: A Vector Space Approach", Wiley, 1983. .. [2] W.B. Smith and R.R. Hocking, "Algorithm AS 53: Wishart Variate Generator", Applied Statistics, vol. 21, pp. 341-345, 1972. Examples -------- >>> import matplotlib.pyplot as plt >>> from scipy.stats import wishart, chi2 >>> x = np.linspace(1e-5, 8, 100) >>> w = wishart.pdf(x, df=3, scale=1); w[:5] array([ 0.00126156, 0.10892176, 0.14793434, 0.17400548, 0.1929669 ]) >>> c = chi2.pdf(x, 3); c[:5] array([ 0.00126156, 0.10892176, 0.14793434, 0.17400548, 0.1929669 ]) >>> plt.plot(x, w) The input quantiles can be any shape of array, as long as the last axis labels the components. cC`s2tt|j|tj|jt|_dS(N(RURRKRRbRQtwishart_docdict_params(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKscC`st|||S(sl Create a frozen Wishart distribution. See `wishart_frozen` for more information. (twishart_frozen(REtdftscaleRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRgscC`sJ|dkrd}ntj|dt}|jdkrR|tjtjf}n|jdkrstj|}nk|jdkr|jd|jdk rtdt |jn%|jdkrtd|jn|jd}|dkr|}n=tj |stdn||kr=td |n|||fS( Ng?R"iiisLArray 'scale' must be square if it is two dimensional, but scale.scale = %s.sBArray 'scale' must be at most two-dimensional, but scale.ndim = %ds$Degrees of freedom must be a scalar.sMDegrees of freedom cannot be less than dimension of scale matrix, but df = %d( R!R%RR1RRtRmRhR;RnRj(RERRRo((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRss,  *     cC`stj|dt}|jdkrV|tj|ddddtjf}n|jdkr|dkr|tjtjddf}qtj|ddddtjf}n|jdkr+|jd|jdkstdt |jn|ddddtjf}nm|jdkrs|jd|jdkstdt |jqn%|jdkrtd |jn|jdd!||fkstd ||f|jdd!fn|S( sl Adjust quantiles array so that last axis labels the components of each data point. R"iNiisGQuantiles must be square if they are two dimensional, but x.shape = %s.iseQuantiles must be square in the first two dimensions if they are three dimensional, but x.shape = %s.snQuantiles must be at most two-dimensional with an additional dimension for multiplecomponents, but x.ndim = %ds=Quantiles have incompatible dimensions: should be %s, got %s.( R%RR1RRlRtRmRhR;Rn(RER3Ro((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRus.2 ".%#cC`s~tj|}|jdkr.|tj}n.|jdkr\tdtt|n|j}t|}||fS(Niis_Size must be an integer or tuple of integers; thus must have dimension <= 1. Got size.ndim = %s(R%RRRtR;Rnttupletprod(RERitnRh((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt _process_sizes  c C`sStj|jd}tj|j}tj|jd} xt|jdD]} |j|dddd| f\} || >> import matplotlib.pyplot as plt >>> from scipy.stats import invwishart, invgamma >>> x = np.linspace(0.01, 1, 100) >>> iw = invwishart.pdf(x, df=6, scale=1) >>> iw[:3] array([ 1.20546865e-15, 5.42497807e-06, 4.45813929e-03]) >>> ig = invgamma.pdf(x, 6/2., scale=1./2) >>> ig[:3] array([ 1.20546865e-15, 5.42497807e-06, 4.45813929e-03]) >>> plt.plot(x, iw) The input quantiles can be any shape of array, as long as the last axis labels the components. cC`s2tt|j|tj|jt|_dS(N(RUR RKRRbRQR(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`st|||S(sw Create a frozen inverse Wishart distribution. See `invwishart_frozen` for more information. (tinvwishart_frozen(RERRRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRg sc C`sKtj|jd}tj|j}|dkrAt|n d|}tj|jd}xt|jdD]} tjj |dddd| fdt \} } dtj tj | j || >> from scipy.stats import multinomial >>> rv = multinomial(8, [0.3, 0.2, 0.5]) >>> rv.pmf([1, 3, 4]) 0.042000000000000072 The multinomial distribution for :math:`k=2` is identical to the corresponding binomial distribution (tiny numerical differences notwithstanding): >>> from scipy.stats import binom >>> multinomial.pmf([3, 4], n=7, p=[0.4, 0.6]) 0.29030399999999973 >>> binom.pmf(3, 7, 0.4) 0.29030400000000012 The functions ``pmf``, ``logpmf``, ``entropy``, and ``cov`` support broadcasting, under the convention that the vector parameters (``x`` and ``p``) are interpreted as if each row along the last axis is a single object. For instance: >>> multinomial.pmf([[3, 4], [3, 5]], n=[7, 8], p=[.3, .7]) array([0.2268945, 0.25412184]) Here, ``x.shape == (2, 2)``, ``n.shape == (2,)``, and ``p.shape == (2,)``, but following the rules mentioned above they behave as if the rows ``[3, 4]`` and ``[3, 5]`` in ``x`` and ``[.3, .7]`` in ``p`` were a single object, and as if we had ``x.shape = (2,)``, ``n.shape = (2,)``, and ``p.shape = ()``. To obtain the individual elements without broadcasting, we would do this: >>> multinomial.pmf([3, 4], n=7, p=[.3, .7]) 0.2268945 >>> multinomial.pmf([3, 5], 8, p=[.3, .7]) 0.25412184 This broadcasting also works for ``cov``, where the output objects are square matrices of size ``p.shape[-1]``. For example: >>> multinomial.cov([4, 5], [[.3, .7], [.4, .6]]) array([[[ 0.84, -0.84], [-0.84, 0.84]], [[ 1.2 , -1.2 ], [-1.2 , 1.2 ]]]) In this example, ``n.shape == (2,)`` and ``p.shape == (2, 2)``, and following the rules above, these broadcast as if ``p.shape == (2,)``. Thus the result should also be of shape ``(2,)``, but since each output is a :math:`2 \times 2` matrix, the result in fact has shape ``(2, 2, 2)``, where ``result[0]`` is equal to ``multinomial.cov(n=4, p=[.3, .7])`` and ``result[1]`` is equal to ``multinomial.cov(n=5, p=[.4, .6])``. See also -------- scipy.stats.binom : The binomial distribution. numpy.random.multinomial : Sampling from the multinomial distribution. cC`s2tt|j|tj|jt|_dS(N(RURRKRRbRQtmultinomial_docdict_params(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`st|||S(ss Create a frozen multinomial distribution. See `multinomial_frozen` for more information. (tmultinomial_frozen(RERtpRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRg scC`stj|dtjdt}d|dddfjdd|d ...jki(RsR%RtteinsumRRhR!tnan(RERRR&tnnRR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRf> s +cC`s|j||\}}}tjdtj|d!}|tjt|dd}|t|d8}|dtjf}t|j|j|jd}|j d|7_ tjt j |||t|dddd|f}|j |||tj S(s  Compute the entropy of the multinomial distribution. The entropy is computed using this expression: .. math:: f(x) = - \log n! - n\sum_{i=1}^k p_i \log p_i + \sum_{i=1}^k \sum_{x=0}^n \binom n x p_i^x(1-p_i)^{n-x} \log x! Parameters ---------- %(_doc_default_callparams)s Returns ------- h : scalar Entropy of the Multinomial distribution Notes ----- %(_doc_callparams_note)s iRvi.(i(RsR%RR(RAR RRtRRhRR+R!R-(RERRR&R3tterm1tnew_axes_neededtterm2((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRV s )cC`s=|j||\}}}|j|}|j|||S(s Draw random samples from a Multinomial distribution. Parameters ---------- %(_doc_default_callparams)s size : integer or iterable of integers, optional Number of samples to draw (default 1). %(_doc_random_state)s Returns ------- rvs : ndarray or scalar Random variates of shape (`size`, `len(p)`) Notes ----- %(_doc_callparams_note)s (RsRYR(RERRRiRXR&((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR~ sN(RORPRQR!RKRgRsRuR!R"R*R+ReRfRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRI sd      !    (RcB`sVeZdZd dZdZdZdZdZdZ dd dZ RS( s Create a frozen Multinomial distribution. Parameters ---------- n : int number of trials p: array_like probability of a trial falling into each category; should sum to 1 seed : None or int or np.random.RandomState instance, optional This parameter defines the RandomState object to use for drawing random variates. If None (or np.random), the global np.random state is used. If integer, it is used to seed the local RandomState instance Default is None. c`sUt|_jj||\___fd}|j_dS(Nc`sjjjfS(N(RRR&(RR(RE(s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs s(RR\RsRRR&(RERRRWRs((REs8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK s'cC`s|jj||j|jS(N(R\R*RR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR* scC`s|jj||j|jS(N(R\R+RR(RER3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR+ scC`s|jj|j|jS(N(R\ReRR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRe scC`s|jj|j|jS(N(R\RfRR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRf scC`s|jj|j|jS(N(R\RRR(RE((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR sicC`s|jj|j|j||S(N(R\RRR(RERiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR sN( RORPRQR!RKR*R+ReRfRR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR s     R*R+Rftspecial_ortho_group_gencB`sAeZdZddZdddZdZdddZRS(s  A matrix-valued SO(N) random variable. Return a random rotation matrix, drawn from the Haar distribution (the only uniform distribution on SO(n)). The `dim` keyword specifies the dimension N. Methods ------- ``rvs(dim=None, size=1, random_state=None)`` Draw random samples from SO(N). Parameters ---------- dim : scalar Dimension of matrices Notes ---------- This class is wrapping the random_rot code from the MDP Toolkit, https://github.com/mdp-toolkit/mdp-toolkit Return a random rotation matrix, drawn from the Haar distribution (the only uniform distribution on SO(n)). The algorithm is described in the paper Stewart, G.W., "The efficient generation of random orthogonal matrices with an application to condition estimators", SIAM Journal on Numerical Analysis, 17(3), pp. 403-409, 1980. For more information see https://en.wikipedia.org/wiki/Orthogonal_matrix#Randomization See also the similar `ortho_group`. Examples -------- >>> from scipy.stats import special_ortho_group >>> x = special_ortho_group.rvs(3) >>> np.dot(x, x.T) array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16], [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16], [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]]) >>> import scipy.linalg >>> scipy.linalg.det(x) 1.0 This generates one random matrix from SO(3). It is orthogonal and has a determinant of 1. cC`s/tt|j|tj|j|_dS(N(RUR2RKRRbRQ(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`st|d|S(sv Create a frozen SO(N) distribution. See `special_ortho_group_frozen` for more information. RW(tspecial_ortho_group_frozen(RERoRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRg scC`sM|dks:tj| s:|dks:|t|krItdn|S(sG Dimension N must be specified; it cannot be inferred. ismDimension of rotation must be specified, and must be a scalar greater than 1.N(R!R%RjRR;(RERo((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs s: ic C`s|j|}t|}|dkrbtjgt|D]!}|j|ddd|^q:S|j|}tj|}tj|f}xt|dD]}|j d||f}|ddkrtj |dnd||<|dc||tj ||j 7RARRLRRM( RERoRiRXRtHtDRR3tHxtmat((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR! s&   2-+5 N(RORPRQR!RKRgRsR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR2 s 4  R3cB`s&eZdddZdddZRS(cC`s(t||_|jj||_dS(s Create a frozen SO(N) distribution. Parameters ---------- dim : scalar Dimension of matrices seed : None or int or np.random.RandomState instance, optional This parameter defines the RandomState object to use for drawing random variates. If None (or np.random), the global np.random state is used. If integer, it is used to seed the local RandomState instance Default is None. Examples -------- >>> from scipy.stats import special_ortho_group >>> g = special_ortho_group(5) >>> x = g.rvs() N(R2R\RsRo(RERoRW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKP sicC`s|jj|j||S(N(R\RRo(RERiRX((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRi sN(RORPR!RKR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR3O stortho_group_gencB`s2eZdZddZdZdddZRS(s A matrix-valued O(N) random variable. Return a random orthogonal matrix, drawn from the O(N) Haar distribution (the only uniform distribution on O(N)). The `dim` keyword specifies the dimension N. Methods ------- ``rvs(dim=None, size=1, random_state=None)`` Draw random samples from O(N). Parameters ---------- dim : scalar Dimension of matrices Notes ---------- This class is closely related to `special_ortho_group`. Some care is taken to avoid numerical error, as per the paper by Mezzadri. References ---------- .. [1] F. Mezzadri, "How to generate random matrices from the classical compact groups", :arXiv:`math-ph/0609050v2`. Examples -------- >>> from scipy.stats import ortho_group >>> x = ortho_group.rvs(3) >>> np.dot(x, x.T) array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16], [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16], [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]]) >>> import scipy.linalg >>> np.fabs(scipy.linalg.det(x)) 1.0 This generates one random matrix from O(3). It is orthogonal and has a determinant of +1 or -1. cC`s/tt|j|tj|j|_dS(N(RUR:RKRRbRQ(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`sM|dks:tj| s:|dks:|t|krItdn|S(sG Dimension N must be specified; it cannot be inferred. isLDimension of rotation must be specified,and must be a scalar greater than 1.N(R!R%RjRR;(RERo((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs s:ic C`s{|j|}t|}|dkrbtjgt|D]!}|j|ddd|^q:S|j|}tj|}xt|D]}|jd||f}|ddkrtj |dnd}|dc|tj ||j 7<| tj||dtj ||||j } tj|} | | |d|dfRARRL( RERoRiRXRR6RR3R7R8R9((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR s    2)':N(RORPRQR!RKRsR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR:m s/  trandom_correlation_gencB`sGeZdZddZdZdZdZddddZRS( s A random correlation matrix. Return a random correlation matrix, given a vector of eigenvalues. The `eigs` keyword specifies the eigenvalues of the correlation matrix, and implies the dimension. Methods ------- ``rvs(eigs=None, random_state=None)`` Draw random correlation matrices, all with eigenvalues eigs. Parameters ---------- eigs : 1d ndarray Eigenvalues of correlation matrix. Notes ---------- Generates a random correlation matrix following a numerically stable algorithm spelled out by Davies & Higham. This algorithm uses a single O(N) similarity transformation to construct a symmetric positive semi-definite matrix, and applies a series of Givens rotations to scale it to have ones on the diagonal. References ---------- .. [1] Davies, Philip I; Higham, Nicholas J; "Numerically stable generation of correlation matrices and their factors", BIT 2000, Vol. 40, No. 4, pp. 640 651 Examples -------- >>> from scipy.stats import random_correlation >>> np.random.seed(514) >>> x = random_correlation.rvs((.5, .8, 1.2, 1.5)) >>> x array([[ 1. , -0.20387311, 0.18366501, -0.04953711], [-0.20387311, 1. , -0.24351129, 0.06703474], [ 0.18366501, -0.24351129, 1. , 0.38530195], [-0.04953711, 0.06703474, 0.38530195, 1. ]]) >>> import scipy.linalg >>> e, v = scipy.linalg.eigh(x) >>> e array([ 0.5, 0.8, 1.2, 1.5]) cC`s/tt|j|tj|j|_dS(N(RUR;RKRRbRQ(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRK scC`stj|dt}|j}|jdksL|jd|ksL|dkr[tdntjtj|||krtdnx*|D]"}|| krtdqqW||fS(NR"iis7Array 'eigs' must be a vector of length greater than 1.s-Sum of eigenvalues must equal dimensionality.s%All eigenvalues must be non-negative.( R%RR1RiRRhR;tfabsRA(REteigsttolRoR3((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRss ."  c C`s|d}|d}|dkr$dStjt|d||d}|tj|||}dtjd||}|dkrd} n ||} || fS(s&Computes a 2x2 Givens matrix to put 1's on the diagonal. The input matrix is a 2x2 symmetric matrix M = [ aii aij ; aij ajj ]. The output matrix g is a 2x2 anti-symmetric matrix of the form [ c s ; -s c ]; the elements c and s are returned. Applying the output matrix to the input matrix (as b=g.T M g) results in a matrix with bii=1, provided tr(M) - det(M) >= 1 and floating point issues do not occur. Otherwise, some other valid rotation is returned. When tr(M)==2, also bjj=1. g?igi(gg?(tmathR>R(tcopysign( REtaiitajjtaijtaiidtajjdtddR-tcRH((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt _givens_to_1!s   $   cC`s|jjo5|jtjko5|jd|jdksDtn|jd}xt|dD]t}|||fdkrqbn|||fdkrxpt|d|D] }|||fdkrPqqWn8x5t|d|D] }|||fdkrPqqW|j|||f|||f|||f\}}|j }t |||| d|d||ddd||dddt d t t |||| d|d|d|d|d|dt d t qbW|S( s Given a psd matrix m, rotate to put one's on the diagonal, turning it into a correlation matrix. This also requires the trace equal the dimensionality. Note: modifies input matrix iiRtoffxtincxtoffytincyt overwrite_xt overwrite_y( tflagst c_contiguousR"R%RRhR;RRHtravelR RR(REtmR RtjRGRHtmv((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyt_to_corrBs.   9   gvIh%<=gHz>cC`s|j|d|\}}|j|}tj|d|}tjtj|tj||j}|j|}t |j dj |krt dn|S(s Draw random correlation matrices Parameters ---------- eigs : 1d ndarray Eigenvalues of correlation matrix tol : float, optional Tolerance for input parameter checks diag_tol : float, optional Tolerance for deviation of the diagonal of the resulting matrix. Default: 1e-7 Raises ------ RuntimeError Floating point error prevented generating a valid correlation matrix. Returns ------- rvs : ndarray or scalar Random size N-dimensional matrices, dimension (size, dim, dim), each having eigenvalues eigs. R>RXis-Failed to generate a valid correlation matrix( RsRYRRR%RLRmRMRUR)RR(t RuntimeError(RER=RXR>tdiag_tolRoRR((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRks*"N( RORPRQR!RKRsRHRUR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyR; s 2   ! )tunitary_group_gencB`s2eZdZddZdZdddZRS(s A matrix-valued U(N) random variable. Return a random unitary matrix. The `dim` keyword specifies the dimension N. Methods ------- ``rvs(dim=None, size=1, random_state=None)`` Draw random samples from U(N). Parameters ---------- dim : scalar Dimension of matrices Notes ---------- This class is similar to `ortho_group`. References ---------- .. [1] F. Mezzadri, "How to generate random matrices from the classical compact groups", arXiv:math-ph/0609050v2. Examples -------- >>> from scipy.stats import unitary_group >>> x = unitary_group.rvs(3) >>> np.dot(x, x.conj().T) array([[ 1.00000000e+00, 1.13231364e-17, -2.86852790e-16], [ 1.13231364e-17, 1.00000000e+00, -1.46845020e-16], [ -2.86852790e-16, -1.46845020e-16, 1.00000000e+00]]) This generates one random matrix from U(3). The dot product confirms that it is unitary up to machine precision. cC`s/tt|j|tj|j|_dS(N(RURXRKRRbRQ(RERW((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRKscC`sM|dks:tj| s:|dks:|t|krItdn|S(sG Dimension N must be specified; it cannot be inferred. isLDimension of rotation must be specified,and must be a scalar greater than 1.N(R!R%RjRR;(RERo((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRss:ic C`s|j|}t|}|dkrbtjgt|D]!}|j|ddd|^q:S|j|}dtjd|j d||fd|j d||f}t j j |\}}|j }||t|9}|S(sr Draw random samples from U(N). Parameters ---------- dim : integer Dimension of space (N). size : integer, optional Number of samples to draw (default 1). Returns ------- rvs : ndarray or scalar Random size N-dimensional matrices, dimension (size, dim, dim) iRiRXiy?(RYRR%R0RRRsR?R>RR7R8tqrRR)( RERoRiRXRtztqtrR ((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRs   2% N(RORPRQR!RKRsR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pyRXs(  (qt __future__RRRR?tnumpyR%RRt scipy.linalgR7t scipy.miscRt scipy.specialRRRR R tscipy._lib._utilR tscipy.linalg.blasR tscipy.linalg.miscR tscipy.linalg.lapackRt_discrete_distnsRR]Rt__all__RBRRxRt_LOG_PIR`RR!R/R4tobjectR5RTR[R^R_t_mvn_doc_frozen_callparamst_mvn_doc_frozen_callparams_noteRctmvn_docdict_noparamsRaRRdtnamet__dict__tmethodt method_frozenRbRQRRt_matnorm_doc_frozen_callparamst#_matnorm_doc_frozen_callparams_noteRtmatnorm_docdict_noparamsRRRRt _dirichlet_doc_frozen_callparamst%_dirichlet_doc_frozen_callparams_noteRtdirichlet_docdict_noparamsRRRRRRt_wishart_doc_default_callparamst_wishart_doc_callparams_notet_wishart_doc_frozen_callparamst#_wishart_doc_frozen_callparams_noteRtwishart_docdict_noparamsRRRRRR R RRt#_multinomial_doc_default_callparamst _multinomial_doc_callparams_notet"_multinomial_doc_frozen_callparamst'_multinomial_doc_frozen_callparams_noteRtmultinomial_docdict_noparamsRRRR2RR3R:RR;RRXR(((s8lib/python2.7/site-packages/scipy/stats/_multivariate.pytsH   (  # G"   \        ? 2       ,       Y :"   ? G    O 0   y h  [