CoCalc Public Filespymc3.ipynb
Author: Harald Schilly
Views : 104
Compute Environment: Ubuntu 18.04 (Deprecated)

# Pymc3: Probabilistic Programming in Python: Bayesian Modeling and Probabilistic Machine Learning

Available in the "Python 3 (Ubuntu Linux)" kernel.

http://docs.pymc.io/intro.html

In [1]:
import pymc3 as pm
pm.__version__

/usr/local/lib/python3.5/dist-packages/h5py/__init__.py:36: FutureWarning: Conversion of the second argument of issubdtype from float to np.floating is deprecated. In future, it will be treated as np.float64 == np.dtype(float).type. from ._conv import register_converters as _register_converters
'3.3'
In [2]:
import numpy as np
import matplotlib.pyplot as plt

# Initialize 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.random.randn(size)
X2 = np.random.randn(size) * 0.2

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

In [3]:
basic_model = pm.Model()

with basic_model:

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

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

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


In [4]:
map_estimate = pm.find_MAP(model=basic_model)
map_estimate

logp = -149.58, ||grad|| = 12.242: 100%|██████████| 19/19 [00:00<00:00, 792.13it/s] 
{'alpha': array(0.90660093), 'beta': array([0.94848596, 2.60711845]), 'sigma': array(0.96298858), 'sigma_log__': array(-0.03771373)}
In [5]:
with basic_model:
# draw 500 posterior samples
trace = pm.sample()

Auto-assigning NUTS sampler... Initializing NUTS using jitter+adapt_diag... Sequential sampling (2 chains in 1 job) NUTS: [sigma_log__, beta, alpha] 100%|██████████| 1000/1000 [00:01<00:00, 630.44it/s] 100%|██████████| 1000/1000 [00:01<00:00, 581.02it/s] 
In [6]:
trace['alpha'][-5:]

array([0.9572314 , 0.93554263, 0.86782858, 0.97563444, 0.70462067])
In [7]:
_ = pm.traceplot(trace)

In [8]:
pm.summary(trace)

mean sd mc_error hpd_2.5 hpd_97.5 n_eff Rhat
alpha 0.904619 0.097804 0.002686 0.720830 1.092549 1000.0 0.999163
beta__0 0.947115 0.092629 0.002394 0.744485 1.128060 1000.0 1.001146
beta__1 2.608587 0.538713 0.016015 1.602050 3.685764 1000.0 0.999511
sigma 0.989664 0.072976 0.001916 0.857547 1.140489 1000.0 0.999216
In [ ]: