Merge pull request #685 from pymc-devs/refactor_glm

Refactor glm

Author: twiecki

pymc3/examples/glm_linear.py

@@ -1,13 +1,7 @@
-from __future__ import print_function
-
-import numpy as np
 import sys
 
-try:
-    import statsmodels.api as sm
-except ImportError:
-    print("Example requires statsmodels")
-    sys.exit(0)
+import numpy as np
+import scipy.optimize as opt
 
 from pymc3 import *
 
@@ -31,13 +25,12 @@ def run(n=2000):
     import matplotlib.pyplot as plt
 
     with model:
-        trace = sample(n, Slice())
+        start = find_MAP(fmin=opt.fmin_powell)
+        trace = sample(n, Slice(), start=start)
 
     plt.plot(x, y, 'x')
     glm.plot_posterior_predictive(trace)
     # plt.show()
 
 if __name__ == '__main__':
     run()
-
-

pymc3/examples/glm_robust.py

@@ -1,14 +1,6 @@
-from __future__ import print_function
-
 import numpy as np
 import sys
 
-try:
-    import statsmodels.api as sm
-except ImportError:
-    print("Example requires statsmodels")
-    sys.exit(0)
-
 from pymc3 import *
 
 # Generate data
@@ -26,9 +18,9 @@
 data_outlier = dict(x=x, y=y)
 
 with Model() as model:
-    family = glm.families.T(link=glm.links.Identity,
-			                priors={'nu': 1.5,
-				                    'lam': ('sigma', Uniform.dist(0, 20))})
+    family = glm.families.T(#link=glm.families.identity,
+                            priors={'nu': 1.5,
+                                    'lam': Uniform.dist(0, 20)})
     glm.glm('y ~ x', data_outlier, family=family)
 
 def run(n=2000):
@@ -45,5 +37,3 @@ def run(n=2000):
 
 if __name__ == '__main__':
     run()
-
-

pymc3/glm/__init__.py

@@ -1,7 +1,5 @@
 try:
-    import statsmodels.api as sm
-    from . import links
     from . import families
     from .glm import *
 except ImportError:
-    print("Warning: statsmodels and/or patsy not found, not importing glm submodule.")
+    print("Warning: patsy not found, not importing glm submodule.")

pymc3/glm/families.py

@@ -1,24 +1,29 @@
 import numbers
 from copy import copy
 
-try:
-    from statsmodels.genmod.families.family import (Gaussian, Binomial, Poisson)
-except ImportError:
-    Gaussian = None
-    Binomial = None
-    Poisson = None
-
-from .links import *
+import theano.tensor
 from ..model import modelcontext
-from ..distributions import Normal, T, Uniform, Bernoulli, Poisson
+from .. import distributions as pm_dists
 
 __all__ = ['Normal', 'T', 'Binomial', 'Poisson']
 
+# Define link functions
+
+# Hack as assigning a function in the class definition automatically binds it as a method.
+class Identity():
+    def __call__(self, x):
+        return x
+
+identity = Identity()
+logit = theano.tensor.nnet.sigmoid
+inverse = theano.tensor.inv
+log = theano.tensor.log
+
 class Family(object):
     """Base class for Family of likelihood distribution and link functions.
     """
     priors = {}
-    link = Identity
+    link = None
 
     def __init__(self, **kwargs):
         # Overwrite defaults
@@ -29,9 +34,6 @@ def __init__(self, **kwargs):
             else:
                 setattr(self, key, val)
 
-        # Instantiate link function
-        self.link_func = self.link()
-
     def _get_priors(self, model=None):
         """Return prior distributions of the likelihood.
 
@@ -45,7 +47,7 @@ def _get_priors(self, model=None):
             if isinstance(val, numbers.Number):
                 priors[key] = val
             else:
-                priors[key] = model.Var(val[0], val[1])
+                priors[key] = model.Var(key, val)
 
         return priors
 
@@ -61,50 +63,38 @@ def create_likelihood(self, y_est, y_data, model=None):
         """
         priors = self._get_priors(model=model)
         # Wrap y_est in link function
-        priors[self.parent] = self.link_func.theano(y_est)
+        priors[self.parent] = self.link(y_est)
         return self.likelihood('y', observed=y_data, **priors)
 
-    def create_statsmodel_family(self):
-        """Instantiate and return statsmodel family object.
-        """
-        if self.sm_family is None:
-            return None
-        else:
-            return self.sm_family(self.link.sm)
-
     def __repr__(self):
         return """Family {klass}:
     Likelihood   : {likelihood}({parent})
     Priors       : {priors}
     Link function: {link}.""".format(klass=self.__class__, likelihood=self.likelihood.__name__, parent=self.parent, priors=self.priors, link=self.link)
 
 
-class Normal(Family):
-    sm_family = Gaussian
-    link = Identity
-    likelihood = Normal
-    parent = 'mu'
-    priors = {'sd': ('sigma', Uniform.dist(0, 100))}
-
-
 class T(Family):
-    sm_family = Gaussian
-    link = Identity
-    likelihood = T
+    link = identity
+    likelihood = pm_dists.T
     parent = 'mu'
