Book a Demo!
CoCalc Logo Icon
StoreFeaturesDocsShareSupportNewsAboutPoliciesSign UpSign In
Download
353 views
Tweet
Share
Share
Kernel: Python 3 (Ubuntu Linux)

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

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'
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

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)


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)}
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] 
trace['alpha'][-5:]
array([0.9572314 , 0.93554263, 0.86782858, 0.97563444, 0.70462067])
_ = pm.traceplot(trace)
Image in a Jupyter notebook
pm.summary(trace)