Help on code comprehension from an example project of pymc

[Robert]

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()