CoCalc Public Filespymc3.ipynbOpen with one click!
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 [ ]: