J'ai essayé d'exécuter pymc

environnement

mac python3.5

Installation

pip install pymc

Essayez de déplacer la partie 1

mymodel.py

# Import relevant modules
import pymc
import numpy as np

# Some data
n = 5 * np.ones(4, dtype=int)
x = np.array([-.86, -.3, -.05, .73])

# Priors on unknown parameters
alpha = pymc.Normal('alpha', mu=0, tau=.01)
beta = pymc.Normal('beta', mu=0, tau=.01)

# Arbitrary deterministic function of parameters
@pymc.deterministic
def theta(a=alpha, b=beta):
    """theta = logit^{-1}(a+b)"""
    return pymc.invlogit(a + b * x)

# Binomial likelihood for data
d = pymc.Binomial('d', n=n, p=theta, value=np.array([0., 1., 3., 5.]),
                  observed=True)

main.py

import pymc
import mymodel

S = pymc.MCMC(mymodel, db='pickle')
S.sample(iter=10000, burn=5000, thin=2)
pymc.Matplot.plot(S)

python main.py

[-----------------65%----              ] 6510 of 10000 complete in 0.5 sec
[-----------------100%-----------------] 10000 of 10000 complete in 0.8 sec
Plotting beta
Plotting alpha
Plotting theta_0
Plotting theta_1
Plotting theta_2
Plotting theta_3

Essayez de déplacer la partie 2

pip install triangle-plot
pip install corner

jupyter notebook

Commencez par créer de fausses données. Utilisez un modèle de régression linéaire avec deux prédicteurs (x, y) et une variable de résultat (z).

python

from numpy import *
Nobs = 20
x_true = random.uniform(0,10, size=Nobs)
y_true = random.uniform(-1,1, size=Nobs)
alpha_true = 0.5
beta_x_true = 1.0
beta_y_true = 10.0
eps_true = 0.5
z_true = alpha_true + beta_x_true*x_true + beta_y_true*y_true
z_obs = z_true + random.normal(0, eps_true, size=Nobs)

python

%matplotlib inline
from matplotlib import pyplot as plt
plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.scatter(x_true, z_obs, c=y_true, marker='o')
plt.colorbar()
plt.xlabel('X')
plt.ylabel('Z')
plt.subplot(1,2,2)
plt.scatter(y_true, z_obs, c=x_true, marker='o')
plt.colorbar()
plt.xlabel('Y')
plt.ylabel('Z')

python

import pymc
# define the parameters with their associated priors
alpha = pymc.Uniform('alpha', -100,100, value=median(z_obs))
betax = pymc.Uniform('betax', -100,100, value=std(z_obs)/std(x_true))
betay = pymc.Uniform('betay', -100,100, value=std(z_obs)/std(y_true))
eps = pymc.Uniform('eps', 0, 100, value=0.01)

# Now define the model
@pymc.deterministic
def model(alpha=alpha, betax=betax, betay=betay, x=x_true, y=y_true):
    return alpha + betax*x + betay*y

# pymc parametrizes the width of the normal distribution by tau=1/sigma**2
@pymc.deterministic
def tau(eps=eps):
    return power(eps, -2)

# Lastly relate the model/parameters to the data
data = pymc.Normal('data', mu=model, tau=tau, value=z_obs, observed=True)

Une fois la définition en place, envoyez un objet pymc à pymc.MCMC pour créer et exécuter un échantillonneur. Nous appliquons ici l'échantillonneur Metropolis-Hastings en utilisant use_step_method. Il s'agit d'un échantillonneur plus intelligent que le Metropolis-Hastings par défaut. Pour échantillonner plus efficacement l'espace des paramètres, nous utiliserons l'échantillon précédent pour construire une matrice de covariance.

python

sampler = pymc.MCMC([alpha,betax,betay,eps,model,tau,z_obs,x_true,y_true])
sampler.use_step_method(pymc.AdaptiveMetropolis, [alpha,betax,betay,eps],
                        scales={alpha:0.1, betax:0.1, betay:1.0, eps:0.1})
sampler.sample(iter=10000)

[100%**] 10000 of 10000 complete

Terrain.

python

pymc.Matplot.plot(sampler)

Je vais coller uniquement l'alpha. Pour chaque probabilité, pymc trace les traces (panneau supérieur gauche), l'autocorrélation (panneau inférieur gauche) et les histogrammes d'échantillons (panneau droit). Vous pouvez voir que rien ne se passe au début. Après cela, l'échantillonneur adaptatif MH trouvera une meilleure matrice de covariance.

Ensuite, obtenons la valeur médiane de chaque variable de l'échantillon. L'échantillon MCMC examine la trace de chaque variable.

python

m_alpha = median(alpha.trace())
m_betax = median(betax.trace())
m_betay = median(betay.trace())
m_eps = median(eps.trace())

La médiane est plus robuste aux distributions asymétriques (par exemple, les distributions eps). Cette fois, nous modifierons le prédicteur y lors de l'inversion du tracé et de x. De cette façon, vous pouvez voir si vous avez trouvé la bonne tendance. Vous pouvez également tracer des zones ombrées qui représentent une variation intrinsèque.

python

plt.figure(figsize=(12,6))
plt.subplot(1,2,1)
plt.plot(x_true, z_obs-m_alpha-m_betay*y_true, 'o')
plt.xlabel('X')
plt.ylabel('Z - alpha - beta_y y')
# Now plot the model
xx = array([x_true.min(), x_true.max()])
plt.plot(xx, xx*m_betax)
plt.plot(xx, xx*m_betax + m_eps, '--', color='k')
plt.plot(xx, xx*m_betax - m_eps, '--', color='k')
plt.subplot(1,2,2)
plt.plot(y_true, z_obs-m_alpha-m_betax*x_true, 'o')
plt.xlabel('Y')
plt.ylabel('Z - alpha - beta_x x')
yy = array([y_true.min(), y_true.max()])
plt.plot(yy, yy*m_betay)
plt.plot(yy, yy*m_betay + m_eps, '--', color='k')
plt.plot(yy, yy*m_betay - m_eps, '--', color='k')

Examinez la covariance de divers paramètres dans le package de tracé triangle_plot. Pour ce faire, il semble qu'une chaîne de Markov de tableaux bidimensionnels indexés par [i, p] soit nécessaire. Concaténez tous les éléments en échantillons et transposez.

python

samples = array([alpha.trace(),betax.trace(),betay.trace(),eps.trace()]).T
print((samples.shape))
import triangle
triangle.corner(samples[:,:], labels=['alpha','betax','betay','eps'], truths=[alpha_true, beta_x_true, beta_y_true, eps_true])

Exemple d'utilisation

https://github.com/QuantScientist/kaggle-seizure-prediction-challenge-2016/blob/master/jupyter/bayesian-logistic-regression/ieegPymc3Version1.ipynb

référence

https://users.obs.carnegiescience.edu/cburns/ipynbs/PyMC.html https://pymc-devs.github.io/pymc/INSTALL.html#dependencies