Added base class for variational methods

pymc3/data.py

@@ -1,7 +1,10 @@
+import itertools
 import pkgutil
 import io
 
-__all__ = ['get_data_file']
+import theano.tensor as tt
+
+__all__ = ['get_data_file', 'DataGenerator']
 
 
 def get_data_file(pkg, path):
@@ -19,3 +22,28 @@ def get_data_file(pkg, path):
     """
 
     return io.BytesIO(pkgutil.get_data(pkg, path))
+
+
+class DataGenerator(object):
+    """
+    Helper class that helps to infer data type of generator with looking
+    at the first item, preserving the order of the resulting generator
+    """
+    def __init__(self, generator):
+        self.test_value = next(generator)
+        self.gen = itertools.chain([self.test_value], generator)
+        self.tensortype = tt.TensorType(
+            self.test_value.dtype,
+            ((False, ) * self.test_value.ndim))
+
+    def __next__(self):
+        return next(self.gen)
+
+    def __iter__(self):
+        return self
+
+    def __eq__(self, other):
+        return id(self) == id(other)
+
+    def __hash__(self):
+        return hash(id(self))

pymc3/distributions/dist_math.py

@@ -10,6 +10,9 @@
 
 from .special import gammaln, multigammaln
 
+c = - 0.5 * np.log(2 * np.pi)
+
+
 def bound(logp, *conditions, **kwargs):
     """
     Bounds a log probability density with several conditions.
@@ -95,3 +98,64 @@ def i1(x):
                      x**9 / 1474560 + x**11 / 176947200 + x**13 / 29727129600,
                      np.e**x / (2 * np.pi * x)**0.5 * (1 - 3 / (8 * x) + 15 / (128 * x**2) + 315 / (3072 * x**3)
                                                        + 14175 / (98304 * x**4)))
+
+
+def sd2rho(sd):
+    """
+    `sd -> rho` theano converter
+
+    :math:`mu + sd*e = mu + log(1+exp(rho))*e`"""
+    return tt.log(tt.exp(sd) - 1)
+
+
+def rho2sd(rho):
+    """
+    `rho -> sd` theano converter
+
+    :math:`mu + sd*e = mu + log(1+exp(rho))*e`"""
+    return tt.log1p(tt.exp(rho))
+
+
+def log_normal(x, mean, **kwargs):
+    """
+    Calculate logarithm of normal distribution at point `x`
+    with given `mean` and `std`
+
+    Parameters
+    ----------
+    x : Tensor
+        point of evaluation
+    mean : Tensor
+        mean of normal distribution
+    kwargs : one of parameters `{sd, tau, w, rho}`
+
+    Notes
+    -----
+    There are four variants for density parametrization.
+    They are:
+        1) standard deviation - `std`
+        2) `w`, logarithm of `std` :math:`w = log(std)`
+        3) `rho` that follows this equation :math:`rho = log(exp(std) - 1)`
+        4) `tau` that follows this equation :math:`tau = std^{-1}`
+    ----
+    """
+    sd = kwargs.get('sd')
+    w = kwargs.get('w')
+    rho = kwargs.get('rho')
+    tau = kwargs.get('tau')
+    eps = kwargs.get('eps', 0.0)
+    check = sum(map(lambda a: a is not None, [sd, w, rho, tau]))
+    if check > 1:
+        raise ValueError('more than one required kwarg is passed')
+    if check == 0:
+        raise ValueError('none of required kwarg is passed')
+    if sd is not None:
+        std = sd
+    elif w is not None:
+        std = tt.exp(w)
+    elif rho is not None:
+        std = rho2sd(rho)
+    else:
+        std = tau**(-1)
+    std += eps
+    return c - tt.log(tt.abs_(std)) - (x - mean) ** 2 / (2 * std ** 2)

pymc3/distributions/distribution.py

@@ -8,7 +8,7 @@
 from .dist_math import bound
 
 
-__all__ = ['DensityDist', 'Distribution', 'Continuous', 'Bound', 
+__all__ = ['DensityDist', 'Distribution', 'Continuous', 'Bound',
            'Discrete', 'NoDistribution', 'TensorType', 'draw_values']
 
 
@@ -30,8 +30,9 @@ def __new__(cls, name, *args, **kwargs):
 
         if isinstance(name, string_types):
             data = kwargs.pop('observed', None)
+            total_size = kwargs.pop('total_size', None)
             dist = cls.dist(*args, **kwargs)
-            return model.Var(name, dist, data)
+            return model.Var(name, dist, data, total_size)
         else:
             raise TypeError("Name needs to be a string but got: %s" % name)
 
@@ -410,7 +411,7 @@ def __init__(self, distribution, lower, upper, transform='infer', *args, **kwarg
                 self.transform = transforms.upperbound(upper)
                 if default >= upper:
                     self.testval = upper - 1
-                    
+
         if issubclass(distribution, Discrete):
             self.transform = None
 

pymc3/math.py

@@ -22,3 +22,7 @@ def invlogit(x, eps=sys.float_info.epsilon):
 
 def logit(p):
     return tt.log(p / (1 - p))
+
+
+def flatten_list(tensors):
+    return tt.concatenate([var.ravel() for var in tensors])

pymc3/model.py

@@ -8,7 +8,7 @@
 
 import pymc3 as pm
 from .memoize import memoize
-from .theanof import gradient, hessian, inputvars
+from .theanof import gradient, hessian, inputvars, generator
 from .vartypes import typefilter, discrete_types, continuous_types
 from .blocking import DictToArrayBijection, ArrayOrdering
 
@@ -458,7 +458,7 @@ def cont_vars(self):
         """All the continuous variables in the model"""
         return list(typefilter(self.vars, continuous_types))
 
-    def Var(self, name, dist, data=None):
+    def Var(self, name, dist, data=None, total_size=None):
         """Create and add (un)observed random variable to the model with an
         appropriate prior distribution.
 
@@ -469,6 +469,8 @@ def Var(self, name, dist, data=None):
         data : array_like (optional)
            If data is provided, the variable is observed. If None,
            the variable is unobserved.
+        total_size : scalar
+            upscales logp of variable with :math:`coef = total_size/var.shape[0]`
 
         Returns
         -------
@@ -477,11 +479,13 @@ def Var(self, name, dist, data=None):
         name = self.name_for(name)
         if data is None:
             if getattr(dist, "transform", None) is None:
-                var = FreeRV(name=name, distribution=dist, model=self)
+                var = FreeRV(name=name, distribution=dist, model=self,
+                             total_size=total_size)
                 self.free_RVs.append(var)
             else:
                 var = TransformedRV(name=name, distribution=dist, model=self,
-                                    transform=dist.transform)
+                                    transform=dist.transform,
+                                    total_size=total_size)
                 pm._log.debug('Applied {transform}-transform to {name}'
                               ' and added transformed {orig_name} to model.'.format(
                                 transform=dist.transform.name,
@@ -491,7 +495,7 @@ def Var(self, name, dist, data=None):
                 return var
         elif isinstance(data, dict):
             var = MultiObservedRV(name=name, data=data, distribution=dist,
-                                  model=self)
+                                  model=self, total_size=total_size)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs += var.missing_values
@@ -500,7 +504,8 @@ def Var(self, name, dist, data=None):
                     self.named_vars[v.name] = v
         else:
             var = ObservedRV(name=name, data=data,
-                             distribution=dist, model=self)
+                             distribution=dist, model=self,
+                             total_size=total_size)
             self.observed_RVs.append(var)
             if var.missing_values:
                 self.free_RVs.append(var.missing_values)
@@ -717,15 +722,18 @@ class FreeRV(Factor, TensorVariable):
     """Unobserved random variable that a model is specified in terms of."""
 
     def __init__(self, type=None, owner=None, index=None, name=None,
-                 distribution=None, model=None):
+                 distribution=None, model=None, total_size=None):
         """
         Parameters
         ----------
         type : theano type (optional)
         owner : theano owner (optional)
         name : str
         distribution : Distribution
-        model : Model"""
+        model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
+        """
         if type is None:
             type = distribution.type
         super(FreeRV, self).__init__(type, owner, index, name)
@@ -736,7 +744,14 @@ def __init__(self, type=None, owner=None, index=None, name=None,
             self.distribution = distribution
             self.tag.test_value = np.ones(
                 distribution.shape, distribution.dtype) * distribution.default()
-            self.logp_elemwiset = distribution.logp(self)
+            logp_elemwiset = distribution.logp(self)
+            if total_size is None:
+                coef = tt.as_tensor(1)
+            else:
+                assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                                  'needs to be at least 1 dimensional')
+                coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
+            self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
 
             incorporate_methods(source=distribution, destination=self,
@@ -759,6 +774,8 @@ def pandas_to_array(data):
         return data
     elif isinstance(data, theano.gof.graph.Variable):
         return data
+    elif hasattr(data, '__next__'):
+        return generator(data)
     else:
         return np.asarray(data)
 
@@ -792,7 +809,7 @@ class ObservedRV(Factor, TensorVariable):
     """
 
     def __init__(self, type=None, owner=None, index=None, name=None, data=None,
-                 distribution=None, model=None):
+                 distribution=None, model=None, total_size=None):
         """
         Parameters
         ----------
@@ -801,6 +818,8 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
         name : str
         distribution : Distribution
         model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
         """
         from .distributions import TensorType
         if type is None:
@@ -814,7 +833,14 @@ def __init__(self, type=None, owner=None, index=None, name=None, data=None,
             data = as_tensor(data, name, model, distribution)
             self.missing_values = data.missing_values
 
-            self.logp_elemwiset = distribution.logp(data)
+            logp_elemwiset = distribution.logp(data)
+            if total_size is None:
+                coef = tt.as_tensor(1)
+            else:
+                assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                                  'needs to be at least 1 dimensional')
+                coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
+            self.logp_elemwiset = logp_elemwiset * coef
             self.model = model
             self.distribution = distribution
 
@@ -835,7 +861,7 @@ class MultiObservedRV(Factor):
     Potentially partially observed.
     """
 
-    def __init__(self, name, data, distribution, model):
+    def __init__(self, name, data, distribution, model, total_size=None):
         """
         Parameters
         ----------
@@ -844,14 +870,23 @@ def __init__(self, name, data, distribution, model):
         name : str
         distribution : Distribution
         model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
         """
         self.name = name
         self.data = {name: as_tensor(data, name, model, distribution)
                      for name, data in data.items()}
 
         self.missing_values = [datum.missing_values for datum in self.data.values()
                                if datum.missing_values is not None]
-        self.logp_elemwiset = distribution.logp(**self.data)
+        logp_elemwiset = distribution.logp(**self.data)
+        if total_size is None:
+            coef = tt.as_tensor(1)
+        else:
+            assert logp_elemwiset.ndim >= 1, ('Variable with scaled density '
+                                              'needs to be at least 1 dimensional')
+            coef = tt.as_tensor(total_size) / logp_elemwiset.shape[0]
+        self.logp_elemwiset = logp_elemwiset * coef
         self.model = model
         self.distribution = distribution
 
@@ -896,17 +931,20 @@ def Potential(name, var, model=None):
 class TransformedRV(TensorVariable):
 
     def __init__(self, type=None, owner=None, index=None, name=None,
-                 distribution=None, model=None, transform=None):
+                 distribution=None, model=None, transform=None,
+                 total_size=None):
         """
         Parameters
         ----------
 
         type : theano type (optional)
         owner : theano owner (optional)
-
         name : str
         distribution : Distribution
-        model : Model"""
+        model : Model
+        total_size : scalar Tensor (optional)
+            needed for upscaling logp
+        """
         if type is None:
             type = distribution.type
         super(TransformedRV, self).__init__(type, owner, index, name)
@@ -916,7 +954,7 @@ def __init__(self, type=None, owner=None, index=None, name=None,
 
             transformed_name = "{}_{}_".format(name, transform.name)
             self.transformed = model.Var(
-                transformed_name, transform.apply(distribution))
+                transformed_name, transform.apply(distribution), total_size=total_size)
 
             normalRV = transform.backward(self.transformed)
 

pymc3/tests/test_model.py

@@ -1,5 +1,5 @@
 import unittest
-import theano.tensor as tt
+from theano import theano, tensor as tt
 import pymc3 as pm
 from pymc3.distributions import HalfCauchy, Normal
 from pymc3 import Potential, Deterministic
@@ -118,3 +118,12 @@ def test_model_root(self):
             self.assertTrue(model is model.root)
             with pm.Model() as sub:
                 self.assertTrue(model is sub.root)
+
+    def test_density_scaling(self):
+        with pm.Model() as model1:
+            Normal('n', observed=[[1]], total_size=1)
+            p1 = theano.function([], model1.logpt)
+        with pm.Model() as model2:
+            Normal('n', observed=[[1]], total_size=2)
+            p2 = theano.function([], model2.logpt)
+        self.assertEqual(p1() * 2, p2())