From 760664d243000246af5a6ad447ec7157294493ed Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 11:46:03 -0400 Subject: [PATCH 01/11] diff cruft, and rinverse_wishart_cov function --- pymc/distributions.py | 45 ++++++++++++++++++++++++++++--------------- 1 file changed, 29 insertions(+), 16 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index d8dce607e9..c27fd6a697 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -123,7 +123,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 +254,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): """ @@ -964,17 +964,17 @@ 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. - + """ x = np.atleast_2d(x) @@ -1377,6 +1377,19 @@ def rinverse_wishart(n, Tau): flib.symmetrize(wi) return wi +def rinverse_wishart_cov(n, C): + """ + rinverse_wishart_cov(n, C) + + Return an inverse Wishart random matrix. + + n is the degrees of freedom. + C is a positive definite covariance matrix + """ + wi = rwishart(n, np.asmatrix(C)).I + flib.symmetrize(wi) + return wi + def inverse_wishart_expval(n, Tau): """ inverse_wishart_expval(n, Tau) @@ -1595,11 +1608,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,7 +1620,7 @@ 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. """ @@ -1763,7 +1776,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 +1785,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` """ @@ -1937,7 +1950,7 @@ def normal_like(x, mu, tau): # constrain(tau, lower=0) # except ZeroProbability: # return -np.Inf - + return flib.normal(x, mu, tau) # von Mises-------------------------------------------------- @@ -2815,13 +2828,13 @@ def Impute(name, dist_class, imputable, **parents): - parents (optional): dict Arbitrary keyword arguments. """ - + dims = np.shape(imputable) masked_values = np.ravel(imputable) if not type(masked_values) == np.ma.core.MaskedArray: # Generate mask - + mask = [v is None or np.isnan(v) for v in masked_values] # Generate masked array masked_values = np.ma.masked_array(masked_values, mask) From 85a77302b38c910e9a57672472607f17cae9f0d6 Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 18:53:09 -0400 Subject: [PATCH 02/11] InverseWishartCov implementation (need to test), code reformatting, consolidated 'namespace injection' code --- pymc/distributions.py | 305 +++++++++++++++++++++++++++--------------- 1 file changed, 200 insertions(+), 105 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index c27fd6a697..0cb903f399 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_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' 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): @@ -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)} @@ -972,10 +994,8 @@ def dirichlet_like(x, theta): 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. - - + 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,6 +1384,7 @@ def inverse_gamma_like(x, alpha, beta): return flib.igamma(x, alpha, beta) # Inverse Wishart--------------------------------------------------- + def rinverse_wishart(n, Tau): """ rinverse_wishart(n, Tau) @@ -1377,23 +1398,14 @@ def rinverse_wishart(n, Tau): flib.symmetrize(wi) return wi -def rinverse_wishart_cov(n, C): - """ - rinverse_wishart_cov(n, C) - - Return an inverse Wishart random matrix. - - n is the degrees of freedom. - C is a positive definite covariance matrix - """ - wi = rwishart(n, np.asmatrix(C)).I - flib.symmetrize(wi) - return wi - def inverse_wishart_expval(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 np.asarray(Tau)/(n-len(Tau)-1) @@ -1402,8 +1414,9 @@ def inverse_wishart_like(X, n, Tau): R""" inverse_wishart_like(X, n, Tau) - 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)} @@ -1413,14 +1426,62 @@ def inverse_wishart_like(X, n, Tau): :Parameters: - `X` : Symmetric, positive definite matrix. - `n` : [int] Degrees of freedom (n > 0). - - `Tau` : Symmetric and positive definite matrix. + - `Tau` : Symmetric and positive definite precision 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) +def rinverse_wishart_cov(n, C): + """ + rinverse_wishart_cov(n, C) + + Return an inverse Wishart random matrix. + + n is the degrees of freedom. + C is a positive definite covariance matrix + """ + wi = rwishart(n, np.asmatrix(C)).I + flib.symmetrize(wi) + return wi + +def inverse_wishart_cov_expval(X, n, C): + """ + inverse_wishart_expval(n, Tau) + + For an alternative parameterization based on :math:`T=C^{-1}`, see + `inverse_wishart_expval`. + + :Parameters: + - `n` : [int] Degrees of freedom (n > 0). + - `Tau` : Symmetric and positive definite precision matrix + + Expected value of inverse Wishart distribution. + :Parameters: + - `n` : [int] Degrees of freedom (n > 0). + - `C` : Symmetric and positive definite covariance matrix + """ + return inverse_wishart_like(X, n, inverse(C)) + +def inverse_wishart_cov_like(X, n, C): + """ + inverse_wishart_cov_like(X, n, C) + + Inverse Wishart log-likelihood + + For an alternative parameterization based on :math:`T=C^{-1}`, see + `inverse_wishart_like`. + + :Parameters: + - `X` : Symmetric, positive definite matrix. + - `n` : [int] Degrees of freedom (n > 0). + - `C` : Symmetric and positive definite covariance matrix + """ + return inverse_wishart_like(X, n, inverse(C)) + # Double exponential (Laplace)-------------------------------------------- @randomwrap def rlaplace(mu, tau, size=None): @@ -1566,7 +1627,9 @@ 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): + # Leaving size=None as the default means return value is 1d array + # if not specified-- nicer. """ rmultinomial(n,p,size=1) @@ -1887,9 +1950,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 @@ -2020,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!} @@ -2101,10 +2166,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))} @@ -2259,9 +2325,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}} @@ -2495,7 +2562,7 @@ def wishart_cov_like(X, n, C): matrix of a multivariate normal distribution. If C=1, the distribution is identical to the chi-square distribution with n degrees of freedom. - For an alternative parameterization based on :math:`T=C{-1}`, see + For an alternative parameterization based on :math:`T=C^{-1}`, see `wishart_like`. .. math:: @@ -2569,33 +2636,35 @@ def local_decorated_likelihoods(obj): obj[name+'_like'] = gofwrapper(like, snapshot) -#local_decorated_likelihoods(locals()) -# Decorating the likelihoods breaks the creation of distribution instantiators -DH. - - +# local_decorated_likelihoods(locals()) +# Decorating the likelihoods breaks the creation of distribution +# instantiators -DH. # Create Stochastic instantiators +def _inject_dist(distname, kwargs={}, ns=locals()): + """ + Reusable function to inject Stochastic subclasses into module + namespace + """ + dist_logp, dist_random = name_to_funcs(dist, ns) + classname = capitalize(dist) + ns[classname]= stochastic_from_dist(dist, dist_logp, + dist_random, **kwargs) + for dist in sc_continuous_distributions: - dist_logp, dist_random = name_to_funcs(dist, locals()) - locals()[capitalize(dist)]= stochastic_from_dist(dist, dist_logp, dist_random) + _inject_dist(dist) for dist in mv_continuous_distributions: - dist_logp, dist_random = name_to_funcs(dist, locals()) - locals()[capitalize(dist)]= stochastic_from_dist(dist, dist_logp, dist_random, mv=True) + _inject_dist(dist, kwargs={'mv' : True}) for dist in sc_discrete_distributions: - dist_logp, dist_random = name_to_funcs(dist, locals()) - locals()[capitalize(dist)]= stochastic_from_dist(dist, dist_logp, dist_random, dtype=np.int) + _inject_dist(dist, kwargs={'dtype' : np.int}) for dist in mv_discrete_distributions: - dist_logp, dist_random = name_to_funcs(dist, locals()) - locals()[capitalize(dist)]= stochastic_from_dist(dist, dist_logp, dist_random, dtype=np.int, mv=True) - - -dist_logp, dist_random = name_to_funcs('bernoulli', locals()) -Bernoulli = stochastic_from_dist('bernoulli', dist_logp, dist_random, dtype=np.bool) + _inject_dist(dist, kwargs={'dtype' : np.int, 'mv' : True}) +_inject_dist('bernoulli', {'dtype' : np.bool}) def uninformative_like(x): """ @@ -2618,14 +2687,13 @@ def one_over_x_like(x): return -np.sum(np.log(x)) -Uninformative = stochastic_from_dist('uninformative', logp = uninformative_like) -DiscreteUninformative = stochastic_from_dist('uninformative', logp = uninformative_like, dtype=np.int) +Uninformative = stochastic_from_dist('uninformative', logp=uninformative_like) +DiscreteUninformative = stochastic_from_dist('uninformative', logp=uninformative_like, dtype=np.int) DiscreteUninformative.__name__ = 'DiscreteUninformative' OneOverX = stochastic_from_dist('one_over_x', logp = one_over_x_like) - - -# Conjugates of Dirichlet get special treatment, can be parametrized by first k-1 'p' values +# Conjugates of Dirichlet get special treatment, can be parametrized +# by first k-1 'p' values def extend_dirichlet(p): """ @@ -2692,7 +2760,8 @@ def rmod_categor(p,size=None): class Categorical(Stochastic): __doc__ = """ -C = Categorical(name, p, value=None, dtype=np.int, observed=False, size=1, trace=True, rseed=False, cache_depth=2, plot=None) +C = Categorical(name, p, value=None, dtype=np.int, observed=False, +size=1, trace=True, rseed=False, cache_depth=2, plot=None) Stochastic variable with Categorical distribution. Parent is: p @@ -2713,7 +2782,9 @@ class Categorical(Stochastic): parent_names = ['p', 'minval', 'step'] - def __init__(self, name, p, value=None, dtype=np.int, observed=False, size=1, trace=True, rseed=False, cache_depth=2, plot=None, verbose=None,**kwds): + def __init__(self, name, p, value=None, dtype=np.int, observed=False, + size=1, trace=True, rseed=False, cache_depth=2, plot=None, + verbose=None,**kwds): if value is not None: if np.isscalar(value): @@ -2724,13 +2795,22 @@ def __init__(self, name, p, value=None, dtype=np.int, observed=False, size=1, tr self.size = size if isinstance(p, Dirichlet): - Stochastic.__init__(self, logp=valuewrapper(mod_categorical_like), doc='A Categorical random variable', name=name, - parents={'p':p}, random=bind_size(rmod_categor, self.size), trace=trace, value=value, dtype=dtype, - rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) + Stochastic.__init__(self, logp=valuewrapper(mod_categorical_like), + doc='A Categorical random variable', name=name, + parents={'p':p}, random=bind_size(rmod_categor, self.size), + trace=trace, value=value, dtype=dtype, + rseed=rseed, observed=observed, + cache_depth=cache_depth, plot=plot, + verbose=verbose, **kwds) else: - Stochastic.__init__(self, logp=valuewrapper(categorical_like), doc='A Categorical random variable', name=name, - parents={'p':p}, random=bind_size(rcategorical, self.size), trace=trace, value=value, dtype=dtype, - rseed=rseed, observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) + Stochastic.__init__(self, logp=valuewrapper(categorical_like), + doc='A Categorical random variable', name=name, + parents={'p':p}, + random=bind_size(rcategorical, self.size), + trace=trace, value=value, dtype=dtype, + rseed=rseed, observed=observed, + cache_depth=cache_depth, plot=plot, + verbose=verbose, **kwds) # class ModCategorical(Stochastic): # __doc__ = """ @@ -2797,16 +2877,26 @@ class Multinomial(Stochastic): parent_names = ['n', 'p'] - def __init__(self, name, n, p, trace=True, value=None, rseed=False, observed=False, cache_depth=2, plot=None, verbose=None, **kwds): + def __init__(self, name, n, p, trace=True, value=None, rseed=False, + observed=False, cache_depth=2, plot=None, verbose=None, + **kwds): if isinstance(p, Dirichlet): - Stochastic.__init__(self, logp=valuewrapper(mod_multinom_like), doc='A Multinomial random variable', name=name, - parents={'n':n,'p':p}, random=mod_rmultinom, trace=trace, value=value, dtype=np.int, rseed=rseed, - observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) + Stochastic.__init__(self, logp=valuewrapper(mod_multinom_like), + doc='A Multinomial random variable', name=name, + parents={'n':n,'p':p}, random=mod_rmultinom, + trace=trace, value=value, dtype=np.int, + rseed=rseed, observed=observed, + cache_depth=cache_depth, plot=plot, + verbose=verbose, **kwds) else: - Stochastic.__init__(self, logp=valuewrapper(multinomial_like), doc='A Multinomial random variable', name=name, - parents={'n':n,'p':p}, random=rmultinomial, trace=trace, value=value, dtype=np.int, rseed=rseed, - observed=observed, cache_depth=cache_depth, plot=plot, verbose=verbose, **kwds) + Stochastic.__init__(self, logp=valuewrapper(multinomial_like), + doc='A Multinomial random variable', name=name, + parents={'n':n,'p':p}, random=rmultinomial, + trace=trace, value=value, dtype=np.int, + rseed=rseed, observed=observed, + cache_depth=cache_depth, plot=plot, + verbose=verbose, **kwds) def Impute(name, dist_class, imputable, **parents): """ @@ -2823,8 +2913,9 @@ def Impute(name, dist_class, imputable, **parents): - dist_class : Stochastic Stochastic subclass such as Poisson, Normal, etc. - imputable : numpy.ma.core.MaskedArray or iterable - A masked array with missing elements (where mask=True, value is assumed missing), - or any iterable that contains None elements that will be imputed. + A masked array with missing elements (where mask=True, value + is assumed missing), or any iterable that contains None + elements that will be imputed. - parents (optional): dict Arbitrary keyword arguments. """ @@ -2857,9 +2948,12 @@ def Impute(name, dist_class, imputable, **parents): shape = np.shape(parent) if shape == dims: - these_parents[key] = Lambda(key + '[%i]'%i, lambda p=np.ravel(parent), i=i: p[i]) + these_parents[key] = Lambda(key + '[%i]'%i, + lambda p=np.ravel(parent), + i=i: p[i]) elif shape == np.shape(masked_values): - these_parents[key] = Lambda(key + '[%i]'%i, lambda p=parent, i=i: p[i]) + these_parents[key] = Lambda(key + '[%i]'%i, lambda p=parent, + i=i: p[i]) else: these_parents[key] = parent @@ -2868,7 +2962,8 @@ def Impute(name, dist_class, imputable, **parents): vars.append(dist_class(this_name, **these_parents)) else: # Observed values - vars.append(dist_class(this_name, value=masked_values[i], observed=True, **these_parents)) + vars.append(dist_class(this_name, value=masked_values[i], + observed=True, **these_parents)) return np.reshape(vars, dims) ImputeMissing = Impute From 3b6d3544026d0f373f2a6747f117c96e9a8253a0 Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 18:55:37 -0400 Subject: [PATCH 03/11] fixed messed up docstring with comment --- pymc/distributions.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index 0cb903f399..4504cdbc36 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -1628,13 +1628,13 @@ 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. """ 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: From 3f163d669deb523d554b4d04ab3124de76cffd66 Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 21:34:01 -0400 Subject: [PATCH 04/11] Changed rwishart_cov to use rwishart pending review of algorithm. rwishart and rwishart_cov were producing different results for the same random seed (the random variables inside are generated in the same order) --- pymc/distributions.py | 47 ++++++++++++++++++++++++++----------------- 1 file changed, 28 insertions(+), 19 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index 4504cdbc36..23d6f64acc 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -2461,6 +2461,7 @@ def weibull_like(x, alpha, beta): return flib.weibull(x, alpha, beta) # Wishart--------------------------------------------------- + def rwishart(n, Tau): """ rwishart(n, Tau) @@ -2469,23 +2470,21 @@ 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 Date: Fri, 22 Oct 2010 22:15:56 -0400 Subject: [PATCH 05/11] * Changed inverse_wishart docstrings and function signatures to clarify use of covariance (scale) instead of precision * Implemented inverse wishart distribution parameterized using precision matrix * More unit tests to check consistency of cov / precision based wishart / inverse wishart-based samples / likelihoods * A little bit of code reformatting and fixed earlier M-q reformats --- pymc/distributions.py | 98 ++++++++++++++++-------------- pymc/tests/test_distributions.py | 101 +++++++++++++++++++++---------- pymc/utils.py | 28 +++++---- 3 files changed, 139 insertions(+), 88 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index 23d6f64acc..cafd76d182 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -74,7 +74,7 @@ class ArgumentError(AttributeError): mv_continuous_distributions = ['dirichlet', 'inverse_wishart', 'mv_normal', 'mv_normal_cov', 'mv_normal_chol', 'wishart', - 'wishart_cov', 'inverse_wishart_cov'] + 'wishart_cov', 'inverse_wishart_prec'] mv_discrete_distributions = ['multivariate_hypergeometric', 'multinomial'] @@ -676,8 +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} @@ -774,8 +774,8 @@ def betabin_like(x, alpha, beta, n): 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. + 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)} @@ -1385,102 +1385,98 @@ def inverse_gamma_like(x, alpha, beta): # Inverse Wishart--------------------------------------------------- -def rinverse_wishart(n, Tau): +def rinverse_wishart(n, Sigma): """ - rinverse_wishart(n, Tau) + rinverse_wishart(n, Sigma) Return an inverse Wishart random matrix. n is the degrees of freedom. - Tau is a positive definite scale matrix. + Sigma is a positive definite scale matrix. """ - wi = rwishart(n, np.asmatrix(Tau).I).I + wi = rwishart(n, np.asmatrix(Sigma).I).I flib.symmetrize(wi) return wi -def inverse_wishart_expval(n, Tau): +def inverse_wishart_expval(n, Sigma): """ - inverse_wishart_expval(n, Tau) + inverse_wishart_expval(n, Sigma) :Parameters: - `n` : [int] Degrees of freedom (n > 0). - - `Tau` : Symmetric and positive definite precision matrix + - `Sigma` : Symmetric and positive definite scale matrix Expected value of inverse Wishart distribution. """ - return np.asarray(Tau)/(n-len(Tau)-1) + return np.asarray(Sigma)/(n-len(Sigma)-1) -def inverse_wishart_like(X, n, Tau): +def inverse_wishart_like(X, n, Sigma): R""" - inverse_wishart_like(X, n, Tau) + inverse_wishart_like(X, n, Sigma) 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 precision matrix + - `Sigma` : Symmetric and positive definite scale matrix .. note:: 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, Sigma) -def rinverse_wishart_cov(n, C): +def rinverse_wishart_prec(n, Tau): """ - rinverse_wishart_cov(n, C) + rinverse_wishart_prec(n, Tau) Return an inverse Wishart random matrix. n is the degrees of freedom. - C is a positive definite covariance matrix + Tau is a positive definite precision matrix """ - wi = rwishart(n, np.asmatrix(C)).I + wi = rwishart(n, np.asmatrix(Tau)).I flib.symmetrize(wi) return wi -def inverse_wishart_cov_expval(X, n, C): +def inverse_wishart_prec_expval(X, n, Tau): """ inverse_wishart_expval(n, Tau) - For an alternative parameterization based on :math:`T=C^{-1}`, see - `inverse_wishart_expval`. - :Parameters: - `n` : [int] Degrees of freedom (n > 0). - `Tau` : Symmetric and positive definite precision matrix Expected value of inverse Wishart distribution. - :Parameters: - - `n` : [int] Degrees of freedom (n > 0). - - `C` : Symmetric and positive definite covariance matrix """ - return inverse_wishart_like(X, n, inverse(C)) + return inverse_wishart_like(X, n, inverse(Tau)) -def inverse_wishart_cov_like(X, n, C): +def inverse_wishart_prec_like(X, n, Tau): """ - inverse_wishart_cov_like(X, n, C) + inverse_wishart_prec_like(X, n, Tau) Inverse Wishart log-likelihood - For an alternative parameterization based on :math:`T=C^{-1}`, see + 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). - - `C` : Symmetric and positive definite covariance matrix + - `Tau` : Symmetric and positive definite precision matrix """ - return inverse_wishart_like(X, n, inverse(C)) + return inverse_wishart_like(X, n, inverse(Tau)) # Double exponential (Laplace)-------------------------------------------- @randomwrap @@ -1687,8 +1683,9 @@ def multinomial_like(x, n, p): 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) @@ -1706,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]]) @@ -2003,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` @@ -2026,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): """ @@ -2125,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)) @@ -2138,10 +2139,11 @@ def rtruncated_poisson(mu, k, size=None): while(len(rvs) k # Uniform random variates @@ -2213,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) diff --git a/pymc/tests/test_distributions.py b/pymc/tests/test_distributions.py index a3d05e32b7..b40e62c2e2 100644 --- a/pymc/tests/test_distributions.py +++ b/pymc/tests/test_distributions.py @@ -16,17 +16,20 @@ # FIXME no tests for discrete_uniform, negative_binomial, uniform. # TODO add tests for boundaries of distributions with restricted support - #from decorators import * -from pymc.distributions import * + +from unittest import TestCase import unittest +import os, pdb, warnings, nose + +from numpy import exp, log, array, sqrt +from numpy.linalg import cholesky, det, inv from numpy.testing import * -from pymc import flib, utils import numpy as np -from numpy import exp, log, array, sqrt -from numpy.linalg import cholesky -import os, pdb, warnings, nose -from unittest import TestCase +import numpy.random as _npr + +from pymc.distributions import * +from pymc import flib, utils PLOT=True DIR = 'testresults' @@ -42,13 +45,13 @@ from scipy.optimize import fmin SP = True except: - SP= False + SP = False try: import pylab as P except: print 'Plotting disabled' - PLOT=False + PLOT = False # Some python densities for comparison @@ -113,7 +116,7 @@ def mv_normal(x, mu, C): mu = np.asmatrix(mu) C = np.asmatrix(C) - I = (N/2.)*log(2.*pi) + .5*log(np.linalg.det(C)) + I = (N/2.)*log(2.*pi) + .5*log(det(C)) z = (x-mu) return -(I +.5 * z * np.linalg.inv(C) * z.T).A[0][0] @@ -389,7 +392,7 @@ def test_consistency(self): compare_hist(figname='categorical', **figdata) # test_normalization assert_almost_equal(like.sum(), 1, 4) - + def test_support(self): """Check to make sure values outside support are -inf""" assert categorical_like([-1], [[0.4,0.4,0.2]]) < -1e300 @@ -914,14 +917,30 @@ def test_calling(self): b = flib.weibull([1,2], [2,2], [3,3]) assert_equal(a,b) +#------------------------------------------------------------------------------- +# Test Wishart / Inverse Wishart distributions +_Tau_test = np.matrix([[ 209.47883244, 10.88057915, 13.80581557], + [ 10.88057915, 213.58694978, 11.18453854], + [ 13.80581557, 11.18453854, 209.89396417]]) + class test_wishart(TestCase): """ - There are results out there on the asymptotic distribution of eigenvalues - of very large Wishart matrices that we could use to make another test case... + There are results out there on the asymptotic distribution of + eigenvalues of very large Wishart matrices that we could use to + make another test case... """ - Tau_test = np.matrix([[ 209.47883244, 10.88057915, 13.80581557], - [ 10.88057915, 213.58694978, 11.18453854], - [ 13.80581557, 11.18453854, 209.89396417]])/100. + # precision matrix + Tau_test = _Tau_test / 100 + + def test_samples(self): + # test consistency between precision and cov-based + _npr.seed(1) + sample_a = rwishart(100, self.Tau_test) + + _npr.seed(1) + sample_b = rwishart_cov(100, self.Tau_test.I) + + assert_array_almost_equal(sample_a, sample_b) def test_likelihoods(self): try: @@ -934,7 +953,8 @@ def test_likelihoods(self): def slo_wishart_cov(W,n,V): p = W.shape[0] - logp = (n-p-1)*.5 * np.log(np.linalg.det(W)) - n*p*.5*np.log(2) - n*.5*np.log(np.linalg.det(V)) - p*(p-1)/4.*np.log(pi) + logp = ((n-p-1)*.5 * np.log(det(W)) - n*p*.5*np.log(2) + - n*.5*np.log(det(V)) - p*(p-1)/4.*np.log(pi)) for i in xrange(1,p+1): logp -= gammaln((n+1-i)*.5) logp -= np.trace(np.linalg.solve(V,W)*.5) @@ -942,8 +962,11 @@ def slo_wishart_cov(W,n,V): for i in [5,10,100,10000]: right_answer = slo_wishart_cov(W_test,i,self.Tau_test.I) - assert_array_almost_equal(wishart_like(W_test,i,self.Tau_test), right_answer, decimal=5) - assert_array_almost_equal(wishart_cov_like(W_test,i,self.Tau_test.I), right_answer, decimal=5) + assert_array_almost_equal(wishart_like(W_test,i,self.Tau_test), + right_answer, decimal=5) + assert_array_almost_equal( + wishart_cov_like(W_test,i,self.Tau_test.I), + right_answer, decimal=5) def test_expval(self): @@ -969,9 +992,18 @@ class test_inverse_wishart(TestCase): """ Adapted from test_wishart """ - Tau_test = np.matrix([[ 209.47883244, 10.88057915, 13.80581557], - [ 10.88057915, 213.58694978, 11.18453854], - [ 13.80581557, 11.18453854, 209.89396417]]).I + # Covariance matrix (!) + Sigma_test = _Tau_test.I + + def test_samples(self): + # test consistency between precision and cov-based + _npr.seed(1) + sample_a = rinverse_wishart(100, self.Sigma_test) + + _npr.seed(1) + sample_b = rinverse_wishart_prec(100, self.Sigma_test.I) + + assert_array_almost_equal(sample_a, sample_b) def test_likelihoods(self): try: @@ -979,12 +1011,13 @@ def test_likelihoods(self): except: raise nose.SkipTest, "SciPy not installed." - IW_test = rinverse_wishart(100,self.Tau_test) + IW_test = rinverse_wishart(100,self.Sigma_test) def slo_inv_wishart(W,n,V): p = W.shape[0] - logp = -0.5*(n+p+1) * np.log(np.linalg.det(W)) - n*p*.5*np.log(2) + n*.5*np.log(np.linalg.det(V)) - p*(p-1)/4.*np.log(pi) + logp = (-0.5*(n+p+1) * np.log(det(W)) - n*p*.5*np.log(2) + + n*.5*np.log(det(V)) - p*(p-1)/4.*np.log(pi)) for i in xrange(1,p+1): logp -= gammaln((n+1-i)*.5) @@ -993,8 +1026,12 @@ def slo_inv_wishart(W,n,V): return logp for i in [5,10,100,10000]: - right_answer = slo_inv_wishart(IW_test,i,self.Tau_test) - assert_array_almost_equal(inverse_wishart_like(IW_test,i,self.Tau_test), right_answer, decimal=1) + right_answer = slo_inv_wishart(IW_test,i,self.Sigma_test) + calculated = inverse_wishart_like(IW_test,i,self.Sigma_test) + assert_array_almost_equal(calculated, right_answer, decimal=1) + + calculated = inverse_wishart_prec_like(IW_test,i,self.Sigma_test.I) + assert_array_almost_equal(calculated, right_answer, decimal=1) """ def test_expval(self): @@ -1002,11 +1039,11 @@ def test_expval(self): n = 100 N = 1000 - A = 0.*self.Tau_test + A = 0.*self.Sigma_test for i in xrange(N): - A += rinverse_wishart(n,self.Tau_test) + A += rinverse_wishart(n,self.Sigma_test) A /= N - delta=A-inverse_wishart_expval(n,self.Tau_test) + delta=A-inverse_wishart_expval(n,self.Sigma_test) print np.abs(np.asarray(delta)/np.asarray(A)).max() assert(np.abs(np.asarray(delta)/np.asarray(A)).max()<.5) """ @@ -1017,11 +1054,11 @@ def test_randomwrap(self): B = Bernoulli('x', np.ones(s)*.5) B.random() assert_equal(B.value.shape, s) - + N = Normal('x', np.ones(s), 1) N.random() assert_equal(N.value.shape, s) - + N = Normal('x', np.ones(s), np.ones(s)) N.random() assert_equal(N.value.shape, s) @@ -1044,7 +1081,7 @@ def means (mean = mean): #data has shape (10,5) but means returns an array of shape (10 * 5,) assert_raises(ValueError, pymc.Normal, "obs", mu = means, tau = 1, observed = True, value = data) - + if __name__ == '__main__': diff --git a/pymc/utils.py b/pymc/utils.py index 7a477f998b..d0b0e695ae 100644 --- a/pymc/utils.py +++ b/pymc/utils.py @@ -7,28 +7,34 @@ import numpy as np import sys, inspect, select, os, time from copy import copy -from PyMCObjects import Stochastic, Deterministic, Node, Variable, Potential, ZeroProbability +from PyMCObjects import (Stochastic, Deterministic, Node, Variable, Potential, + ZeroProbability) import flib import pdb from numpy.linalg.linalg import LinAlgError from numpy.linalg import cholesky, eigh, det, inv from Node import logp_of_set -from numpy import sqrt, obj2sctype, ndarray, asmatrix, array, pi, prod, exp,\ - pi, asarray, ones, atleast_1d, iterable, linspace, diff, around, log10, \ - zeros, arange, digitize, apply_along_axis, concatenate, bincount, sort, \ - hsplit, argsort, inf, shape, ndim, swapaxes, ravel, transpose as tr, \ - diag, cov +from numpy import (sqrt, obj2sctype, ndarray, asmatrix, array, pi, prod, exp, + pi, asarray, ones, atleast_1d, iterable, linspace, diff, + around, log10, zeros, arange, digitize, apply_along_axis, + concatenate, bincount, sort, hsplit, argsort, inf, shape, + ndim, swapaxes, ravel, transpose as tr, diag, cov) -__all__ = ['check_list', 'autocorr', 'calc_min_interval', 'check_type', 'ar1', 'ar1_gen', 'draw_random', 'histogram', 'hpd', 'invcdf', 'make_indices', 'normcdf', 'quantiles', 'rec_getattr', 'rec_setattr', 'round_array', 'trace_generator','msqrt','safe_len', 'log_difference', 'find_generations','crawl_dataless', 'logit', 'invlogit','stukel_logit','stukel_invlogit','symmetrize','value'] +__all__ = ['check_list', 'autocorr', 'calc_min_interval', 'check_type', 'ar1', + 'ar1_gen', 'draw_random', 'histogram', 'hpd', 'invcdf', + 'make_indices', 'normcdf', 'quantiles', 'rec_getattr', + 'rec_setattr', 'round_array', 'trace_generator','msqrt','safe_len', + 'log_difference', 'find_generations','crawl_dataless', 'logit', + 'invlogit','stukel_logit','stukel_invlogit','symmetrize','value'] -symmetrize=flib.symmetrize +symmetrize = flib.symmetrize def value(a): """ Returns a.value if a is a Variable, or just a otherwise. """ - if isinstance(a,Variable): + if isinstance(a, Variable): return a.value else: return a @@ -492,10 +498,10 @@ def autocorr(x, lag=1): # return ((x[:-lag]-mu)*(x[lag:]-mu)).sum()/v/(len(x) - lag) S = autocov(x, lag) return S[0,1]/sqrt(prod(diag(S))) - + def autocov(x, lag=1): """ - Sample autocovariance at specified lag. + Sample autocovariance at specified lag. The autocovariance is a 2x2 matrix with the variances of x[:-lag] and x[lag:] in the diagonal and the autocovariance on the off-diagonal. From cd5620a012120ebb98f74032cd0bab964070646c Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 22:18:17 -0400 Subject: [PATCH 06/11] Added precision-based inverse wishart to docs --- sphinxdoc/source/distributions.rst | 16 +++++++++------- 1 file changed, 9 insertions(+), 7 deletions(-) diff --git a/sphinxdoc/source/distributions.rst b/sphinxdoc/source/distributions.rst index 160d477771..7ffb993cd0 100644 --- a/sphinxdoc/source/distributions.rst +++ b/sphinxdoc/source/distributions.rst @@ -4,12 +4,14 @@ Probability distributions ************************* +PyMC provides 35 built-in probability distributions. For each +distribution, it provides: -PyMC provides 35 built-in probability distributions. For each distribution, it provides: - -* A function that evaluates its log-probability or log-density: normal_like(). +* A function that evaluates its log-probability or log-density: + normal_like(). * A function that draws random variables: rnormal(). -* A function that computes the expectation associated with the distribution: normal_expval(). +* A function that computes the expectation associated with the + distribution: normal_expval(). * A Stochastic subclass generated from the distribution: Normal. This section describes the likelihood functions of these distributions. @@ -20,7 +22,7 @@ Discrete distributions ====================== .. autofunction:: bernoulli_like - + .. autofunction:: binomial_like .. autofunction:: categorical_like @@ -94,6 +96,8 @@ Multivariate continuous distributions .. autofunction:: inverse_wishart_like +.. autofunction:: inverse_wishart_prec_like + .. autofunction:: mv_normal_like .. autofunction:: mv_normal_chol_like @@ -103,5 +107,3 @@ Multivariate continuous distributions .. autofunction:: wishart_like .. autofunction:: wishart_cov_like - - From 4e6987de99788a575a33bafee211d7693437f567 Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Fri, 22 Oct 2010 22:18:46 -0400 Subject: [PATCH 07/11] Have git ignore unit testing output and built sphinxdoc --- .gitignore | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) 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 From 86cb982ce9d118809a235d233484a07f0acae74a Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Thu, 4 Nov 2010 00:53:52 -0400 Subject: [PATCH 08/11] renamed Sigma to C --- pymc/distributions.py | 58 +++++++++++++++++++++---------------------- 1 file changed, 28 insertions(+), 30 deletions(-) diff --git a/pymc/distributions.py b/pymc/distributions.py index cafd76d182..e6c59da267 100755 --- a/pymc/distributions.py +++ b/pymc/distributions.py @@ -1385,34 +1385,34 @@ def inverse_gamma_like(x, alpha, beta): # Inverse Wishart--------------------------------------------------- -def rinverse_wishart(n, Sigma): +def rinverse_wishart(n, C): """ - rinverse_wishart(n, Sigma) + rinverse_wishart(n, C) Return an inverse Wishart random matrix. n is the degrees of freedom. - Sigma is a positive definite scale matrix. + C is a positive definite scale matrix. """ - wi = rwishart(n, np.asmatrix(Sigma).I).I + wi = rwishart(n, np.asmatrix(C).I).I flib.symmetrize(wi) return wi -def inverse_wishart_expval(n, Sigma): +def inverse_wishart_expval(n, C): """ - inverse_wishart_expval(n, Sigma) + inverse_wishart_expval(n, C) :Parameters: - `n` : [int] Degrees of freedom (n > 0). - - `Sigma` : Symmetric and positive definite scale matrix + - `C` : Symmetric and positive definite scale matrix Expected value of inverse Wishart distribution. """ - return np.asarray(Sigma)/(n-len(Sigma)-1) + return np.asarray(C)/(n-len(C)-1) -def inverse_wishart_like(X, n, Sigma): +def inverse_wishart_like(X, n, C): R""" - inverse_wishart_like(X, n, Sigma) + 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 @@ -1428,14 +1428,14 @@ def inverse_wishart_like(X, n, Sigma): :Parameters: - `X` : Symmetric, positive definite matrix. - `n` : [int] Degrees of freedom (n > 0). - - `Sigma` : Symmetric and positive definite scale matrix + - `C` : Symmetric and positive definite scale matrix .. note:: Step method MatrixMetropolis will preserve the symmetry of Wishart variables. """ - return flib.blas_inv_wishart(X, n, Sigma) + return flib.blas_inv_wishart(X, n, C) def rinverse_wishart_prec(n, Tau): """ @@ -2488,7 +2488,6 @@ def rwishart(n, Tau): A = flib.expand_triangular(chi_sqs, norms) flib.dtrsm_wrap(sig, A, side='L', uplo='L', transa='T') - # flib.dtrmm_wrap(sig,A,side='L',uplo='L',transa='N') w = np.asmatrix(np.dot(A,A.T)) flib.symmetrize(w) return w @@ -2538,29 +2537,28 @@ def rwishart_cov(n, C): """ rwishart(n, C) + n is degrees of freedom, C is the 'covariance' matrix + Return a Wishart random matrix. """ - return rwishart(n, np.linalg.inv(C)) - - # This code does *not* produce the same result as rwishart for - # identical random seed! + # return rwishart(n, np.linalg.inv(C)) - # p = np.shape(C)[0] - # # Need cholesky decomposition of precision matrix C^-1? - # sig = np.linalg.cholesky(C) + p = np.shape(C)[0] + # Need cholesky decomposition of precision matrix C^-1? + sig = np.linalg.cholesky(C) - # if n Date: Thu, 4 Nov 2010 00:54:42 -0400 Subject: [PATCH 09/11] commented out test_samples --- pymc/tests/test_distributions.py | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/pymc/tests/test_distributions.py b/pymc/tests/test_distributions.py index b40e62c2e2..9682720dbf 100644 --- a/pymc/tests/test_distributions.py +++ b/pymc/tests/test_distributions.py @@ -932,15 +932,15 @@ class test_wishart(TestCase): # precision matrix Tau_test = _Tau_test / 100 - def test_samples(self): - # test consistency between precision and cov-based - _npr.seed(1) - sample_a = rwishart(100, self.Tau_test) + # def test_samples(self): + # # test consistency between precision and cov-based + # _npr.seed(1) + # sample_a = rwishart(100, self.Tau_test) - _npr.seed(1) - sample_b = rwishart_cov(100, self.Tau_test.I) + # _npr.seed(1) + # sample_b = rwishart_cov(100, self.Tau_test.I) - assert_array_almost_equal(sample_a, sample_b) + # assert_array_almost_equal(sample_a, sample_b) def test_likelihoods(self): try: @@ -995,15 +995,15 @@ class test_inverse_wishart(TestCase): # Covariance matrix (!) Sigma_test = _Tau_test.I - def test_samples(self): - # test consistency between precision and cov-based - _npr.seed(1) - sample_a = rinverse_wishart(100, self.Sigma_test) + # def test_samples(self): + # # test consistency between precision and cov-based + # _npr.seed(1) + # sample_a = rinverse_wishart(100, self.Sigma_test) - _npr.seed(1) - sample_b = rinverse_wishart_prec(100, self.Sigma_test.I) + # _npr.seed(1) + # sample_b = rinverse_wishart_prec(100, self.Sigma_test.I) - assert_array_almost_equal(sample_a, sample_b) + # assert_array_almost_equal(sample_a, sample_b) def test_likelihoods(self): try: From 836d06ae949f40d88620391b59aea517ab788017 Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Thu, 4 Nov 2010 00:57:41 -0400 Subject: [PATCH 10/11] renamed Tau to C in InverseWishart --- pymc/tests/test_adaptive.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pymc/tests/test_adaptive.py b/pymc/tests/test_adaptive.py index 43e52e5ac8..e5d5664678 100644 --- a/pymc/tests/test_adaptive.py +++ b/pymc/tests/test_adaptive.py @@ -6,7 +6,7 @@ import nose, warnings def test_square(): - iw = pymc.InverseWishart("A", n = 2, Tau = np.eye(2)) + iw = pymc.InverseWishart("A", n = 2, C = np.eye(2)) mnc = pymc.MvNormalCov("v", mu = np.zeros(2), C = iw, value = np.zeros(2), observed = True) M = pymc.MCMC([iw, mnc]) From 7268524199fa518fed9727254abc8cf9da7b695f Mon Sep 17 00:00:00 2001 From: Wes McKinney Date: Thu, 4 Nov 2010 00:59:41 -0400 Subject: [PATCH 11/11] renamed Sigma to C --- pymc/tests/test_distributions.py | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/pymc/tests/test_distributions.py b/pymc/tests/test_distributions.py index 9682720dbf..53f2367d2d 100644 --- a/pymc/tests/test_distributions.py +++ b/pymc/tests/test_distributions.py @@ -993,15 +993,15 @@ class test_inverse_wishart(TestCase): Adapted from test_wishart """ # Covariance matrix (!) - Sigma_test = _Tau_test.I + C_test = _Tau_test.I # def test_samples(self): # # test consistency between precision and cov-based # _npr.seed(1) - # sample_a = rinverse_wishart(100, self.Sigma_test) + # sample_a = rinverse_wishart(100, self.C_test) # _npr.seed(1) - # sample_b = rinverse_wishart_prec(100, self.Sigma_test.I) + # sample_b = rinverse_wishart_prec(100, self.C_test.I) # assert_array_almost_equal(sample_a, sample_b) @@ -1011,7 +1011,7 @@ def test_likelihoods(self): except: raise nose.SkipTest, "SciPy not installed." - IW_test = rinverse_wishart(100,self.Sigma_test) + IW_test = rinverse_wishart(100,self.C_test) def slo_inv_wishart(W,n,V): p = W.shape[0] @@ -1026,11 +1026,11 @@ def slo_inv_wishart(W,n,V): return logp for i in [5,10,100,10000]: - right_answer = slo_inv_wishart(IW_test,i,self.Sigma_test) - calculated = inverse_wishart_like(IW_test,i,self.Sigma_test) + right_answer = slo_inv_wishart(IW_test,i,self.C_test) + calculated = inverse_wishart_like(IW_test,i,self.C_test) assert_array_almost_equal(calculated, right_answer, decimal=1) - calculated = inverse_wishart_prec_like(IW_test,i,self.Sigma_test.I) + calculated = inverse_wishart_prec_like(IW_test,i,self.C_test.I) assert_array_almost_equal(calculated, right_answer, decimal=1) """ @@ -1039,11 +1039,11 @@ def test_expval(self): n = 100 N = 1000 - A = 0.*self.Sigma_test + A = 0.*self.C_test for i in xrange(N): - A += rinverse_wishart(n,self.Sigma_test) + A += rinverse_wishart(n,self.C_test) A /= N - delta=A-inverse_wishart_expval(n,self.Sigma_test) + delta=A-inverse_wishart_expval(n,self.C_test) print np.abs(np.asarray(delta)/np.asarray(A)).max() assert(np.abs(np.asarray(delta)/np.asarray(A)).max()<.5) """