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
%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');
from pymc3 import Model, Normal, HalfNormal
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)
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) 543 input.type.broadcastable, 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,).
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
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)}
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])
from pymc3 import traceplot
traceplot(trace);
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
type(basic_model)
pymc3.model.Model
type(trace)
pymc3.backends.base.MultiTrace
dir(trace)
['__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']
trace.varnames
['alpha', 'beta', 'sigma_log', 'sigma']
trace.get_values('alpha')
array([ 1.00911214, 1.00911214, 1.0349891 , ..., 1.17301994, 1.01888131, 0.94337901])
dir(basic_model)
['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__', 'add_random_variable', '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']
basic_model.observed_RVs
[Y_obs]
basic_model.unobserved_RVs
[alpha, beta, sigma_log, sigma]