-    priors = {'lam': ('sigma', Uniform.dist(0, 100)),
+    priors = {'lam': pm_dists.HalfCauchy.dist(beta=10),
               'nu': 1}
 
+class Normal(Family):
+    link = identity
+    likelihood = pm_dists.Normal
+    parent = 'mu'
+    priors = {'sd': pm_dists.HalfCauchy.dist(beta=10)}
 
 class Binomial(Family):
-    link = Logit
-    sm_family = Binomial
-    likelihood = Bernoulli
+    link = logit
+    likelihood = pm_dists.Bernoulli
     parent = 'p'
+    priors = {'p': pm_dists.Beta.dist(alpha=1, beta=1)}
 
 
 class Poisson(Family):
-    link = Log
-    sm_family = Poisson
-    likelihood = Poisson
-    parent = ''
+    link = log
+    likelihood = pm_dists.Poisson
+    parent = 'mu'
+    priors = {'mu': pm_dists.HalfCauchy.dist(beta=10)}

pymc3/glm/glm.py

@@ -6,17 +6,14 @@
 import theano
 import pandas as pd
 from collections import defaultdict
-from statsmodels.formula.api import glm as glm_sm
-import statsmodels.api as sm
 from pandas.tools.plotting import scatter_matrix
 
-from . import links
 from . import families
 
 def linear_component(formula, data, priors=None,
                      intercept_prior=None,
                      regressor_prior=None,
-                     init=True, init_vals=None, family=None,
+                     init_vals=None, family=None,
                      model=None):
     """Create linear model according to patsy specification.
 
@@ -35,9 +32,6 @@ def linear_component(formula, data, priors=None,
     regressor_prior : pymc3 distribution
         Prior to use for all regressor(s).
         Default: Normal.dist(mu=0, tau=1.0E-12)
-    init : bool
-        Whether to set the starting values via statsmodels
-        Default: True
     init_vals : dict
         Set starting values externally: parameter -> value
         Default: None
@@ -70,10 +64,8 @@ def linear_component(formula, data, priors=None,
     _, dmatrix = patsy.dmatrices(formula, data)
     reg_names = dmatrix.design_info.column_names
 
-    if init_vals is None and init:
-        init_vals = glm_sm(formula, data, family=family).fit().params
-    else:
-        init_vals = defaultdict(lambda: None)
+    if init_vals is None:
+        init_vals = {}
 
     # Create individual coefficients
     model = modelcontext(model)
@@ -82,13 +74,15 @@ def linear_component(formula, data, priors=None,
     if reg_names[0] == 'Intercept':
         prior = priors.get('Intercept', intercept_prior)
         coeff = model.Var(reg_names.pop(0), prior)
-        coeff.tag.test_value = init_vals['Intercept']
+        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)
-        coeff.tag.test_value = init_vals[reg_name]
+        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))
@@ -114,15 +108,10 @@ def glm(*args, **kwargs):
     regressor_prior : pymc3 distribution
         Prior to use for all regressor(s).
         Default: Normal.dist(mu=0, tau=1.0E-12)
-    init : bool
-        Whether initialize test values via statsmodels
-        Default: True
     init_vals : dict
         Set starting values externally: parameter -> value
         Default: None
-    find_MAP : bool
-        Whether to call find_MAP on non-initialized nodes.
-    family : statsmodels.family
+    family : Family object
         Distribution of likelihood, see pymc3.glm.families
         (init has to be True).
 
@@ -135,7 +124,7 @@ def glm(*args, **kwargs):
     # Logistic regression
     vars = glm('male ~ height + weight',
                data,
-               family=glm.families.Binomial(link=glm.links.Logit))
+               family=glm.families.Binomial(link=glm.families.logit))
     """
 
     model = modelcontext(kwargs.get('model'))
@@ -147,21 +136,10 @@ def glm(*args, **kwargs):
     data = args[1]
     y_data = np.asarray(patsy.dmatrices(formula, data)[0]).T
 
-    # Create GLM
-    kwargs['family'] = family.create_statsmodel_family()
-
     y_est, coeffs = linear_component(*args, **kwargs)
     family.create_likelihood(y_est, y_data)
 
-    # Find vars we have not initialized yet
-    non_init_vars = set(model.vars).difference(set(coeffs))
-    if len(non_init_vars) != 0 and call_find_map:
-        start = find_MAP(vars=non_init_vars)
-
-        for var in non_init_vars:
-            var.tag.test_value = start[var.name]
-
-    return [y_est] + coeffs + list(non_init_vars)
+    return [y_est] + coeffs
 
 
 def plot_posterior_predictive(trace, eval=None, lm=None, samples=30, **kwargs):

pymc3/glm/links.py

@@ -50,7 +50,7 @@ PyMC is tested on Python 2.7 and 3.3 and depends on Theano, NumPy, SciPy, and Ma
 
 ### Optional
 
-The GLM submodule relies on Pandas, Patsy, Statsmodels.
+The GLM submodule relies on Pandas and Patsy.
 
 [`scikits.sparse`](https://github.com/njsmith/scikits-sparse) enables sparse scaling matrices which are useful for large problems. Installation on Ubuntu is easy: