Merge pull request #730 from pymc-devs/transform

automatic transform application

Author: jsalvatier

pymc3/backends/base.py

@@ -31,6 +31,7 @@ def __init__(self, name, model=None, vars=None):
         self.varnames = [str(var) for var in vars]
         self.fn = model.fastfn(vars)
 
+
         ## Get variable shapes. Most backends will need this
         ## information.
         var_values = list(zip(self.varnames, self.fn(model.test_point)))

pymc3/blocking.py

@@ -47,7 +47,7 @@ def map(self, dpt):
         """
         apt = np.empty(self.ordering.dimensions)
         for var, slc, _ in self.ordering.vmap:
-                apt[slc] = np.ravel(dpt[var])
+                apt[slc] = dpt[var].ravel()
         return apt
 
     def rmap(self, apt):
@@ -61,7 +61,7 @@ def rmap(self, apt):
         dpt = self.dpt.copy()
 
         for var, slc, shp in self.ordering.vmap:
-            dpt[var] = np.reshape(np.atleast_1d(apt)[slc], shp)
+            dpt[var] = apt[slc].reshape(shp)
 
         return dpt
 

pymc3/distributions/continuous.py

@@ -9,12 +9,23 @@
 
 from .dist_math import *
 from numpy.random import uniform as runiform, normal as rnormal
+from .transforms import logtransform, logoddstransform, interval_transform
 
 __all__ = ['Uniform', 'Flat', 'Normal', 'Beta', 'Exponential', 'Laplace',
            'T', 'StudentT', 'Cauchy', 'HalfCauchy', 'Gamma', 'Weibull','Bound',
            'Tpos', 'Lognormal', 'ChiSquared', 'HalfNormal', 'Wald',
            'Pareto', 'InverseGamma']
 
+class PositiveContinuous(Continuous):
+    """Base class for positive continuous distributions"""
+    def __init__(self, transform=logtransform, *args, **kwargs):
+        super(PositiveContinuous, self).__init__(transform=transform, *args, **kwargs)
+
+class UnitContinuous(Continuous):
+    """Base class for continuous distributions on [0,1]"""
+    def __init__(self, transform=logoddstransform, *args, **kwargs):
+        super(UnitContinuous, self).__init__(transform=transform, *args, **kwargs)
+
 def get_tau_sd(tau=None, sd=None):
     """
     Find precision and standard deviation
@@ -70,13 +81,17 @@ class Uniform(Continuous):
     upper : float
         Upper limit (defaults to 1)
     """
-    def __init__(self, lower=0, upper=1, *args, **kwargs):
+    def __init__(self, lower=0, upper=1, transform='interval', *args, **kwargs):
         super(Uniform, self).__init__(*args, **kwargs)
+        
         self.lower = lower
         self.upper = upper
         self.mean = (upper + lower) / 2.
         self.median = self.mean
 
+        if transform is 'interval':
+            self.transform = interval_transform(lower, upper)
+
     def logp(self, value):
         lower = self.lower
         upper = self.upper
@@ -142,7 +157,7 @@ def logp(self, value):
         )
 
 
-class HalfNormal(Continuous):
+class HalfNormal(PositiveContinuous):
     """
     Half-normal log-likelihood, a normal distribution with mean 0 limited
     to the domain :math:`x \in [0, \infty)`.
@@ -206,7 +221,7 @@ def logp(self, value):
 
 
 
-class Beta(Continuous):
+class Beta(UnitContinuous):
     """
     Beta log-likelihood. The conjugate prior for the parameter
     :math:`p` of the binomial distribution.
@@ -268,7 +283,7 @@ def logp(self, value):
             beta > 0)
 
 
-class Exponential(Continuous):
+class Exponential(PositiveContinuous):
     """
     Exponential distribution
 
@@ -320,7 +335,7 @@ def logp(self, value):
         return -log(2 * b) - abs(value - mu) / b
 
 
-class Lognormal(Continuous):
+class Lognormal(PositiveContinuous):
     """
     Log-normal log-likelihood.
 
@@ -409,7 +424,7 @@ def logp(self, value):
 StudentT = T
 
 
-class Pareto(Continuous):
+class Pareto(PositiveContinuous):
     """
     Pareto log-likelihood. The Pareto is a continuous, positive
     probability distribution with two parameters. It is often used
@@ -482,7 +497,7 @@ def logp(self, value):
                                             value - alpha) / beta) ** 2),
             beta > 0)
 
-class HalfCauchy(Continuous):
+class HalfCauchy(PositiveContinuous):
     """
     Half-Cauchy log-likelihood. Simply the absolute value of Cauchy.
 
@@ -510,7 +525,7 @@ def logp(self, value):
             value >= 0)
 
 
-class Gamma(Continuous):
+class Gamma(PositiveContinuous):
     """
     Gamma log-likelihood.
 
@@ -575,7 +590,7 @@ def logp(self, value):
             beta > 0)
 
 
-class InverseGamma(Continuous):
+class InverseGamma(PositiveContinuous):
     """
     Inverse gamma log-likelihood, the reciprocal of the gamma distribution.
 
@@ -634,7 +649,7 @@ def __init__(self, nu, *args, **kwargs):
         super(ChiSquared, self).__init__(alpha=nu/2., beta=0.5, *args, **kwargs)
 
 
-class Weibull(Continuous):
+class Weibull(PositiveContinuous):
     """
     Weibull log-likelihood
 

pymc3/distributions/distribution.py

@@ -31,12 +31,13 @@ def dist(cls, *args, **kwargs):
         dist.__init__(*args, **kwargs)
         return dist
 
-    def __init__(self, shape, dtype, testval=None, defaults=[]):
+    def __init__(self, shape, dtype, testval=None, defaults=[], transform=None):
         self.shape = np.atleast_1d(shape)
         self.dtype = dtype
         self.type = TensorType(self.dtype, self.shape)
         self.testval = testval
         self.defaults = defaults
+        self.transform = transform
 
     def default(self):
         return self.get_test_val(self.testval, self.defaults)
@@ -80,15 +81,6 @@ def __init__(self, logp, shape=(), dtype='float64',testval=0, *args, **kwargs):
         super(DensityDist, self).__init__(shape, dtype, testval, *args, **kwargs)
         self.logp = logp
 
-def TransformedVar(*args, **kwargs):
-    try:
-        model = Model.get_context()
-    except TypeError:
-        raise TypeError("No model on context stack, which is needed to use the Normal('x', 0,1) syntax. Add a 'with model:' block")
-
-    return model.TransformedVar(*args, **kwargs)
-
-
 class MultivariateContinuous(Continuous):
 
     pass

pymc3/distributions/multivariate.py

@@ -3,6 +3,7 @@
 from .dist_math import *
 
 import numpy as np
+from .transforms import simplextransform
 
 from theano.tensor.nlinalg import det, matrix_inverse, trace, eigh
 from theano.tensor import dot, cast, eye, diag, eq, le, ge, gt, all
@@ -65,8 +66,8 @@ class Dirichlet(Continuous):
         Only the first `k-1` elements of `x` are expected. Can be used
         as a parent of Multinomial and Categorical nevertheless.
     """
-    def __init__(self, a, *args, **kwargs):
-        super(Dirichlet, self).__init__(*args, **kwargs)
+    def __init__(self, a, transform=simplextransform, *args, **kwargs):
+        super(Dirichlet, self).__init__(transform=transform, *args, **kwargs)
         self.a = a
         self.k = a.shape[0]
         self.mean = a / sum(a)

pymc3/distributions/transforms.py

@@ -38,7 +38,7 @@ def __init__(self, dist, transform, *args, **kwargs):
         testval = forward(dist.default())
 
         self.dist = dist
-        self.transform = transform
+        self.transform_used = transform
         v = forward(FreeRV(name='v', distribution=dist))
         self.type = v.type
 
@@ -49,12 +49,42 @@ def __init__(self, dist, transform, *args, **kwargs):
 
 
     def logp(self, x):
-        return self.dist.logp(self.transform.backward(x)) + self.transform.jacobian_det(x)
+        return self.dist.logp(self.transform_used.backward(x)) + self.transform_used.jacobian_det(x)
 
 transform = Transform
 
+
 logtransform = transform("log", log, exp, idfn)
 
+
+logistic = t.nnet.sigmoid
+def logistic_jacobian(x):
+    ex = exp(-x)
+    return log(ex/(ex +1)**2)
+
+def logit(x): 
+    return log(x/(1-x))
+logoddstransform = transform("logodds", logit, logistic, logistic_jacobian)
+
+
+def interval_transform(a, b):
+    def interval_real(x):
+        r= log((x-a)/(b-x))
+        return r
+
+    def real_interval(x):
+        r =  (b-a)*exp(x)/(1+exp(x)) + a
+        return r
+
+    def real_interval_jacobian(x):
+        ex = exp(-x)
+        jac = log(ex*(b-a)/(ex + 1)**2)
+        return jac
+
+    return transform("interval", interval_real, real_interval, real_interval_jacobian)
+
+
+
 simplextransform = transform("simplex",
                              lambda p: p[:-1],
                              lambda p: concatenate(

pymc3/examples/ARM12_6.py

@@ -50,17 +50,14 @@ def run(n=3000):
         n = 50
     with model:
         start = {'groupmean': obs_means.mean(),
-                 'groupsd': obs_means.std(),
-                 'sd': data.groupby('group').lradon.std().mean(),
+                 'groupsd_interval': 0,
+                 'sd_interval': 0,
                  'means': np.array(obs_means),
                  'floor_m': 0.,
                  }
 
         start = find_MAP(start, [groupmean, sd, floor_m])
-        H = model.fastd2logp()
-        h = np.diag(H(start))
-
-        step = HamiltonianMC(model.vars, h)
+        step = NUTS(model.vars, scaling=start)
 
         trace = sample(n, step, start)
 

pymc3/examples/ARM12_6uranium.py

@@ -53,13 +53,14 @@ def run(n=3000):
         n = 50
     with model:
 
-        start = {'groupmean': obs_means.mean(),
-                 'groupsd': obs_means.std(),
-                 'sd': data.groupby('group').lradon.std().mean(),
+        start = Point({
+                 'groupmean': obs_means.mean(),
+                 'groupsd_interval': 0,
+                 'sd_interval': 0,
                  'means': np.array(obs_means),
                  'u_m': np.array([.72]),
                  'floor_m': 0.,
-                 }
+                 })
 
         start = find_MAP(start, model.vars[:-1])
         H = model.fastd2logp()

pymc3/examples/dirichlet.py

@@ -7,9 +7,7 @@
     k = 5
     a = constant(np.array([2, 3., 4, 2, 2]))
 
-    p, p_m1 = model.TransformedVar(
-        'p', Dirichlet.dist(a, shape=k),
-        simplextransform)
+    p = Dirichlet('p', a, shape=k)
 
     c = Categorical('c', p, observed=np.random.randint(0, k, 5))
 

pymc3/examples/discrete_find_MAP.py

@@ -40,13 +40,13 @@
 
 with model:
     for i in range(n+1):
-        s = {'p': 0.5, 'surv_sim': i}
+        s = {'p_logodds': 0.5, 'surv_sim': i}
         map_est = mc.find_MAP(start=s, vars=model.vars,
                               fmin=mc.starting.optimize.fmin_bfgs)
         print('surv_sim: %i->%i, p: %f->%f, LogP:%f'%(s['surv_sim'],
                                                       map_est['surv_sim'],
-                                                      s['p'],
-                                                      map_est['p'],
+                                                      s['p_logodds'],
+                                                      map_est['p_logodds'],
                                                       model.logp(map_est)))
 
 # Once again because the gradient of `surv_sim` provides no information to the
@@ -58,12 +58,12 @@
 
 with model:
     for i in range(n+1):
-        s = {'p': 0.5, 'surv_sim': i}
+        s = {'p_logodds': 0.0, 'surv_sim': i}
         map_est = mc.find_MAP(start=s, vars=model.vars)
         print('surv_sim: %i->%i, p: %f->%f, LogP:%f'%(s['surv_sim'],
                                                       map_est['surv_sim'],
-                                                      s['p'],
-                                                      map_est['p'],
+                                                      s['p_logodds'],
+                                                      map_est['p_logodds'],
                                                       model.logp(map_est)))
 
 # For most starting values this converges to the maximum log likelihood of
@@ -82,12 +82,12 @@ def bh(*args, **kwargs):
 
 with model:
     for i in range(n+1):
-        s = {'p': 0.5, 'surv_sim': i}
+        s = {'p_logodds': 0.0, 'surv_sim': i}
         map_est = mc.find_MAP(start=s, vars=model.vars, fmin=bh)
         print('surv_sim: %i->%i, p: %f->%f, LogP:%f'%(s['surv_sim'],
                                                       map_est['surv_sim'],
-                                                      s['p'],
-                                                      map_est['p'],
+                                                      s['p_logodds'],
+                                                      map_est['p_logodds'],
                                                       model.logp(map_est)))
 
 # By default `basinhopping` uses a gradient minimization technique,
@@ -96,13 +96,13 @@ def bh(*args, **kwargs):
 
 with model:
     for i in range(n+1):
-        s = {'p': 0.5, 'surv_sim': i}
+        s = {'p_logodds': 0.0, 'surv_sim': i}
         map_est = mc.find_MAP(start=s, vars=model.vars,
                               fmin=bh, minimizer_kwargs={"method": "Powell"})
         print('surv_sim: %i->%i, p: %f->%f, LogP:%f'%(s['surv_sim'],
                                                       map_est['surv_sim'],
-                                                      s['p'],
-                                                      map_est['p'],
+                                                      s['p_logodds'],
+                                                      map_est['p_logodds'],
                                                       model.logp(map_est)))
 
 # Confident in our MAP estimate we can sample from the posterior, making sure