In [1]:
import numpy as np
import matplotlib.pyplot as plt

# Intialize random number generator
np.random.seed(123)

# True parameter values
alpha, sigma = 1, 1
beta = [1, 2.5]

# Size of dataset
size = 100

# Predictor variable
X1 = np.linspace(0, 1, size)
X2 = np.linspace(0,.2, size)

# Simulate outcome variable
Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

In [2]:
%matplotlib inline

fig, axes = plt.subplots(1, 2, sharex=True, figsize=(10,4))
axes[0].scatter(X1, Y)
axes[1].scatter(X2, Y)
axes[0].set_ylabel('Y'); axes[0].set_xlabel('X1'); axes[1].set_xlabel('X2');

In [3]:
from pymc3 import Model, Normal, HalfNormal

In [4]:
basic_model = Model()

with basic_model:

# Priors for unknown model parameters
alpha = Normal('alpha', mu=0, sd=10)
beta = Normal('beta', mu=0, sd=10, shape=2)
sigma = HalfNormal('sigma', sd=1)

# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2

# Likelihood (sampling distribution) of observations
Y_obs = Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

In [22]:
basic_model = Model()

with basic_model:

# Priors for unknown model parameters
alpha = Normal('alpha', mu=0, sd=10)
beta = Normal('beta', mu=0, sd=10, shape=2)
sigma = HalfNormal('sigma', sd=1)

# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2

# Likelihood (sampling distribution) of observations
Y = Normal('Y_obs', mu=mu, sd=sigma)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-22-c25c4c7ed7fb> in <module>()
12
13     # Likelihood (sampling distribution) of observations
---> 14     Y = Normal('Y_obs', mu=mu, sd=sigma)

/home/downey/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
20             data = kwargs.pop('observed', None)
21             dist = cls.dist(*args, **kwargs)
---> 22             return model.Var(name, dist, data)
23         elif name is None:
24             return object.__new__(cls) #for pickle

/home/downey/anaconda/lib/python2.7/site-packages/pymc3/model.pyc in Var(self, name, dist, data)
224         if data is None:
225             if getattr(dist, "transform", None) is None:
--> 226                 var = FreeRV(name=name, distribution=dist, model=self)
227                 self.free_RVs.append(var)
228             else:

/home/downey/anaconda/lib/python2.7/site-packages/pymc3/model.pyc in __init__(self, type, owner, index, name, distribution, model)
436             self.tag.test_value = np.ones(
437                 distribution.shape, distribution.dtype) * distribution.default()
--> 438             self.logp_elemwiset = distribution.logp(self)
439             self.model = model
440

/home/downey/anaconda/lib/python2.7/site-packages/pymc3/distributions/continuous.pyc in logp(self, value)
171
172         return bound(
--> 173             (-tau * (value - mu) ** 2 + log(tau / pi / 2.)) / 2.,
174             tau > 0,
175             sd > 0

/home/downey/anaconda/lib/python2.7/site-packages/theano/tensor/var.pyc in __sub__(self, other)
145         # and the return value in that case
146         try:
--> 147             return theano.tensor.basic.sub(self, other)
148         except (NotImplementedError, AsTensorError):
149             return NotImplemented

/home/downey/anaconda/lib/python2.7/site-packages/theano/gof/op.pyc in __call__(self, *inputs, **kwargs)
505         """
506         return_list = kwargs.pop('return_list', False)
--> 507         node = self.make_node(*inputs, **kwargs)
508
509         if config.compute_test_value != 'off':

/home/downey/anaconda/lib/python2.7/site-packages/theano/tensor/elemwise.pyc in make_node(self, *inputs)
544                     ['x'] * difference + range(length),
--> 545                     inplace=False)(input))
546         inputs = args
547

/home/downey/anaconda/lib/python2.7/site-packages/theano/gof/op.pyc in __call__(self, *inputs, **kwargs)
515             for i, ins in enumerate(node.inputs):
516                 try:
--> 517                     storage_map[ins] = [self._get_test_value(ins)]
518                     compute_map[ins] = [True]
519                 except AttributeError:

/home/downey/anaconda/lib/python2.7/site-packages/theano/gof/op.pyc in _get_test_value(cls, v)
452             # ensure that the test value is correct
453             try:
--> 454                 ret = v.type.filter(v.tag.test_value)
455             except Exception, e:
456                 # Better error message.

/home/downey/anaconda/lib/python2.7/site-packages/theano/tensor/type.pyc in filter(self, data, strict, allow_downcast)
167             raise TypeError("Wrong number of dimensions: expected %s,"
168                             " got %s with shape %s." % (self.ndim, data.ndim,
--> 169                                                         data.shape))
170         if not data.flags.aligned:
171             try:

TypeError: For compute_test_value, one input test value does not have the requested type.

The error when converting the test value to that variable type:
Wrong number of dimensions: expected 0, got 1 with shape (100,).

In [23]:
alpha = Normal('alpha', mu=0, sd=10)
beta = Normal('beta', mu=0, sd=10, shape=2)
sigma = HalfNormal('sigma', sd=1)

# Expected value of outcome
mu = alpha + beta[0]*X1 + beta[1]*X2

# Likelihood (sampling distribution) of observations
Y = Normal('Y_obs', mu=mu, sd=sigma)

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-23-0fffde8373ce> in <module>()
----> 1 alpha = Normal('alpha', mu=0, sd=10)
2 beta = Normal('beta', mu=0, sd=10, shape=2)
3 sigma = HalfNormal('sigma', sd=1)
4
5 # Expected value of outcome

/home/downey/anaconda/lib/python2.7/site-packages/pymc3/distributions/distribution.pyc in __new__(cls, name, *args, **kwargs)
15             model = Model.get_context()
16         except TypeError:
---> 17             raise TypeError("No model on context stack, which is needed to use the Normal('x', 0,1) syntax. Add a 'with model:' block")
18
19         if isinstance(name, str):

TypeError: No model on context stack, which is needed to use the Normal('x', 0,1) syntax. Add a 'with model:' block

In [5]:
from pymc3 import find_MAP

map_estimate = find_MAP(model=basic_model)

print(map_estimate)

{'alpha': array(1.01366409951285), 'beta': array([ 1.46791595,  0.29358319]), 'sigma_log': array(0.11928764983495602)}

In [7]:
from scipy import optimize
from pymc3 import NUTS, sample

with basic_model:

# obtain starting values via MAP
start = find_MAP(fmin=optimize.fmin_powell)

# instantiate sampler
step = NUTS(scaling=start)

# draw 2000 posterior samples
trace = sample(2000, step, start=start)

[-----------------100%-----------------] 2000 of 2000 complete in 5.4 sec
/home/downey/anaconda/lib/python2.7/site-packages/theano/gof/cmodule.py:293: RuntimeWarning: numpy.ndarray size changed, may indicate binary incompatibility
rval = __import__(module_name, {}, {}, [module_name])

In [8]:
from pymc3 import traceplot

traceplot(trace);

In [9]:
from pymc3 import summary

summary(trace)

alpha:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------

0.999            0.234            0.007            [0.535, 1.424]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

0.543          0.844          1.002          1.153          1.454

beta:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------

1.638            2.107            0.128            [-2.402, 5.976]
-0.436           10.226           0.622            [-21.629, 19.058]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

-2.333         0.188          1.580          2.918          6.062
-21.700        -6.718         -0.462         6.751          19.007

sigma_log:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------

0.141            0.077            0.002            [-0.004, 0.293]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

-0.001         0.087          0.140          0.191          0.300

sigma:

Mean             SD               MC Error         95% HPD interval
-------------------------------------------------------------------

1.155            0.089            0.002            [0.989, 1.332]

Posterior quantiles:
2.5            25             50             75             97.5
|--------------|==============|==============|--------------|

0.999          1.091          1.150          1.210          1.349

In [13]:
type(basic_model)

Out[13]:
pymc3.model.Model

In [14]:
type(trace)

Out[14]:
pymc3.backends.base.MultiTrace

In [15]:
dir(trace)

Out[15]:
['__class__',
'__delattr__',
'__dict__',
'__doc__',
'__format__',
'__getattr__',
'__getattribute__',
'__getitem__',
'__hash__',
'__init__',
'__len__',
'__module__',
'__new__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'_attrs',
'_slice',
'_straces',
'chains',
'get_values',
'nchains',
'point',
'varnames']

In [17]:
trace.varnames

Out[17]:
['alpha', 'beta', 'sigma_log', 'sigma']

In [18]:
trace.get_values('alpha')

Out[18]:
array([ 1.00911214,  1.00911214,  1.0349891 , ...,  1.17301994,
1.01888131,  0.94337901])

In [19]:
dir(basic_model)

Out[19]:
['Var',
'Y_obs',
'__class__',
'__delattr__',
'__dict__',
'__doc__',
'__enter__',
'__exit__',
'__format__',
'__getattribute__',
'__getitem__',
'__hash__',
'__init__',
'__module__',
'__new__',
'__reduce__',
'__reduce_ex__',
'__repr__',
'__setattr__',
'__sizeof__',
'__str__',
'__subclasshook__',
'__weakref__',
'alpha',
'basic_RVs',
'beta',
'cont_vars',
'contexts',
'd2logp',
'deterministics',
'disc_vars',
'dlogp',
'fastd2logp',
'fastdlogp',
'fastfn',
'fastlogp',
'fn',
'free_RVs',
'get_context',
'get_contexts',
'logp',
'logp_elemwise',
'logpt',
'makefn',
'missing_values',
'model',
'named_vars',
'observed_RVs',
'potentials',
'profile',
'sigma_log',
'test_point',
'unobserved_RVs',
'vars']

In [20]:
basic_model.observed_RVs

Out[20]:
[Y_obs]

In [21]:
basic_model.unobserved_RVs

Out[21]:
[alpha, beta, sigma_log, sigma]

In [ ]: