Merge branch 'constant' of github.com:pymc-devs/pymc3 into constant

Author: fonnesbeck

docs/source/notebooks/NUTS_scaling_using_ADVI.ipynb

@@ -490,7 +490,7 @@
    ],
    "source": [
     "with mdl:\n",
-    "    v_params = pm.variational.advi(n=100000, verbose=False) \n",
+    "    v_params = pm.variational.advi(n=100000) \n",
     "\n",
     "_ = plt.plot(-np.log10(-v_params.elbo_vals))"
    ]

pymc3/distributions/transforms.py

@@ -240,7 +240,7 @@ def jacobian_det(self, y_):
 t_stick_breaking = lambda eps: StickBreaking(eps)
 
 
-class Circular(ElemwiseTransform):
+class Circular(Transform):
     """Transforms a linear space into a circular one.
     """
     name = "circular"
@@ -250,5 +250,8 @@ def backward(self, y):
 
     def forward(self, x):
         return x
+        
+    def jacobian_det(self, x):
+        return 0
 
 circular = Circular()

pymc3/glm/families.py

@@ -39,24 +39,26 @@ def __init__(self, **kwargs):
             else:
                 setattr(self, key, val)
 
-    def _get_priors(self, model=None):
+    def _get_priors(self, model=None, name=''):
         """Return prior distributions of the likelihood.
 
         Returns
         -------
         dict : mapping name -> pymc3 distribution
         """
+        if name:
+            name = '{}_'.format(name)
         model = modelcontext(model)
         priors = {}
         for key, val in self.priors.items():
             if isinstance(val, numbers.Number):
                 priors[key] = val
             else:
-                priors[key] = model.Var(key, val)
+                priors[key] = model.Var('{}{}'.format(name, key), val)
 
         return priors
 
-    def create_likelihood(self, y_est, y_data, model=None):
+    def create_likelihood(self, name, y_est, y_data, model=None):
         """Create likelihood distribution of observed data.
 
         Parameters
@@ -66,10 +68,12 @@ def create_likelihood(self, y_est, y_data, model=None):
         y_data : array
             Observed dependent variable
         """
-        priors = self._get_priors(model=model)
+        priors = self._get_priors(model=model, name=name)
         # Wrap y_est in link function
         priors[self.parent] = self.link(y_est)
-        return self.likelihood('y', observed=y_data, **priors)
+        if name:
+            name = '{}_'.format(name)
+        return self.likelihood('{}y'.format(name), observed=y_data, **priors)
 
     def __repr__(self):
         return """Family {klass}:

pymc3/glm/glm.py

@@ -1,23 +1,41 @@
 import numpy as np
 from ..distributions import Normal
-from ..tuning.starting import find_MAP
 from ..model import modelcontext
 import patsy
-import theano
 import pandas as pd
-from collections import defaultdict
-from pandas.tools.plotting import scatter_matrix
+import theano
+from collections import defaultdict, namedtuple
 
 from . import families
 
 __all__ = ['glm', 'linear_component', 'plot_posterior_predictive']
 
 
-def linear_component(formula, data, priors=None,
-                     intercept_prior=None,
-                     regressor_prior=None,
-                     init_vals=None, family=None,
-                     model=None):
+def _xy_to_data_and_formula(X, y):
+    if not isinstance(y, pd.Series):
+        y = pd.Series(y, name='y')
+    else:
+        if not y.name:
+            y.name = 'y'
+    if not isinstance(X, (pd.DataFrame, pd.Series)):
+        if len(X.shape) > 1:
+            cols = ['x%d' % i for i in range(X.shape[1])]
+        else:
+            cols = ['x']
+        X = pd.DataFrame(X, columns=cols)
+    elif isinstance(X, pd.Series):
+        if not X.name:
+            X.name = 'x'
+    # else -> pd.DataFrame -> ok
+    data = pd.concat([y, X], 1)
+    formula = patsy.ModelDesc(
+        [patsy.Term([patsy.LookupFactor(y.name)])],
+        [patsy.Term([patsy.LookupFactor(p)]) for p in X.columns]
+    )
+    return data, formula
+
+
+class linear_component(namedtuple('Estimate', 'y_est,coeffs')):
     """Create linear model according to patsy specification.
 
     Parameters
