diff --git a/.gitignore b/.gitignore index 4771e8a695..e9fb3d510a 100644 --- a/.gitignore +++ b/.gitignore @@ -9,4 +9,6 @@ docs/INSTALL.tex docs/README.tex docs/UserGuide.pdf docs/database.tex -develop \ No newline at end of file +develop +testresults/* +sphinxdoc/build/* \ No newline at end of file diff --git a/pymc/distributions.py b/pymc/distributions.py index d8dce607e9..e6c59da267 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -34,6 +34,7 @@ from PyMCObjects import Stochastic, Deterministic from CommonDeterministics import Lambda from numpy import pi, inf +import itertools import pdb import utils import warnings @@ -41,7 +42,8 @@ def poiscdf(a, x): x = np.atleast_1d(x) a = np.resize(a, x.shape) - return np.array([flib.gammq(b,y) for b,y in zip(a.ravel(), x.ravel())]).reshape(x.shape) + values = np.array([flib.gammq(b,y) for b, y in zip(a.ravel(), x.ravel())]) + return values.reshape(x.shape) # Import utility functions import inspect, types @@ -53,32 +55,48 @@ class ArgumentError(AttributeError): """Incorrect class argument""" pass -sc_continuous_distributions = ['bernoulli', 'beta', 'cauchy', 'chi2', 'degenerate', -'exponential', 'exponweib', 'gamma', 'half_normal', 'hypergeometric', -'inverse_gamma', 'laplace', 'logistic', 'lognormal', 'normal', 't', 'uniform', -'weibull', 'skew_normal', 'truncated_normal', 'von_mises'] - -sc_discrete_distributions = ['binomial', 'geometric', 'poisson', 'negative_binomial', 'categorical', 'discrete_uniform', 'truncated_poisson'] - -sc_nonnegative_distributions = ['bernoulli', 'beta', 'chi2', 'exponential', 'exponweib', 'gamma', 'half_normal', 'hypergeometric', 'inverse_gamma', 'lognormal', 'weibull'] - -mv_continuous_distributions = ['dirichlet','inverse_wishart','mv_normal','mv_normal_cov','mv_normal_chol','wishart','wishart_cov'] - -mv_discrete_distributions = ['multivariate_hypergeometric','multinomial'] - -mv_nonnegative_distributions = ['dirichlet', 'inverse_wishart', 'wishart', 'wishart_cov', 'multivariate_hypergeometric', 'multinomial'] - - -availabledistributions = sc_continuous_distributions + sc_discrete_distributions + mv_continuous_distributions + mv_discrete_distributions - -# Changes lower case, underscore-separated names into "Class style" capitalized names -# For example, 'negative_binomial' becomes 'NegativeBinomial' +sc_continuous_distributions = ['bernoulli', 'beta', 'cauchy', 'chi2', + 'degenerate', 'exponential', 'exponweib', + 'gamma', 'half_normal', 'hypergeometric', + 'inverse_gamma', 'laplace', 'logistic', + 'lognormal', 'normal', 't', 'uniform', + 'weibull', 'skew_normal', 'truncated_normal', + 'von_mises'] + +sc_discrete_distributions = ['binomial', 'geometric', 'poisson', + 'negative_binomial', 'categorical', + 'discrete_uniform', 'truncated_poisson'] + +sc_nonnegative_distributions = ['bernoulli', 'beta', 'chi2', 'exponential', + 'exponweib', 'gamma', 'half_normal', + 'hypergeometric', 'inverse_gamma', 'lognormal', + 'weibull'] + +mv_continuous_distributions = ['dirichlet', 'inverse_wishart', 'mv_normal', + 'mv_normal_cov', 'mv_normal_chol', 'wishart', + 'wishart_cov', 'inverse_wishart_prec'] + +mv_discrete_distributions = ['multivariate_hypergeometric', 'multinomial'] + +mv_nonnegative_distributions = ['dirichlet', 'inverse_wishart', 'wishart', + 'wishart_cov', 'multivariate_hypergeometric', + 'multinomial'] + +availabledistributions = (sc_continuous_distributions + + sc_discrete_distributions + + mv_continuous_distributions + + mv_discrete_distributions) + +# Changes lower case, underscore-separated names into "Class style" +# capitalized names For example, 'negative_binomial' becomes +# 'NegativeBinomial' capitalize = lambda name: ''.join([s.capitalize() for s in name.split('_')]) -# ============================================================================================ -# = User-accessible function to convert a logp and random function to a Stochastic subclass. = -# ============================================================================================ +# ============================================================================== +# User-accessible function to convert a logp and random function to a +# Stochastic subclass. +# ============================================================================== # TODO Document this function def bind_size(randfun, shape): @@ -123,7 +141,7 @@ def new_dist_class(*new_class_args): (dtype, name, parent_names, parents_default, docstr, logp, random, mv) = new_class_args class new_class(Stochastic): - + def __init__(self, *args, **kwds): (dtype, name, parent_names, parents_default, docstr, logp, random, mv) = new_class_args parents=parents_default @@ -254,9 +272,9 @@ def shape_error(): new_class.__doc__ = docstr new_class.mv = mv new_class.raw_fns = {'logp': logp, 'random': random} - + return new_class - + def stochastic_from_dist(name, logp, random=None, dtype=np.float, mv=False): """ @@ -310,7 +328,8 @@ def stochastic_from_dist(name, logp, random=None, dtype=np.float, mv=False): docstr += logp.__doc__ logp=valuewrapper(logp) - return new_dist_class(dtype, name, parent_names, parents_default, docstr, logp, random, mv) + return new_dist_class(dtype, name, parent_names, parents_default, docstr, + logp, random, mv) #------------------------------------------------------------- @@ -657,7 +676,8 @@ def beta_like(x, alpha, beta): R""" beta_like(x, alpha, beta) - Beta log-likelihood. The conjugate prior for the parameter :math: `p` of the binomial distribution. + Beta log-likelihood. The conjugate prior for the parameter + :math:`p` of the binomial distribution. .. math:: f(x \mid \alpha, \beta) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha) \Gamma(\beta)} x^{\alpha - 1} (1 - x)^{\beta - 1} @@ -753,7 +773,9 @@ def betabin_like(x, alpha, beta, n): R""" betabin_like(x, alpha, beta) - Beta-binomial log-likelihood. Equivalent to binomial random variables with probabilities drawn from a :math: `\texttt{Beta}(\alpha,\beta)` distribution. + Beta-binomial log-likelihood. Equivalent to binomial random + variables with probabilities drawn from a + :math:`\texttt{Beta}(\alpha,\beta)` distribution. .. math:: f(x \mid \alpha, \beta, n) = \frac{\Gamma(\alpha + \beta)}{\Gamma(\alpha)} \frac{\Gamma(n+1)}{\Gamma(x+1)\Gamma(n-x+1)} \frac{\Gamma(\alpha + x)\Gamma(n+\beta-x)}{\Gamma(\alpha+\beta+n)} @@ -964,18 +986,16 @@ def dirichlet_like(x, theta): \cdot\left(1-\sum_{i=1}^{k-1}x_i\right)^\theta_k :Parameters: - x : (n, k-1) array - Array of shape (n, k-1) where `n` is the number of samples - and `k` the dimension. + x : (n, k-1) array + Array of shape (n, k-1) where `n` is the number of samples + and `k` the dimension. :math:`0 < x_i < 1`, :math:`\sum_{i=1}^{k-1} x_i < 1` theta : array An (n,k) or (1,k) array > 0. - - .. note:: - Only the first `k-1` elements of `x` are expected. Can be used as a parent of Multinomial and Categorical - nevertheless. - + .. note:: + Only the first `k-1` elements of `x` are expected. Can be used + as a parent of Multinomial and Categorical nevertheless. """ x = np.atleast_2d(x) theta = np.atleast_2d(theta) @@ -1364,49 +1384,99 @@ def inverse_gamma_like(x, alpha, beta): return flib.igamma(x, alpha, beta) # Inverse Wishart--------------------------------------------------- -def rinverse_wishart(n, Tau): + +def rinverse_wishart(n, C): """ - rinverse_wishart(n, Tau) + rinverse_wishart(n, C) Return an inverse Wishart random matrix. n is the degrees of freedom. - Tau is a positive definite scale matrix. + C is a positive definite scale matrix. """ - wi = rwishart(n, np.asmatrix(Tau).I).I + wi = rwishart(n, np.asmatrix(C).I).I flib.symmetrize(wi) return wi -def inverse_wishart_expval(n, Tau): +def inverse_wishart_expval(n, C): """ - inverse_wishart_expval(n, Tau) + inverse_wishart_expval(n, C) + + :Parameters: + - `n` : [int] Degrees of freedom (n > 0). + - `C` : Symmetric and positive definite scale matrix Expected value of inverse Wishart distribution. """ - return np.asarray(Tau)/(n-len(Tau)-1) + return np.asarray(C)/(n-len(C)-1) -def inverse_wishart_like(X, n, Tau): +def inverse_wishart_like(X, n, C): R""" - inverse_wishart_like(X, n, Tau) + inverse_wishart_like(X, n, C) - Inverse Wishart log-likelihood. The inverse Wishart distribution is the conjugate - prior for the covariance matrix of a multivariate normal distribution. + Inverse Wishart log-likelihood. The inverse Wishart distribution + is the conjugate prior for the covariance matrix of a multivariate + normal distribution. .. math:: - f(X \mid n, T) = \frac{{\mid T \mid}^{n/2}{\mid X \mid}^{(n-k-1)/2} \exp\left\{ -\frac{1}{2} Tr(TX^{-1}) \right\}}{2^{nk/2} \Gamma_p(n/2)} + f(X \mid n, T) = \frac{{\mid T \mid}^{n/2}{\mid X + \mid}^{(n-k-1)/2} \exp\left\{ -\frac{1}{2} Tr(TX^{-1}) + \right\}}{2^{nk/2} \Gamma_p(n/2)} where :math:`k` is the rank of X. :Parameters: - `X` : Symmetric, positive definite matrix. - `n` : [int] Degrees of freedom (n > 0). - - `Tau` : Symmetric and positive definite matrix. + - `C` : Symmetric and positive definite scale matrix .. note:: - Step method MatrixMetropolis will preserve the symmetry of Wishart variables. + Step method MatrixMetropolis will preserve the symmetry of + Wishart variables. """ - return flib.blas_inv_wishart(X,n,Tau) + return flib.blas_inv_wishart(X, n, C) + +def rinverse_wishart_prec(n, Tau): + """ + rinverse_wishart_prec(n, Tau) + + Return an inverse Wishart random matrix. + + n is the degrees of freedom. + Tau is a positive definite precision matrix + """ + wi = rwishart(n, np.asmatrix(Tau)).I + flib.symmetrize(wi) + return wi + +def inverse_wishart_prec_expval(X, n, Tau): + """ + inverse_wishart_expval(n, Tau) + + :Parameters: + - `n` : [int] Degrees of freedom (n > 0). + - `Tau` : Symmetric and positive definite precision matrix + + Expected value of inverse Wishart distribution. + """ + return inverse_wishart_like(X, n, inverse(Tau)) + +def inverse_wishart_prec_like(X, n, Tau): + """ + inverse_wishart_prec_like(X, n, Tau) + + Inverse Wishart log-likelihood + + For an alternative parameterization based on :math:`C=Tau^{-1}`, see + `inverse_wishart_like`. + + :Parameters: + - `X` : Symmetric, positive definite matrix. + - `n` : [int] Degrees of freedom (n > 0). + - `Tau` : Symmetric and positive definite precision matrix + """ + return inverse_wishart_like(X, n, inverse(Tau)) # Double exponential (Laplace)-------------------------------------------- @randomwrap @@ -1553,12 +1623,14 @@ def lognormal_like(x, mu, tau): # Multinomial---------------------------------------------- #@randomwrap -def rmultinomial(n,p,size=None): # Leaving size=None as the default means return value is 1d array if not specified-- nicer. +def rmultinomial(n,p,size=None): """ rmultinomial(n,p,size=1) Random multinomial variates. """ + # Leaving size=None as the default means return value is 1d array + # if not specified-- nicer. # Single value for p: if len(np.shape(p))==1: @@ -1595,11 +1667,11 @@ def multinomial_like(x, n, p): :Parameters: x : (ns, k) int - Random variable indicating the number of time outcome i is + Random variable indicating the number of time outcome i is observed. :math:`\sum_{i=1}^k x_i=n`, :math:`x_i \ge 0`. n : int Number of trials. - p : (k,) + p : (k,) Probability of each one of the different outcomes. :math:`\sum_{i=1}^k p_i = 1)`, :math:`p_i \ge 0`. @@ -1607,12 +1679,13 @@ def multinomial_like(x, n, p): - :math:`E(X_i)=n p_i` - :math:`Var(X_i)=n p_i(1-p_i)` - :math:`Cov(X_i,X_j) = -n p_i p_j` - - If :math: `\sum_i p_i < 0.999999` a log-likelihood value of -inf + - If :math: `\sum_i p_i < 0.999999` a log-likelihood value of -inf will be returned. """ - - x = np.atleast_2d(x) #flib expects 2d arguments. Do we still want to support multiple p values along realizations ? + # flib expects 2d arguments. Do we still want to support multiple p + # values along realizations ? + x = np.atleast_2d(x) p = np.atleast_2d(p) return flib.multinomial(x, n, p) @@ -1630,9 +1703,11 @@ def rmultivariate_hypergeometric(n, m, size=None): urn = np.repeat(np.arange(N), m) if size: - draw = np.array([[urn[i] for i in np.random.permutation(len(urn))[:n]] for j in range(size)]) + draw = np.array([[urn[i] for i in np.random.permutation(len(urn))[:n]] + for j in range(size)]) - r = [[np.sum(draw[j]==i) for i in range(len(m))] for j in range(size)] + r = [[np.sum(draw[j]==i) for i in range(len(m))] + for j in range(size)] else: draw = np.array([urn[i] for i in np.random.permutation(len(urn))[:n]]) @@ -1763,7 +1838,7 @@ def mv_normal_cov_like(x, mu, C): R""" mv_normal_cov_like(x, mu, C) - Multivariate normal log-likelihood parameterized by a covariance + Multivariate normal log-likelihood parameterized by a covariance matrix. .. math:: @@ -1772,7 +1847,7 @@ def mv_normal_cov_like(x, mu, C): :Parameters: - `x` : (n,k) - `mu` : (k) Location parameter. - - `C` : (k,k) Positive definite covariance matrix. + - `C` : (k,k) Positive definite covariance matrix. .. seealso:: :func:`mv_normal_like`, :func:`mv_normal_chol_like` """ @@ -1874,9 +1949,10 @@ def negative_binomial_like(x, mu, alpha): R""" negative_binomial_like(x, mu, alpha) - Negative binomial log-likelihood. The negative binomial distribution describes a - Poisson random variable whose rate parameter is gamma distributed. PyMC's chosen - parameterization is based on this mixture interpretation. + Negative binomial log-likelihood. The negative binomial + distribution describes a Poisson random variable whose rate + parameter is gamma distributed. PyMC's chosen parameterization is + based on this mixture interpretation. .. math:: f(x \mid \mu, \alpha) = \frac{\Gamma(x+\alpha)}{x! \Gamma(\alpha)} (\alpha/(\mu+\alpha))^\alpha (\mu/(\mu+\alpha))^x @@ -1926,7 +2002,8 @@ def normal_like(x, mu, tau): :Parameters: - `x` : Input data. - `mu` : Mean of the distribution. - - `tau` : Precision of the distribution, which corresponds to :math:`1/\sigma^2` (tau > 0). + - `tau` : Precision of the distribution, which corresponds to + :math:`1/\sigma^2` (tau > 0). .. note:: - :math:`E(X) = \mu` @@ -1937,7 +2014,7 @@ def normal_like(x, mu, tau): # constrain(tau, lower=0) # except ZeroProbability: # return -np.Inf - + return flib.normal(x, mu, tau) # von Mises-------------------------------------------------- @@ -1949,7 +2026,7 @@ def rvon_mises(mu, kappa, size=None): Random von Mises variates. """ # TODO: Just return straight from numpy after release 1.3 - return (np.random.mtrand.vonmises( mu, kappa, size) + np.pi)%(2.*np.pi)-np.pi + return (np.random.mtrand.vonmises(mu, kappa, size) + np.pi)%(2.*np.pi)-np.pi def von_mises_expval(mu, kappa): """ @@ -2007,10 +2084,11 @@ def poisson_like(x,mu): R""" poisson_like(x,mu) - Poisson log-likelihood. The Poisson is a discrete probability distribution. - It is often used to model the number of events occurring in a fixed period of - time when the times at which events occur are independent. The Poisson - distribution can be derived as a limiting case of the binomial distribution. + Poisson log-likelihood. The Poisson is a discrete probability + distribution. It is often used to model the number of events + occurring in a fixed period of time when the times at which events + occur are independent. The Poisson distribution can be derived as + a limiting case of the binomial distribution. .. math:: f(x \mid \mu) = \frac{e^{-\mu}\mu^x}{x!} @@ -2047,7 +2125,8 @@ def rtruncated_poisson(mu, k, size=None): except TypeError, ValueError: # More than one mu k=np.array(k)-1 - return np.transpose(np.array([rtruncated_poisson(x, i, size) for x,i in zip(mu, np.resize(k, np.size(mu)))])) + return np.array([rtruncated_poisson(x, i, size) + for x,i in zip(mu, np.resize(k, np.size(mu)))]).T # Calculate constant for acceptance probability C = np.exp(flib.factln(k+1)-flib.factln(k+1-m)) @@ -2060,10 +2139,11 @@ def rtruncated_poisson(mu, k, size=None): while(len(rvs) k # Uniform random variates @@ -2088,10 +2168,11 @@ def truncated_poisson_like(x,mu,k): R""" truncpoisson_like(x,mu,k) - Truncated Poisson log-likelihood. The Truncated Poisson is a discrete - probability distribution that is arbitrarily truncated to be greater than some - minimum value k. For example, zero-truncated Poisson distributions can be used - to model counts that are constrained to be non-negative. + Truncated Poisson log-likelihood. The Truncated Poisson is a + discrete probability distribution that is arbitrarily truncated to + be greater than some minimum value k. For example, zero-truncated + Poisson distributions can be used to model counts that are + constrained to be non-negative. .. math:: f(x \mid \mu, k) = \frac{e^{-\mu}\mu^x}{x!(1-F(k|\mu))} @@ -2134,7 +2215,13 @@ def truncated_normal_expval(mu, tau, a, b): """Expectation value of the truncated normal distribution. .. math:: - E(X)=\mu + \frac{\sigma(\varphi_1-\varphi_2)}{T}, where T=\Phi\left(\frac{B-\mu}{\sigma}\right)-\Phi\left(\frac{A-\mu}{\sigma}\right) and \varphi_1 = \varphi\left(\frac{A-\mu}{\sigma}\right) and \varphi_2 = \varphi\left(\frac{B-\mu}{\sigma}\right), where \varphi is the probability density function of a standard normal random variable and tau is 1/sigma**2. """ + E(X)=\mu + \frac{\sigma(\varphi_1-\varphi_2)}{T}, where + T=\Phi\left(\frac{B-\mu}{\sigma}\right)-\Phi + \left(\frac{A-\mu}{\sigma}\right) and \varphi_1 = + \varphi\left(\frac{A-\mu}{\sigma}\right) and \varphi_2 = + \varphi\left(\frac{B-\mu}{\sigma}\right), where \varphi is the + probability density function of a standard normal random + variable and tau is 1/sigma**2.""" phia = np.exp(normal_like(a, mu, tau)) phib = np.exp(normal_like(b, mu, tau)) sigma = 1./np.sqrt(tau) @@ -2246,9 +2333,10 @@ def rt(nu, size=None): def t_like(x, nu): R"""t_like(x, nu) - Student's T log-likelihood. Describes a zero-mean normal variable whose precision is - gamma distributed. Alternatively, describes the mean of several zero-mean normal - random variables divided by their sample standard deviation. + Student's T log-likelihood. Describes a zero-mean normal variable + whose precision is gamma distributed. Alternatively, describes the + mean of several zero-mean normal random variables divided by their + sample standard deviation. .. math:: f(x \mid \nu) = \frac{\Gamma(\frac{\nu+1}{2})}{\Gamma(\frac{\nu}{2}) \sqrt{\nu\pi}} \left( 1 + \frac{x^2}{\nu} \right)^{-\frac{\nu+1}{2}} @@ -2381,6 +2469,7 @@ def weibull_like(x, alpha, beta): return flib.weibull(x, alpha, beta) # Wishart--------------------------------------------------- + def rwishart(n, Tau): """ rwishart(n, Tau) @@ -2389,23 +2478,20 @@ def rwishart(n, Tau): Tau is the inverse of the 'covariance' matrix :math:`C`. """ - p = np.shape(Tau)[0] - # sig = np.linalg.cholesky(np.linalg.inv(Tau)) sig = np.linalg.cholesky(Tau) if n