Groups | Search | Server Info | Keyboard shortcuts | Login | Register [http] [https] [nntp] [nntps]

Groups > comp.lang.python > #100469 > unrolled thread

Help on code comprehension from an example project of pymc

Back to article view | Back to comp.lang.python

Contents

  Help on code comprehension from an example project of pymc Robert <rxjwg98@gmail.com> - 2015-12-15 08:15 -0800
    Re: Help on code comprehension from an example project of pymc Chris Angelico <rosuav@gmail.com> - 2015-12-16 09:34 +1100
    Re: Help on code comprehension from an example project of pymc Terry Reedy <tjreedy@udel.edu> - 2015-12-15 17:45 -0500

#100469 — Help on code comprehension from an example project of pymc

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

[toc] | [next] | [standalone]

#100484

On Wed, Dec 16, 2015 at 3:15 AM, Robert <rxjwg98@gmail.com> wrote:
> 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.
>
> This is quite different from my any other language experience.
>
> Could you help me on this?

What this suggests is that the call you're doing returns something (as
every function call must, unless it raises an exception or doesn't
terminate), but it also has side effects. You're ignoring the return
value (which is why "data = " doesn't affect your code), but you need
the side effects. Unfortunately this is the case with quite a lot of
Python's high end math/stats modules, presumably because their APIs
are lifted directly from other languages; have a look at pyplot, for
instance, where a more Pythonic API would be to instantiate a plot
object every time. (In that particular example, I believe that's an
option; but most examples just use the module-level default.)

It's my guess (based on skimming the docs) that the objects passed in,
including those referenced as function default arguments, are getting
mutated. But I'm no expert on pymc. It's equally likely that
module-level (global) state is being changed.

ChrisA

[toc] | [prev] | [next] | [standalone]

#100485

On 12/15/2015 11:15 AM, Robert wrote:
> 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 the function is called strictly for its side-effect, then it would be 
normal to not keep the 'return' value.  Code checkers will catch this 
and warn.  Just because code is make available, does not mean it follows 
the best style.  Perhaps the programmer though 'data' might be needed 
before writing the rest.

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



-- 
Terry Jan Reedy

[toc] | [prev] | [standalone]

Back to top | Article view | comp.lang.python

csiph-web