@@ -38,10 +56,6 @@ def linear_component(formula, data, priors=None,
     init_vals : dict
         Set starting values externally: parameter -> value
         Default: None
-    family : statsmodels.family
-        Link function to pass to statsmodels (init has to be True).
-    See `statsmodels.api.families`
-        Default: identity
 
     Output
     ------
@@ -50,51 +64,78 @@ def linear_component(formula, data, priors=None,
     Example
     -------
     # Logistic regression
-    y_est, coeffs = glm('male ~ height + weight',
-                        htwt_data,
-                        family=glm.families.Binomial(link=glm.family.logit))
-    y_data = Bernoulli('y', y_est, observed=data.male)
+    y_est, coeffs = linear_component('male ~ height + weight',
+                        htwt_data)
+    probability = glm.families.logit(y_est)
+    y_data = Bernoulli('y', probability, observed=data.male)
     """
-    if intercept_prior is None:
-        intercept_prior = Normal.dist(mu=0, tau=1.0E-12)
-    if regressor_prior is None:
-        regressor_prior = Normal.dist(mu=0, tau=1.0E-12)
-
-    if priors is None:
-        priors = defaultdict(None)
-
-    # Build patsy design matrix and get regressor names.
-    _, dmatrix = patsy.dmatrices(formula, data)
-    reg_names = dmatrix.design_info.column_names
-
-    if init_vals is None:
-        init_vals = {}
-
-    # Create individual coefficients
-    model = modelcontext(model)
-    coeffs = []
-
-    if reg_names[0] == 'Intercept':
-        prior = priors.get('Intercept', intercept_prior)
-        coeff = model.Var(reg_names.pop(0), prior)
-        if 'Intercept' in init_vals:
-            coeff.tag.test_value = init_vals['Intercept']
-        coeffs.append(coeff)
-
-    for reg_name in reg_names:
-        prior = priors.get(reg_name, regressor_prior)
-        coeff = model.Var(reg_name, prior)
-        if reg_name in init_vals:
-            coeff.tag.test_value = init_vals[reg_name]
-        coeffs.append(coeff)
-
-    y_est = theano.dot(np.asarray(dmatrix),
-                       theano.tensor.stack(*coeffs)).reshape((1, -1))
+    __slots__ = ()
 
-    return y_est, coeffs
-
-
-def glm(*args, **kwargs):
+    def __new__(cls, formula, data, priors=None,
+                     intercept_prior=None,
+                     regressor_prior=None,
+                     init_vals=None,
+                     model=None,
+                     name=''):
+        if intercept_prior is None:
+            intercept_prior = Normal.dist(mu=0, tau=1.0E-12)
+        if regressor_prior is None:
+            regressor_prior = Normal.dist(mu=0, tau=1.0E-12)
+
+        if priors is None:
+            priors = defaultdict(None)
+
+        # Build patsy design matrix and get regressor names.
+        _, dmatrix = patsy.dmatrices(formula, data)
+        reg_names = dmatrix.design_info.column_names
+
+        if init_vals is None:
+            init_vals = {}
+
+        # Create individual coefficients
+        model = modelcontext(model)
+        coeffs = []
+        if name:
+            name = '{}_'.format(name)
+        if reg_names[0] == 'Intercept':
+            prior = priors.get('Intercept', intercept_prior)
+            coeff = model.Var('{}{}'.format(name, reg_names.pop(0)), prior)
+            if 'Intercept' in init_vals:
+                coeff.tag.test_value = init_vals['Intercept']
+            coeffs.append(coeff)
+
+        for reg_name in reg_names:
+            prior = priors.get(reg_name, regressor_prior)
+            coeff = model.Var('{}{}'.format(name, reg_name), prior)
+            if reg_name in init_vals:
+                coeff.tag.test_value = init_vals[reg_name]
+            coeffs.append(coeff)
+
+        y_est = theano.dot(np.asarray(dmatrix),
+                           theano.tensor.stack(*coeffs)).reshape((1, -1))
+
+        return super(linear_component, cls).__new__(cls, y_est, coeffs)
+
+    @classmethod
+    def from_xy(cls, X, y,
+                priors=None,
+                intercept_prior=None,
+                regressor_prior=None,
+                init_vals=None,
+                model=None,
+                name=''):
+        data, formula = _xy_to_data_and_formula(X, y)
+        return cls(formula, data,
+                   priors=priors,
+                   intercept_prior=intercept_prior,
+                   regressor_prior=regressor_prior,
+                   init_vals=init_vals,
+                   model=model,
+                   name=name
+                   )
+
+
+class glm(namedtuple('Estimate', 'y_est,coeffs')):
     """Create GLM after Patsy model specification string.
 
     Parameters
@@ -121,29 +162,66 @@ def glm(*args, **kwargs):
 
     Output
     ------
-    vars : List of created random variables (y_est, coefficients etc)
+    (y_est, coeffs) : Estimate for y, list of coefficients
 
     Example
     -------
     # Logistic regression
     vars = glm('male ~ height + weight',
                data,
-               family=glm.families.Binomial(link=glm.families.logit))
+               family=glm.families.Binomial())
     """
-
-    model = modelcontext(kwargs.get('model'))
-
-    family = kwargs.pop('family', families.Normal())
-
-    call_find_map = kwargs.pop('find_MAP', True)
-    formula = args[0]
-    data = args[1]
-    y_data = np.asarray(patsy.dmatrices(formula, data)[0]).T
-
-    y_est, coeffs = linear_component(*args, **kwargs)
-    family.create_likelihood(y_est, y_data)
-
-    return [y_est] + coeffs
+    __slots__ = ()
+
+    def __new__(cls, formula, data, priors=None,
+            intercept_prior=None,
+            regressor_prior=None,
+            init_vals=None,
+            family='normal',
+            model=None,
+            name=''):
+        _families = dict(
+            normal=families.Normal,
+            student=families.StudentT,
+            binomial=families.Binomial,
+            poisson=families.Poisson
+        )
+        if isinstance(family, str):
+            family = _families[family]()
+
+        y_data = np.asarray(patsy.dmatrices(formula, data)[0]).T
+
+        y_est, coeffs = linear_component(
+            formula, data, priors=priors,
+            intercept_prior=intercept_prior,
+            regressor_prior=regressor_prior,
+            init_vals=init_vals,
+            model=model,
+            name=name
+            )
+        family.create_likelihood(name, y_est, y_data, model=model)
+
+        return super(glm, cls).__new__(cls, y_est, coeffs)
+
+    @classmethod
+    def from_xy(cls, X, y,
+                priors=None,
+                intercept_prior=None,
+                regressor_prior=None,
+                init_vals=None,
+                family='normal',
+                model=None,
+                name=''):
+        data, formula = _xy_to_data_and_formula(X, y)
+        return cls(formula, data,
+                   priors=priors,
+                   intercept_prior=intercept_prior,
+                   regressor_prior=regressor_prior,
+                   init_vals=init_vals,
+                   model=model,
+                   family=family,
+                   name=name
+                   )
 
 
 def plot_posterior_predictive(trace, eval=None, lm=None, samples=30, **kwargs):

pymc3/model.py

@@ -169,25 +169,16 @@ def logpt(self):
 
 
 class Model(Context, Factor):
-    """Encapsulates the variables and likelihood factors of a model.
+    """Encapsulates the variables and likelihood factors of a model."""
 
-    Parameters
-    ----------
-    verbose : int
-        Model verbosity setting, determining how much feedback various
-        operations provide. Normal verbosity is verbose=1 (default), silence
-        is verbose=0, high is any value greater than 1.
-    """
-
-    def __init__(self, verbose=1):
+    def __init__(self):
         self.named_vars = {}
         self.free_RVs = []
         self.observed_RVs = []
         self.deterministics = []
         self.potentials = []
         self.missing_values = []
         self.model = self
-        self.verbose = verbose
 
     @property
     @memoize
@@ -287,12 +278,11 @@ def Var(self, name, dist, data=None):
             else:
                 var = TransformedRV(name=name, distribution=dist, model=self,
                                     transform=dist.transform)
-                if self.verbose:
-                    pm._log.info('Applied {transform}-transform to {name}'
-                                 ' and added transformed {orig_name} to model.'.format(
-                                    transform=dist.transform.name,
-                                    name=name,
-                                    orig_name='{}_{}_'.format(name, dist.transform.name)))
+                pm._log.debug('Applied {transform}-transform to {name}'
+                              ' and added transformed {orig_name} to model.'.format(
+                                transform=dist.transform.name,
+                                name=name,
+                                orig_name='{}_{}_'.format(name, dist.transform.name)))
                 self.deterministics.append(var)
                 return var
         elif isinstance(data, dict):
@@ -624,8 +614,8 @@ def __init__(self, name, data, distribution, model):
         self.data = {name: as_tensor(data, name, model, distribution)
                      for name, data in data.items()}
 
-        self.missing_values = [data.missing_values for data in self.data.values()
-                               if data.missing_values is not None]
+        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)
         self.model = model
         self.distribution = distribution