diff --git a/pymc3/distributions/continuous.py b/pymc3/distributions/continuous.py index 76d05f9d20..1be24fa4a4 100644 --- a/pymc3/distributions/continuous.py +++ b/pymc3/distributions/continuous.py @@ -13,7 +13,7 @@ from scipy.interpolate import InterpolatedUnivariateSpline import warnings -from pymc3.theanof import floatX +from pymc3.theanof import floatX, intX from . import transforms from pymc3.util import get_variable_name @@ -45,6 +45,7 @@ def __init__(self, transform=transforms.logodds, *args, **kwargs): super(UnitContinuous, self).__init__( transform=transform, *args, **kwargs) + def assert_negative_support(var, label, distname, value=-1e-6): # Checks for evidence of positive support for a variable if var is None: @@ -136,12 +137,11 @@ def __init__(self, lower=0, upper=1, transform='interval', transform = transforms.interval(lower, upper) super(Uniform, self).__init__(transform=transform, *args, **kwargs) - self.lower = lower = tt.as_tensor_variable(lower) - self.upper = upper = tt.as_tensor_variable(upper) - self.mean = (upper + lower) / 2. + self.lower = floatX(lower) + self.upper = floatX(upper) + self.mean = (self.upper + self.lower) / 2. self.median = self.mean - def random(self, point=None, size=None, repeat=None): lower, upper = draw_values([self.lower, self.upper], point=point) @@ -224,14 +224,14 @@ class Normal(Continuous): def __init__(self, mu=0, sd=None, tau=None, **kwargs): tau, sd = get_tau_sd(tau=tau, sd=sd) - self.sd = tt.as_tensor_variable(sd) - self.tau = tt.as_tensor_variable(tau) + self.sd = sd + self.tau = tau - self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu) + self.mean = self.median = self.mode = self.mu = floatX(mu) self.variance = 1. / self.tau - assert_negative_support(sd, 'sd', 'Normal') - assert_negative_support(tau, 'tau', 'Normal') + assert_negative_support(self.sd, 'sd', 'Normal') + assert_negative_support(self.tau, 'tau', 'Normal') super(Normal, self).__init__(**kwargs) @@ -246,9 +246,9 @@ def logp(self, value): sd = self.sd tau = self.tau mu = self.mu - - return bound((-tau * (value - mu)**2 + tt.log(tau / np.pi / 2.)) / 2., - sd > 0) + return bound( + (-tau * (value - mu)**2 + tt.log(tau / floatX(np.pi * 2.))) / 2., + sd > 0) def _repr_latex_(self, name=None, dist=None): if dist is None: @@ -288,14 +288,14 @@ def __init__(self, sd=None, tau=None, *args, **kwargs): super(HalfNormal, self).__init__(*args, **kwargs) tau, sd = get_tau_sd(tau=tau, sd=sd) - self.sd = sd = tt.as_tensor_variable(sd) - self.tau = tau = tt.as_tensor_variable(tau) + self.sd = floatX(sd) + self.tau = floatX(tau) self.mean = tt.sqrt(2 / (np.pi * self.tau)) self.variance = (1. - 2 / np.pi) / self.tau - assert_negative_support(tau, 'tau', 'HalfNormal') - assert_negative_support(sd, 'sd', 'HalfNormal') + assert_negative_support(self.tau, 'tau', 'HalfNormal') + assert_negative_support(self.sd, 'sd', 'HalfNormal') def random(self, point=None, size=None, repeat=None): sd = draw_values([self.sd], point=point) @@ -317,6 +317,7 @@ def _repr_latex_(self, name=None, dist=None): return r'${} \sim \text{{HalfNormal}}(\mathit{{sd}}={})$'.format(name, get_variable_name(sd)) + class Wald(PositiveContinuous): R""" Wald log-likelihood. @@ -376,19 +377,19 @@ class Wald(PositiveContinuous): def __init__(self, mu=None, lam=None, phi=None, alpha=0., *args, **kwargs): super(Wald, self).__init__(*args, **kwargs) mu, lam, phi = self.get_mu_lam_phi(mu, lam, phi) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.mu = mu = tt.as_tensor_variable(mu) - self.lam = lam = tt.as_tensor_variable(lam) - self.phi = phi = tt.as_tensor_variable(phi) + self.alpha = floatX(alpha) + self.mu = floatX(mu) + self.lam = floatX(lam) + self.phi = floatX(phi) self.mean = self.mu + self.alpha self.mode = self.mu * (tt.sqrt(1. + (1.5 * self.mu / self.lam)**2) - 1.5 * self.mu / self.lam) + self.alpha self.variance = (self.mu**3) / self.lam - assert_negative_support(phi, 'phi', 'Wald') - assert_negative_support(mu, 'mu', 'Wald') - assert_negative_support(lam, 'lam', 'Wald') + assert_negative_support(self.phi, 'phi', 'Wald') + assert_negative_support(self.mu, 'mu', 'Wald') + assert_negative_support(self.lam, 'lam', 'Wald') def get_mu_lam_phi(self, mu, lam, phi): if mu is None: @@ -497,15 +498,15 @@ def __init__(self, alpha=None, beta=None, mu=None, sd=None, super(Beta, self).__init__(*args, **kwargs) alpha, beta = self.get_alpha_beta(alpha, beta, mu, sd) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.beta = beta = tt.as_tensor_variable(beta) + self.alpha = floatX(alpha) + self.beta = floatX(beta) self.mean = self.alpha / (self.alpha + self.beta) self.variance = self.alpha * self.beta / ( (self.alpha + self.beta)**2 * (self.alpha + self.beta + 1)) - assert_negative_support(alpha, 'alpha', 'Beta') - assert_negative_support(beta, 'beta', 'Beta') + assert_negative_support(self.alpha, 'alpha', 'Beta') + assert_negative_support(self.beta, 'beta', 'Beta') def get_alpha_beta(self, alpha=None, beta=None, mu=None, sd=None): if (alpha is not None) and (beta is not None): @@ -568,14 +569,14 @@ class Exponential(PositiveContinuous): def __init__(self, lam, *args, **kwargs): super(Exponential, self).__init__(*args, **kwargs) - self.lam = lam = tt.as_tensor_variable(lam) - self.mean = 1. / self.lam + self.lam = floatX(lam) + self.mean = floatX(1.) / self.lam self.median = self.mean * tt.log(2) - self.mode = tt.zeros_like(self.lam) + self.mode = tt.zeros_like(tt.as_tensor_variable(self.lam)) self.variance = self.lam**-2 - assert_negative_support(lam, 'lam', 'Exponential') + assert_negative_support(self.lam, 'lam', 'Exponential') def random(self, point=None, size=None, repeat=None): lam = draw_values([self.lam], point=point) @@ -594,6 +595,7 @@ def _repr_latex_(self, name=None, dist=None): return r'${} \sim \text{{Exponential}}(\mathit{{lam}}={})$'.format(name, get_variable_name(lam)) + class Laplace(Continuous): R""" Laplace log-likelihood. @@ -619,8 +621,8 @@ class Laplace(Continuous): def __init__(self, mu, b, *args, **kwargs): super(Laplace, self).__init__(*args, **kwargs) - self.b = b = tt.as_tensor_variable(b) - self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu) + self.b = b = floatX(b) + self.mean = self.median = self.mode = self.mu = mu = floatX(mu) self.variance = 2 * self.b**2 @@ -681,9 +683,9 @@ def __init__(self, mu=0, sd=None, tau=None, *args, **kwargs): super(Lognormal, self).__init__(*args, **kwargs) tau, sd = get_tau_sd(tau=tau, sd=sd) - self.mu = mu = tt.as_tensor_variable(mu) - self.tau = tau = tt.as_tensor_variable(tau) - self.sd = sd = tt.as_tensor_variable(sd) + self.mu = mu = floatX(mu) + self.tau = tau = floatX(tau) + self.sd = sd = floatX(sd) self.mean = tt.exp(self.mu + 1. / (2 * self.tau)) self.median = tt.exp(self.mu) @@ -752,14 +754,14 @@ class StudentT(Continuous): def __init__(self, nu, mu=0, lam=None, sd=None, *args, **kwargs): super(StudentT, self).__init__(*args, **kwargs) - self.nu = nu = tt.as_tensor_variable(nu) + self.nu = floatX(nu) lam, sd = get_tau_sd(tau=lam, sd=sd) - self.lam = lam = tt.as_tensor_variable(lam) - self.sd = sd = tt.as_tensor_variable(sd) - self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu) + self.lam = floatX(lam) + self.sd = floatX(sd) + self.mean = self.median = self.mode = self.mu = mu = floatX(mu) - self.variance = tt.switch((nu > 2) * 1, - (1 / self.lam) * (nu / (nu - 2)), + self.variance = tt.switch((self.nu > 2) * 1, + (1 / self.lam) * (self.nu / (self.nu - 2)), np.inf) assert_negative_support(lam, 'lam (sd)', 'StudentT') @@ -824,20 +826,19 @@ class Pareto(PositiveContinuous): def __init__(self, alpha, m, *args, **kwargs): super(Pareto, self).__init__(*args, **kwargs) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.m = m = tt.as_tensor_variable(m) + self.alpha = floatX(alpha) + self.m = floatX(m) - self.mean = tt.switch(tt.gt(alpha, 1), alpha * - m / (alpha - 1.), np.inf) - self.median = m * 2.**(1. / alpha) + self.mean = tt.switch(tt.gt(self.alpha, 1), self.alpha * + self.m / (self.alpha - 1.), np.inf) + self.median = self.m * 2.**(1. / self.alpha) self.variance = tt.switch( - tt.gt(alpha, 2), - (alpha * m**2) / ((alpha - 2.) * (alpha - 1.)**2), + tt.gt(self.alpha, 2), + (self.alpha * self.m**2) / ((self.alpha - 2.) * (self.alpha - 1.)**2), np.inf) - assert_negative_support(alpha, 'alpha', 'Pareto') - assert_negative_support(m, 'm', 'Pareto') - + assert_negative_support(self.alpha, 'alpha', 'Pareto') + assert_negative_support(self.m, 'm', 'Pareto') def _random(self, alpha, m, size=None): u = np.random.uniform(size=size) @@ -895,10 +896,10 @@ class Cauchy(Continuous): def __init__(self, alpha, beta, *args, **kwargs): super(Cauchy, self).__init__(*args, **kwargs) - self.median = self.mode = self.alpha = tt.as_tensor_variable(alpha) - self.beta = tt.as_tensor_variable(beta) + self.median = self.mode = self.alpha = floatX(alpha) + self.beta = floatX(beta) - assert_negative_support(beta, 'beta', 'Cauchy') + assert_negative_support(self.beta, 'beta', 'Cauchy') def _random(self, alpha, beta, size=None): u = np.random.uniform(size=size) @@ -951,11 +952,11 @@ class HalfCauchy(PositiveContinuous): def __init__(self, beta, *args, **kwargs): super(HalfCauchy, self).__init__(*args, **kwargs) - self.mode = tt.as_tensor_variable(0) - self.median = tt.as_tensor_variable(beta) - self.beta = tt.as_tensor_variable(beta) + self.mode = floatX(0) + self.median = floatX(beta) + self.beta = floatX(beta) - assert_negative_support(beta, 'beta', 'HalfCauchy') + assert_negative_support(self.beta, 'beta', 'HalfCauchy') def _random(self, beta, size=None): u = np.random.uniform(size=size) @@ -980,6 +981,7 @@ def _repr_latex_(self, name=None, dist=None): return r'${} \sim \text{{HalfCauchy}}(\mathit{{beta}}={})$'.format(name, get_variable_name(beta)) + class Gamma(PositiveContinuous): R""" Gamma log-likelihood. @@ -1023,14 +1025,14 @@ def __init__(self, alpha=None, beta=None, mu=None, sd=None, *args, **kwargs): super(Gamma, self).__init__(*args, **kwargs) alpha, beta = self.get_alpha_beta(alpha, beta, mu, sd) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.beta = beta = tt.as_tensor_variable(beta) - self.mean = alpha / beta - self.mode = tt.maximum((alpha - 1) / beta, 0) - self.variance = alpha / beta**2 + self.alpha = floatX(alpha) + self.beta = floatX(beta) + self.mean = self.alpha / self.beta + self.mode = tt.maximum((self.alpha - 1) / self.beta, 0) + self.variance = self.alpha / self.beta**2 - assert_negative_support(alpha, 'alpha', 'Gamma') - assert_negative_support(beta, 'beta', 'Gamma') + assert_negative_support(self.alpha, 'alpha', 'Gamma') + assert_negative_support(self.beta, 'beta', 'Gamma') def get_alpha_beta(self, alpha=None, beta=None, mu=None, sd=None): if (alpha is not None) and (beta is not None): @@ -1099,16 +1101,16 @@ class InverseGamma(PositiveContinuous): def __init__(self, alpha, beta=1, *args, **kwargs): super(InverseGamma, self).__init__(*args, **kwargs) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.beta = beta = tt.as_tensor_variable(beta) + self.alpha = floatX(alpha) + self.beta = floatX(beta) self.mean = self._calculate_mean() - self.mode = beta / (alpha + 1.) - self.variance = tt.switch(tt.gt(alpha, 2), - (beta**2) / (alpha * (alpha - 1.)**2), + self.mode = self.beta / (self.alpha + 1.) + self.variance = tt.switch(tt.gt(self.alpha, 2), + (self.beta**2) / (self.alpha * (self.alpha - 1.)**2), np.inf) - assert_negative_support(alpha, 'alpha', 'InverseGamma') - assert_negative_support(beta, 'beta', 'InverseGamma') + assert_negative_support(self.alpha, 'alpha', 'InverseGamma') + assert_negative_support(self.beta, 'beta', 'InverseGamma') def _calculate_mean(self): m = self.beta / (self.alpha - 1.) @@ -1163,7 +1165,7 @@ class ChiSquared(Gamma): """ def __init__(self, nu, *args, **kwargs): - self.nu = nu = tt.as_tensor_variable(nu) + self.nu = nu = intX(nu) super(ChiSquared, self).__init__(alpha=nu / 2., beta=0.5, *args, **kwargs) @@ -1201,15 +1203,15 @@ class Weibull(PositiveContinuous): def __init__(self, alpha, beta, *args, **kwargs): super(Weibull, self).__init__(*args, **kwargs) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.beta = beta = tt.as_tensor_variable(beta) - self.mean = beta * tt.exp(gammaln(1 + 1. / alpha)) - self.median = beta * tt.exp(gammaln(tt.log(2)))**(1. / alpha) - self.variance = (beta**2) * \ - tt.exp(gammaln(1 + 2. / alpha - self.mean**2)) + self.alpha = floatX(alpha) + self.beta = floatX(beta) + self.mean = self.beta * tt.exp(gammaln(1 + 1. / self.alpha)) + self.median = self.beta * tt.exp(gammaln(tt.log(2)))**(1. / self.alpha) + self.variance = (self.beta**2) * \ + tt.exp(gammaln(1 + 2. / self.alpha - self.mean**2)) - assert_negative_support(alpha, 'alpha', 'Weibull') - assert_negative_support(beta, 'beta', 'Weibull') + assert_negative_support(self.alpha, 'alpha', 'Weibull') + assert_negative_support(self.beta, 'beta', 'Weibull') def random(self, point=None, size=None, repeat=None): alpha, beta = draw_values([self.alpha, self.beta], @@ -1295,14 +1297,14 @@ class ExGaussian(Continuous): def __init__(self, mu, sigma, nu, *args, **kwargs): super(ExGaussian, self).__init__(*args, **kwargs) - self.mu = mu = tt.as_tensor_variable(mu) - self.sigma = sigma = tt.as_tensor_variable(sigma) - self.nu = nu = tt.as_tensor_variable(nu) - self.mean = mu + nu - self.variance = (sigma**2) + (nu**2) + self.mu = floatX(mu) + self.sigma = floatX(sigma) + self.nu = floatX(nu) + self.mean = self.mu + self.nu + self.variance = (self.sigma**2) + (self.nu**2) - assert_negative_support(sigma, 'sigma', 'ExGaussian') - assert_negative_support(nu, 'nu', 'ExGaussian') + assert_negative_support(self.sigma, 'sigma', 'ExGaussian') + assert_negative_support(self.nu, 'nu', 'ExGaussian') def random(self, point=None, size=None, repeat=None): mu, sigma, nu = draw_values([self.mu, self.sigma, self.nu], @@ -1370,9 +1372,9 @@ def __init__(self, mu=0.0, kappa=None, transform='circular', if transform == 'circular': transform = transforms.Circular() super(VonMises, self).__init__(transform=transform, *args, **kwargs) - self.mean = self.median = self.mode = self.mu = mu = tt.as_tensor_variable(mu) - self.kappa = kappa = tt.as_tensor_variable(kappa) - self.variance = 1 - i1(kappa) / i0(kappa) + self.mean = self.median = self.mode = self.mu = floatX(mu) + self.kappa = floatX(kappa) + self.variance = 1 - i1(self.kappa) / i0(self.kappa) assert_negative_support(kappa, 'kappa', 'VonMises') @@ -1398,7 +1400,6 @@ def _repr_latex_(self, name=None, dist=None): get_variable_name(kappa)) - class SkewNormal(Continuous): R""" Univariate skew-normal log-likelihood. @@ -1441,17 +1442,17 @@ class SkewNormal(Continuous): def __init__(self, mu=0.0, sd=None, tau=None, alpha=1, *args, **kwargs): super(SkewNormal, self).__init__(*args, **kwargs) tau, sd = get_tau_sd(tau=tau, sd=sd) - self.mu = mu = tt.as_tensor_variable(mu) - self.tau = tt.as_tensor_variable(tau) - self.sd = tt.as_tensor_variable(sd) + self.mu = floatX(mu) + self.tau = floatX(tau) + self.sd = floatX(sd) - self.alpha = alpha = tt.as_tensor_variable(alpha) + self.alpha = floatX(alpha) - self.mean = mu + self.sd * (2 / np.pi)**0.5 * alpha / (1 + alpha**2)**0.5 - self.variance = self.sd**2 * (1 - (2 * alpha**2) / ((1 + alpha**2) * np.pi)) + self.mean = mu + self.sd * (2 / np.pi)**0.5 * self.alpha / (1 + self.alpha**2)**0.5 + self.variance = self.sd**2 * (1 - (2 * self.alpha**2) / ((1 + self.alpha**2) * np.pi)) - assert_negative_support(tau, 'tau', 'SkewNormal') - assert_negative_support(sd, 'sd', 'SkewNormal') + assert_negative_support(self.tau, 'tau', 'SkewNormal') + assert_negative_support(self.sd, 'sd', 'SkewNormal') def random(self, point=None, size=None, repeat=None): mu, tau, _, alpha = draw_values( @@ -1504,9 +1505,9 @@ def __init__(self, lower=0, upper=1, c=0.5, *args, **kwargs): super(Triangular, self).__init__(*args, **kwargs) - self.median = self.mean = self.c = c = tt.as_tensor_variable(c) - self.lower = lower = tt.as_tensor_variable(lower) - self.upper = upper = tt.as_tensor_variable(upper) + self.median = self.mean = self.c = c = floatX(c) + self.lower = floatX(lower) + self.upper = floatX(upper) def random(self, point=None, size=None): c, lower, upper = draw_values([self.c, self.lower, self.upper], @@ -1557,10 +1558,10 @@ class Gumbel(Continuous): """ def __init__(self, mu=0, beta=1.0, **kwargs): - self.mu = tt.as_tensor_variable(mu) - self.beta = tt.as_tensor_variable(beta) + self.mu = floatX(mu) + self.beta = floatX(beta) - assert_negative_support(beta, 'beta', 'Gumbel') + assert_negative_support(self.beta, 'beta', 'Gumbel') self.mean = self.mu + self.beta * np.euler_gamma self.median = self.mu - self.beta * tt.log(tt.log(2)) @@ -1627,7 +1628,7 @@ def __init__(self, x_points, pdf_points, transform='interval', interp = InterpolatedUnivariateSpline(x_points, pdf_points, k=1, ext='zeros') Z = interp.integral(x_points[0], x_points[-1]) - self.Z = tt.as_tensor_variable(Z) + self.Z = floatX(Z) self.interp_op = SplineWrapper(interp) self.x_points = x_points self.pdf_points = pdf_points / Z diff --git a/pymc3/distributions/discrete.py b/pymc3/distributions/discrete.py index a15168b07a..dc6aba9f85 100644 --- a/pymc3/distributions/discrete.py +++ b/pymc3/distributions/discrete.py @@ -4,6 +4,7 @@ import theano.tensor as tt from scipy import stats +from pymc3.theanof import floatX, intX from pymc3.util import get_variable_name from .dist_math import bound, factln, binomln, betaln, logpow from .distribution import Discrete, draw_values, generate_samples, reshape_sampled @@ -42,9 +43,9 @@ class Binomial(Discrete): def __init__(self, n, p, *args, **kwargs): super(Binomial, self).__init__(*args, **kwargs) - self.n = n = tt.as_tensor_variable(n) - self.p = p = tt.as_tensor_variable(p) - self.mode = tt.cast(tround(n * p), self.dtype) + self.n = intX(n) + self.p = floatX(p) + self.mode = tt.cast(tround(self.n * self.p), self.dtype) def random(self, point=None, size=None, repeat=None): n, p = draw_values([self.n, self.p], point=point) @@ -101,10 +102,11 @@ class BetaBinomial(Discrete): def __init__(self, alpha, beta, n, *args, **kwargs): super(BetaBinomial, self).__init__(*args, **kwargs) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.beta = beta = tt.as_tensor_variable(beta) - self.n = n = tt.as_tensor_variable(n) - self.mode = tt.cast(tround(alpha / (alpha + beta)), 'int8') + self.alpha = floatX(alpha) + self.beta = floatX(beta) + self.n = intX(n) + self.mode = tt.cast( + tround(self.alpha / (self.alpha + self.beta)), 'int8') def _random(self, alpha, beta, n, size=None): size = size or 1 @@ -166,8 +168,8 @@ class Bernoulli(Discrete): def __init__(self, p, *args, **kwargs): super(Bernoulli, self).__init__(*args, **kwargs) - self.p = p = tt.as_tensor_variable(p) - self.mode = tt.cast(tround(p), 'int8') + self.p = floatX(p) + self.mode = tt.cast(tround(self.p), 'int8') def random(self, point=None, size=None, repeat=None): p = draw_values([self.p], point=point) @@ -207,8 +209,8 @@ class DiscreteWeibull(Discrete): def __init__(self, q, beta, *args, **kwargs): super(DiscreteWeibull, self).__init__(*args, defaults=['median'], **kwargs) - self.q = q = tt.as_tensor_variable(q) - self.beta = beta = tt.as_tensor_variable(beta) + self.q = floatX(q) + self.beta = floatX(beta) self.median = self._ppf(0.5) @@ -282,8 +284,8 @@ class Poisson(Discrete): def __init__(self, mu, *args, **kwargs): super(Poisson, self).__init__(*args, **kwargs) - self.mu = mu = tt.as_tensor_variable(mu) - self.mode = tt.floor(mu).astype('int32') + self.mu = floatX(mu) + self.mode = tt.floor(intX(self.mu)) def random(self, point=None, size=None, repeat=None): mu = draw_values([self.mu], point=point) @@ -336,9 +338,9 @@ class NegativeBinomial(Discrete): def __init__(self, mu, alpha, *args, **kwargs): super(NegativeBinomial, self).__init__(*args, **kwargs) - self.mu = mu = tt.as_tensor_variable(mu) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.mode = tt.floor(mu).astype('int32') + self.mu = floatX(mu) + self.alpha = floatX(alpha) + self.mode = tt.floor(intX(self.mu)) def random(self, point=None, size=None, repeat=None): mu, alpha = draw_values([self.mu, self.alpha], point=point) @@ -394,7 +396,7 @@ class Geometric(Discrete): def __init__(self, p, *args, **kwargs): super(Geometric, self).__init__(*args, **kwargs) - self.p = p = tt.as_tensor_variable(p) + self.p = floatX(p) self.mode = 1 def random(self, point=None, size=None, repeat=None): @@ -438,17 +440,17 @@ class DiscreteUniform(Discrete): def __init__(self, lower, upper, *args, **kwargs): super(DiscreteUniform, self).__init__(*args, **kwargs) - self.lower = tt.floor(lower).astype('int32') - self.upper = tt.floor(upper).astype('int32') + self.lower = intX(np.floor(lower)) + self.upper = intX(np.floor(upper)) self.mode = tt.maximum( - tt.floor((upper + lower) / 2.).astype('int32'), self.lower) + intX(np.floor((upper + lower) / 2.)), self.lower) def _random(self, lower, upper, size=None): # This way seems to be the only to deal with lower and upper # as array-like. samples = stats.uniform.rvs(lower, upper - lower - np.finfo(float).eps, size=size) - return np.floor(samples).astype('int32') + return intX(np.floor(samples)) def random(self, point=None, size=None, repeat=None): lower, upper = draw_values([self.lower, self.upper], point=point) @@ -498,9 +500,9 @@ def __init__(self, p, *args, **kwargs): self.k = tt.shape(p)[-1].tag.test_value except AttributeError: self.k = tt.shape(p)[-1] - self.p = p = tt.as_tensor_variable(p) + p = floatX(p) self.p = (p.T / tt.sum(p, -1)).T - self.mode = tt.argmax(p) + self.mode = tt.argmax(self.p) def random(self, point=None, size=None, repeat=None): def random_choice(k, *args, **kwargs): @@ -556,7 +558,7 @@ class Constant(Discrete): def __init__(self, c, *args, **kwargs): super(Constant, self).__init__(*args, **kwargs) - self.mean = self.median = self.mode = self.c = c = tt.as_tensor_variable(c) + self.mean = self.median = self.mode = self.c = floatX(c) def random(self, point=None, size=None, repeat=None): c = draw_values([self.c], point=point) @@ -616,8 +618,8 @@ class ZeroInflatedPoisson(Discrete): def __init__(self, psi, theta, *args, **kwargs): super(ZeroInflatedPoisson, self).__init__(*args, **kwargs) - self.theta = theta = tt.as_tensor_variable(theta) - self.psi = psi = tt.as_tensor_variable(psi) + self.theta = floatX(theta) + self.psi = floatX(psi) self.pois = Poisson.dist(theta) self.mode = self.pois.mode @@ -682,9 +684,9 @@ class ZeroInflatedBinomial(Discrete): def __init__(self, psi, n, p, *args, **kwargs): super(ZeroInflatedBinomial, self).__init__(*args, **kwargs) - self.n = n = tt.as_tensor_variable(n) - self.p = p = tt.as_tensor_variable(p) - self.psi = psi = tt.as_tensor_variable(psi) + self.n = n = intX(n) + self.p = p = floatX(p) + self.psi = psi = floatX(psi) self.bin = Binomial.dist(n, p) self.mode = self.bin.mode @@ -756,10 +758,10 @@ class ZeroInflatedNegativeBinomial(Discrete): def __init__(self, psi, mu, alpha, *args, **kwargs): super(ZeroInflatedNegativeBinomial, self).__init__(*args, **kwargs) - self.mu = mu = tt.as_tensor_variable(mu) - self.alpha = alpha = tt.as_tensor_variable(alpha) - self.psi = psi = tt.as_tensor_variable(psi) - self.nb = NegativeBinomial.dist(mu, alpha) + self.mu = floatX(mu) + self.alpha = floatX(alpha) + self.psi = floatX(psi) + self.nb = NegativeBinomial.dist(self.mu, self.alpha) self.mode = self.nb.mode def random(self, point=None, size=None, repeat=None): diff --git a/pymc3/distributions/dist_math.py b/pymc3/distributions/dist_math.py index 42a0603856..816b904687 100644 --- a/pymc3/distributions/dist_math.py +++ b/pymc3/distributions/dist_math.py @@ -12,7 +12,7 @@ from .special import gammaln from ..math import logdet as _logdet -from pymc3.theanof import floatX +from pymc3.theanof import floatX, intX f = floatX c = - .5 * np.log(2. * np.pi) @@ -44,7 +44,7 @@ def bound(logp, *conditions, **kwargs): else: alltrue = alltrue_scalar - return tt.switch(alltrue(conditions), logp, -np.inf) + return tt.switch(alltrue(conditions), logp, -f(np.inf)) def alltrue_elemwise(vals): @@ -62,12 +62,13 @@ def logpow(x, m): """ Calculates log(x**m) since m*log(x) will fail when m, x = 0. """ + m = f(m) # return m * log(x) return tt.switch(tt.eq(x, 0), -np.inf, m * tt.log(x)) def factln(n): - return gammaln(n + 1) + return gammaln(n + intX(1)) def binomln(n, k): @@ -256,7 +257,7 @@ def MvNormalLogp(): result += f(2) * n * tt.sum(tt.log(diag)) result += (delta_trans ** f(2)).sum() result = f(-.5) * result - logp = tt.switch(ok, result, -np.inf) + logp = tt.switch(ok, result, -f(np.inf)) def dlogp(inputs, gradients): g_logp, = gradients diff --git a/pymc3/distributions/mixture.py b/pymc3/distributions/mixture.py index 19953bca54..2712b0e312 100644 --- a/pymc3/distributions/mixture.py +++ b/pymc3/distributions/mixture.py @@ -1,6 +1,7 @@ import numpy as np import theano.tensor as tt +from pymc3.theanof import floatX from pymc3.util import get_variable_name from ..math import logsumexp from .dist_math import bound @@ -42,7 +43,7 @@ class Mixture(Distribution): def __init__(self, w, comp_dists, *args, **kwargs): shape = kwargs.pop('shape', ()) - self.w = w = tt.as_tensor_variable(w) + self.w = floatX(w) self.comp_dists = comp_dists defaults = kwargs.pop('defaults', []) @@ -53,7 +54,7 @@ def __init__(self, w, comp_dists, *args, **kwargs): dtype = kwargs.pop('dtype', 'float64') try: - self.mean = (w * self._comp_means()).sum(axis=-1) + self.mean = (self.w * self._comp_means()).sum(axis=-1) if 'mean' not in defaults: defaults.append('mean') @@ -63,7 +64,7 @@ def __init__(self, w, comp_dists, *args, **kwargs): try: comp_modes = self._comp_modes() comp_mode_logps = self.logp(comp_modes) - self.mode = comp_modes[tt.argmax(w * comp_mode_logps, axis=-1)] + self.mode = comp_modes[tt.argmax(self.w * comp_mode_logps, axis=-1)] if 'mode' not in defaults: defaults.append('mode') @@ -169,9 +170,9 @@ class NormalMixture(Mixture): def __init__(self, w, mu, *args, **kwargs): _, sd = get_tau_sd(tau=kwargs.pop('tau', None), sd=kwargs.pop('sd', None)) - self.mu = mu = tt.as_tensor_variable(mu) - self.sd = sd = tt.as_tensor_variable(sd) - super(NormalMixture, self).__init__(w, Normal.dist(mu, sd=sd), + self.mu = floatX(mu) + self.sd = floatX(sd) + super(NormalMixture, self).__init__(w, Normal.dist(self.mu, sd=self.sd), *args, **kwargs) def _repr_latex_(self, name=None, dist=None): diff --git a/pymc3/distributions/multivariate.py b/pymc3/distributions/multivariate.py index 25ab38dc81..aec6a9ed3d 100755 --- a/pymc3/distributions/multivariate.py +++ b/pymc3/distributions/multivariate.py @@ -14,7 +14,7 @@ import pymc3 as pm from pymc3.math import tround -from pymc3.theanof import floatX +from pymc3.theanof import floatX, intX from . import transforms from pymc3.util import get_variable_name from .distribution import Continuous, Discrete, draw_values, generate_samples @@ -105,8 +105,7 @@ def __init__(self, mu, cov=None, tau=None, chol=None, lower=True, raise ValueError('Incompatible parameterization. ' 'Specify exactly one of tau, cov, ' 'or chol.') - mu = tt.as_tensor_variable(mu) - self.mean = self.median = self.mode = self.mu = mu + self.mean = self.median = self.mode = self.mu = floatX(mu) self.solve_lower = tt.slinalg.Solve(A_structure="lower_triangular") # Step methods and advi do not catch LinAlgErrors at the # moment. We work around that by using a cholesky op @@ -117,7 +116,7 @@ def __init__(self, mu, cov=None, tau=None, chol=None, lower=True, if cov is not None: self.k = cov.shape[0] self._cov_type = 'cov' - cov = tt.as_tensor_variable(cov) + cov = floatX(cov) if cov.ndim != 2: raise ValueError('cov must be two dimensional.') self.chol_cov = cholesky(cov) @@ -125,7 +124,7 @@ def __init__(self, mu, cov=None, tau=None, chol=None, lower=True, elif tau is not None: self.k = tau.shape[0] self._cov_type = 'tau' - tau = tt.as_tensor_variable(tau) + tau = floatX(tau) if tau.ndim != 2: raise ValueError('tau must be two dimensional.') self.chol_tau = cholesky(tau) @@ -135,7 +134,7 @@ def __init__(self, mu, cov=None, tau=None, chol=None, lower=True, self._cov_type = 'chol' if chol.ndim != 2: raise ValueError('chol must be two dimensional.') - self.chol_cov = tt.as_tensor_variable(chol) + self.chol_cov = floatX(chol) def random(self, point=None, size=None): if size is None: @@ -247,7 +246,7 @@ def _repr_latex_(self, name=None, dist=None): mu = dist.mu try: cov = dist.cov - except AttributeErrir: + except AttributeError: cov = dist.chol_cov return r'${} \sim \text{{MvNormal}}(\mathit{{mu}}={}, \mathit{{cov}}={})$'.format(name, get_variable_name(mu), @@ -291,9 +290,9 @@ class MvStudentT(Continuous): def __init__(self, nu, Sigma, mu=None, *args, **kwargs): super(MvStudentT, self).__init__(*args, **kwargs) - self.nu = nu = tt.as_tensor_variable(nu) - mu = tt.zeros(Sigma.shape[0]) if mu is None else tt.as_tensor_variable(mu) - self.Sigma = Sigma = tt.as_tensor_variable(Sigma) + self.nu = floatX(nu) + mu = tt.zeros(Sigma.shape[0]) if mu is None else floatX(mu) + self.Sigma = floatX(Sigma) self.mean = self.median = self.mode = self.mu = mu @@ -372,12 +371,12 @@ def __init__(self, a, transform=transforms.stick_breaking, kwargs.setdefault("shape", shape) super(Dirichlet, self).__init__(transform=transform, *args, **kwargs) - self.k = tt.as_tensor_variable(shape) - self.a = a = tt.as_tensor_variable(a) + self.k = floatX(shape) + self.a = floatX(a) self.mean = a / tt.sum(a) - self.mode = tt.switch(tt.all(a > 1), - (a - 1) / tt.sum(a - 1), + self.mode = tt.switch(tt.all(self.a > 1), + (a - 1) / tt.sum(self.a - 1), np.nan) def random(self, point=None, size=None): @@ -457,8 +456,8 @@ def __init__(self, n, p, *args, **kwargs): self.n = tt.shape_padright(n) self.p = p if p.ndim == 2 else tt.shape_padleft(p) else: - self.n = tt.as_tensor_variable(n) - self.p = tt.as_tensor_variable(p) + self.n = intX(n) + self.p = floatX(p) self.mean = self.n * self.p self.mode = tt.cast(tround(self.mean), 'int32') @@ -594,12 +593,12 @@ def __init__(self, nu, V, *args, **kwargs): 'on the issues surrounding the Wishart see here: ' 'https://github.com/pymc-devs/pymc3/issues/538.', UserWarning) - self.nu = nu = tt.as_tensor_variable(nu) - self.p = p = tt.as_tensor_variable(V.shape[0]) - self.V = V = tt.as_tensor_variable(V) + self.nu = intX(nu) + self.p = floatX(V.shape[0]) + self.V = floatX(V) self.mean = nu * V - self.mode = tt.switch(1 * (nu >= p + 1), - (nu - p - 1) * V, + self.mode = tt.switch(1 * (self.nu >= self.p + 1), + (self.nu - self.p - 1) * self.V, np.nan) def logp(self, X): @@ -840,7 +839,7 @@ class LKJCholeskyCov(Continuous): http://math.stackexchange.com/q/130026 """ def __init__(self, eta, n, sd_dist, *args, **kwargs): - self.n = n + self.n = intX(n) self.eta = eta if 'transform' in kwargs: @@ -943,13 +942,13 @@ def __init__(self, eta=None, n=None, p=None, transform='interval', *args, **kwar 'dimension parameter p -> n. Please update your code. ' 'Automatically re-assigning parameters for backwards compatibility.', DeprecationWarning) - self.n = p - self.eta = n + self.n = floatX(p) + self.eta = intX(n) eta = self.eta n = self.n elif (n is not None) and (eta is not None) and (p is None): - self.n = n - self.eta = eta + self.n = intX(n) + self.eta = floatX(eta) else: raise ValueError('Invalid parameter: please use eta as the shape parameter and ' 'n as the dimension parameter.') diff --git a/pymc3/tests/test_distributions.py b/pymc3/tests/test_distributions.py index ac0eedbddc..42c4cf6a7f 100644 --- a/pymc3/tests/test_distributions.py +++ b/pymc3/tests/test_distributions.py @@ -855,10 +855,11 @@ def test_repr_latex_(): x2 = GaussianRandomWalk('Timeseries', mu=x1, sd=1., shape=2) x3 = MvStudentT('Multivariate', nu=5, mu=x2, Sigma=np.diag(np.ones(2)), shape=2) x4 = NormalMixture('Mixture', w=np.array([.5, .5]), mu=x3, sd=x0) + assert x0._repr_latex_()=='$Discrete \\sim \\text{Binomial}(\\mathit{n}=10, \\mathit{p}=0.5)$' assert x1._repr_latex_()=='$Continuous \\sim \\text{Normal}(\\mathit{mu}=0.0, \\mathit{sd}=1.0)$' assert x2._repr_latex_()=='$Timeseries \\sim \\text{GaussianRandomWalk}(\\mathit{mu}=Continuous, \\mathit{sd}=1.0)$' - assert x3._repr_latex_()=='$Multivariate \\sim \\text{MvStudentT}(\\mathit{nu}=5, \\mathit{mu}=Timeseries, \\mathit{Sigma}=array)$' + assert x3._repr_latex_()=='$Multivariate \\sim \\text{MvStudentT}(\\mathit{nu}=5.0, \\mathit{mu}=Timeseries, \\mathit{Sigma}=array)$' assert x4._repr_latex_()=='$Mixture \\sim \\text{NormalMixture}(\\mathit{w}=array, \\mathit{mu}=Multivariate, \\mathit{sigma}=f(Discrete))$' diff --git a/pymc3/theanof.py b/pymc3/theanof.py index b8b5cfa75b..1092852f30 100644 --- a/pymc3/theanof.py +++ b/pymc3/theanof.py @@ -69,6 +69,23 @@ def floatX(X): return np.asarray(X, dtype=theano.config.floatX) +def intX(X): + """ + Convert a theano tensor or numpy array to int16 or int64 + depending on theano.config.floatX type. + """ + if theano.config.floatX == 'float32': + inttype = 'int16' + elif theano.config.floatX == 'float64': + inttype = 'int64' + + try: + return X.astype(inttype) + except AttributeError: + # Scalar passed + return np.asarray(X, dtype=inttype) + + def smartfloatX(x): """ Convert non int types to floatX diff --git a/pymc3/util.py b/pymc3/util.py index 0cb7e875e1..9ad7461807 100644 --- a/pymc3/util.py +++ b/pymc3/util.py @@ -80,7 +80,11 @@ def get_variable_name(variable): """Returns the variable data type if it is a constant, otherwise returns the argument name. """ - name = variable.name + try: + name = variable.name + except AttributeError: + name =None + if name is None: if hasattr(variable, 'get_parents'): try: @@ -88,7 +92,11 @@ def get_variable_name(variable): return 'f(%s)' % ','.join([n for n in names if isinstance(n, str)]) except IndexError: pass - value = variable.eval() + try: + value = variable.eval() + except AttributeError: + value = variable + if not value.shape: return asscalar(value) return 'array'