Hi,
I find the useful small code project for me:
#https://users.obs.carnegiescience.edu/cburns/ipynbs/PyMC.html
It runs as expected.
When I review the code, I find 'data' in the original line:
data = pymc.Normal('data', mu=model, tau=tau, value=z_obs, observed=True)
has not been referenced thereafter.
If I comment out the line as:
#data = pymc.Normal('data', mu=model, tau=tau, value=z_obs, observed=True)
the result is ugly different from the original.
If I change it to:
pymc.Normal('data', mu=model, tau=tau, value=z_obs, observed=True)
it still runs as the original. Reading the last half part
(after data= line), I cannot see anything uses 'data', though I suspect
below line may use it:
sampler = pymc.MCMC([alpha,betax,betay,eps,model,tau,z_obs,x_true,y_true])
This is quite different from my any other language experience.
Could you help me on this?
Thanks,
--------------
#https://users.obs.carnegiescience.edu/cburns/ipynbs/PyMC.html
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)
import pymc
# define the parameters with their associated priors
#alpha = pymc.Normal('alpha',mu=0,tau=.01)
#beta = pymc.Normal('beta',mu=0,tau=.01)
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)
pymc.Normal('data0', mu=model, tau=tau, value=z_obs, observed=True)
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)
pymc.Matplot.plot(sampler)
sampler.sample(iter=10000)
alpha.summary()
m_alpha = median(alpha.trace())
m_betax = median(betax.trace())
m_betay = median(betay.trace())
m_eps = median(eps.trace())
import matplotlib.pyplot as plt
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')
plt.show()
--
https://mail.python.org/mailman/listinfo/python-list
