Querying different results from simple model in PyMC vs BUGS

[Steve]

Hi,

This is my first attempt at using PyMC. I have translated some WinBUGS
code for a simple population ecology density feedback model into PyMC.
The model runs successfully, however the resulting posterior
distributions for parameters 'K' and 'sigma' are highly negatiely
biased relative to the results i get from BUGS (and the MLE's), and
the precision of all estimated parameters is far too high.

I have obviously messed up the sampling in some way but have been
unable to see the problem. If someone were prepared to look over the
PyMC and WinBUGS code (below) and point me in the right direction it
would be very much appreciated.

Thanks for your time.
Cheers,
Steve

# The model
N(i+1) = N(i) * exp(rm * (1 - N[i]/K) + ev[i])
where ev[i] ~ N(0, sigma**2)

# WinBUGS code
"""
model {
for(i in 1:lngTS) { # for each year
# Environmental stochasticity in population growth rate (drawn
from gaussian)
ev[i] ~ dnorm(0, tau)
# Per capita growth rate is density dependent - Ricker-type
Er[i] <- exp(rm * (1 - N[i]/K) + ev[i])
# Pmean equal to number this year times per capita rate
Pmean[i] <- N[i] * Er[i]
# Number next year drawn from Poisson with mean Pmean -
demographic stochasticity
N[i + 1] ~ dpois(Pmean[i])
}
#PRIORS
# rm is maximum exponential growth rate
rm ~ dunif(0, 2)
# K is carrying capacity
K ~ dunif(0, max(N))
# SD in growth rate due to environmental stochasticity
sigma ~ dunif(0, 1)
tau <- 1/(sigma*sigma)
}
"""

######################################################################################

## Ricker density feedback model. Original WinBUGS model code above.

######################################################################################

#------------------------------------------------------------------
MODULES
from numpy import *
from pymc import *
import pdb

#------------------------------------------------------------------
DATA
# Number of individuals at time i
Ni = [354, 356, 426, 405, 628, 618, 603, 465, 299, 334, 321, 385, 348,
373, 598, 592, 425, 539, 456, 360]

# Number of observations -1
lngTS = len(Ni) - 1

# Initial value for enviromental stochasticity
ev_inits = ones(lngTS) * 0.1

#------------------------------------------------------------------
PRIORS
# rm is maximum exponential growth rate
rm = Uniform('rm', lower = 1e-8, upper = 2, value = 0.2)

# K is carrying capacity
K = Uniform('K', lower = 1e-8, upper = max(Ni), value = 500.0)

# sigma is st. dev. in growth rate due to env. stoch.
sigma = Uniform('sigma', lower = 1e-8, upper = 2.0, value = 0.2)

# Convert sigma to precision scale
itau = Lambda('itau', lambda sigma=sigma: 1/(sigma**2.0))

#------------------------------------------------------------------
MODEL
Pmean = []
# Recursive step
for i in range(lngTS):
# Environmental stochasticity in population growth rate (drawn
from normal)
evi = Normal('evi', 0, itau, value = ev_inits[i])
# Per capita growth rate is density dependent - Ricker
Eri = Lambda("Eri", lambda Ni = Ni[i], evi = evi, rm = rm, K = K:
exp(rm * (1 - (Ni/K)) + evi))
# Pmean equal to number this year times per capita rate
Pmeani = Lambda("Pmean_%i"%i, lambda Ni = Ni[i], Eri = Eri: Ni *
Eri)
Pmean.append(Pmeani)

# Demographic stochasticity - number next year drawn from Poisson with
mean 'Pmean'
# (i.e. est. value for shape of Poisson distribution given number of
individuals observed)
@stochastic(observed = True)
def y(value = Ni[1:], mu = Pmean):
return poisson_like(value, mu)

#------------------------------------------------------------------
RUN
if __name__ == '__main__':

var_list = [rm, K, sigma, y]

# Create MCMC object
M = MCMC(var_list)

# Sample
M.isample(20000, burn = 2000, thin = 25, verbose = 2)

M.stats()

[Chris Fonnesbeck]

On Sep 29, 9:59 pm, Steve <jsdel...@gmail.com> wrote:
> Hi,
>
> This is my first attempt at using PyMC. I have translated some WinBUGS
> code for a simple population ecology density feedback model into PyMC.
> The model runs successfully, however the resulting posterior
> distributions for parameters 'K' and 'sigma' are highly negatiely
> biased relative to the results i get from BUGS (and the MLE's), and
> the precision of all estimated parameters is far too high.

Hi Steve,

Let me have a look later in the day. I have ported a few models to
PyMC from BUGS (see our homepage) and have gotten good results, so
perhaps there is some nuance in your model that is causing the
difference.

cf

[Chris Fonnesbeck]

On Sep 30, 9:19 am, Chris Fonnesbeck <fonnesb...@maths.otago.ac.nz>
wrote:

> On Sep 29, 9:59 pm, Steve <jsdel...@gmail.com> wrote:
>
> > Hi,
>
> > This is my first attempt at using PyMC. I have translated some WinBUGS
> > code for a simple population ecology density feedback model into PyMC.
> > The model runs successfully, however the resulting posterior
> > distributions for parameters 'K' and 'sigma' are highly negatiely
> > biased relative to the results i get from BUGS (and the MLE's), and
> > the precision of all estimated parameters is far too high.

Steve,

The main issue with your code is that you were not dealing with the
elements of ev and Er correctly. They should be stashed away in lists,
just as you appropriately did with Pmean. Its a subtle and confusing
point, to be sure, but it will certainly mess up your results. This is
one reason that I like to handle things in vectors rather than in
loops with PyMC, though you had little choice in an explicitly
Markovian model like this. Here is what your loop should look like:

# MODEL
Pmean = []
ev = []
Er = []

# Recursive step
for i in range(lngTS):
# Environmental stochasticity in population growth rate (drawn from
normal)

evi = Normal("evi_%i"%i, 0, itau, value = ev_inits[i])
ev.append(evi)

# Per capita growth rate is density dependent - Ricker

Eri = Lambda("Eri_%i"%i, lambda Ni = Ni[i], evi = evi, rm = rm, K =

K: exp(rm * (1 - (Ni/K)) + evi))

Er.append(Eri)

# Pmean equal to number this year times per capita rate
Pmeani = Lambda("Pmean_%i"%i, lambda Ni = Ni[i], Eri = Eri: Ni *

Eri, plot=False)
Pmean.append(Pmeani)

Here are posterior distributions for sigma, K and rm, which compare
nicely to what I get in WinBUGS:

http://bit.ly/3CPJOA
http://bit.ly/2hnGkc
http://bit.ly/zDciM

One other thing -- your list of variables in var_list did not include
your ev stochastics. I really recommend putting your specified model
into one file and control that model with PyMC from a separate file,
say "runmodel.py":

import ricker
from pymc import MCMC, Matplot

M = MCMC(ricker)
M.sample(10000, 9000, verbose=2)
Matplot.plot(M)


Hope that clears things up,
cf

[Steve]

Hi Chris,

Thank you very much for your help. I see the differences and will now
spend some more time on the documentation to try to better understand
them. I have also taken your advice and am using separate .py files
for specifying and running the model.

A follow up question if I may, whilst the code now yields comparable
results with those from WinBUGS, the computation time using PyMC is
about 40 times greater than for WinBUGS (run through R2WinBUGS with
same #iterations/burnin/thinning). I am guessing this is a consequence
of having to use a loop to calculate individual terms of the objective
function. Did you also note the substantial difference in computation
time when you did the comparison? My goal was to use PyMC to fit such
models to many simulated data sets but without significant reduction
in computation time this may not be possible. Do you have any thoughts
on the speed differences?

Thanks again for your time and for this excellent resource.
Regards,
Steve

On Oct 2, 11:58 am, Chris Fonnesbeck <fonnesb...@maths.otago.ac.nz>

> http://bit.ly/3CPJOAhttp://bit.ly/2hnGkchttp://bit.ly/zDciM

[Chris Fonnesbeck]

On Oct 2, 9:07 pm, Steve <jsdel...@gmail.com> wrote:
> A follow up question if I may, whilst the code now yields comparable
> results with those from WinBUGS, the computation time using PyMC is
> about 40 times greater than for WinBUGS (run through R2WinBUGS with
> same #iterations/burnin/thinning). I am guessing this is a consequence
> of having to use a loop to calculate individual terms of the objective
> function. Did you also note the substantial difference in computation
> time when you did the comparison? My goal was to use PyMC to fit such
> models to many simulated data sets but without significant reduction
> in computation time this may not be possible. Do you have any thoughts
> on the speed differences?

Steve,

Looking at your model again, it appears that it can be re-expressed to
remove the costly loops -- that is what is slowing it down. Have a
look at what I came up with below. It runs much faster on my machine
(Mac Powerbook):

from numpy import *
from pymc import *
import pdb

#------------------------------------------------------------------

# DATA

# Number of individuals at time i

N = [354, 356, 426, 405, 628, 618, 603, 465, 299, 334, 321, 385, 348,

373, 598, 592, 425, 539, 456, 360]

# Number of observations -1

lngTS = len(N) - 1

#------------------------------------------------------------------
# MODEL

# rm is maximum exponential growth rate
rm = Uniform('rm', lower = 1e-8, upper = 2, value = 0.2)

# K is carrying capacity

K = Uniform('K', lower = 1e-8, upper = max(N), value = 500.0)

# sigma is st. dev. in growth rate due to env. stoch.

sigma = Uniform('sigma', lower = 1e-8, upper = 1.0, value = 0.2)

# Convert sigma to precision scale
itau = Lambda('itau', lambda sigma=sigma: 1/(sigma**2.0))

# Environmental stochasticity in population growth rate (drawn from
normal)

ev = Normal('ev', mu=0, tau=itau, value=ones(lngTS) * 0.1)

# Per capita growth rate is density dependent - Ricker

Er = Lambda("Eri", lambda N=N[:-1], ev=ev, rm=rm, K=K: exp(rm * (1 -
(N/K)) + ev))

# Pmean equal to number this year times per capita rate

Pmean = Lambda('Pmean', lambda N=N, Er=Er: [n*er for n,er in zip
(N,Er)])

# Demographic stochasticity - number next year drawn from Poisson with
mean 'Pmean'
# (i.e. est. value for shape of Poisson distribution given number of
individuals observed)

y = Poisson('y', mu=Pmean, value=N[1:], observed=True)

[Steve]

Hi Chris,

Thank you very much, this is fantastic. Your help is much appreciated,
the computation time is now very similar to that under BUGS. I am
pleased to say that i had made some progress towards a similar
solution, but was using decorator classes and could not quite get the
right inputs for the @stochastic variable.

I did have to make one change to the code you sent - i had to coerce N
[:-1] to an array before it would execute in the line
Er = Lambda("Eri", lambda N=array(N[:-1]), ev=ev, rm=rm, K=K: exp(rm *

(1 - (N/K)) + ev))

The error message given (without making the above change) included
reference to lazy function evaluation (error history appended at end
of message). So I was wondering whether this is a consequence of my
install failing the LazyFunction test as:

FAIL: test (pymc.tests.test_LazyFunction.test_LazyFunction)
----------------------------------------------------------------------
Traceback (most recent call last):
File "/usr/lib/python2.5/site-packages/pymc/tests/
test_LazyFunction.py", line 55, in test
assert(L.ultimate_args.value[1] is C.value)

Many thanks for your time once again,
Regards,
Steve

## Error issued from original code sent
In [3]: ERROR: An unexpected error occurred while tokenizing input
The following traceback may be corrupted or invalid
The error message is: ('EOF in multi-line statement', (34, 0))

---------------------------------------------------------------------------
TypeError Traceback (most recent call
last)

/home/sduni/Python/PyMC/examples/<ipython console> in <module>()

/home/sduni/Python/PyMC/examples/cf_rickerMod.py in <module>()
32
33 # Per capita growth rate is density dependent - Ricker
---> 34 Er = Lambda("Eri", lambda N=N[:-1], ev=ev, rm=rm, K=K: exp(rm

* (1 - (N/K)) + ev))

35
36 # Pmean equal to number this year times per capita rate

/usr/lib/python2.5/site-packages/pymc/CommonDeterministics.pyc in
__init__(self, name, lam_fun, doc, *args, **kwds)
64 parents = dict(zip(parent_names[-len
(parent_values):], parent_values))
65
---> 66 pm.Deterministic.__init__(self, eval=lam_fun,
name=name, parents=parents, doc=doc, *args, **kwds)
67
68 def lambda_deterministic(*args, **kwargs):

/usr/lib/python2.5/site-packages/pymc/PyMCObjects.pyc in __init__
(self, eval, doc, name, parents, dtype, trace, cache_depth, plot,
verbose)
341 trace=trace,
342 plot=plot,
--> 343 verbose=verbose)
344
345 # self._value.force_compute()

/usr/lib/python2.5/site-packages/pymc/Node.pyc in __init__(self, doc,
name, parents, cache_depth, trace, dtype, plot, verbose)
148 self.extended_children = set()
149
--> 150 Node.__init__(self, doc, name, parents, cache_depth,
verbose=verbose)
151
152 if self.dtype is None:

/usr/lib/python2.5/site-packages/pymc/Node.pyc in __init__(self, doc,
name, parents, cache_depth, verbose)
67
68 # Initialize
---> 69 self.parents = parents
70
71 # New lazy function

/usr/lib/python2.5/site-packages/pymc/Node.pyc in _set_parents(self,
new_parents)
90
91 # Get new lazy function
---> 92 self.gen_lazy_function()
93
94 parents = property(_get_parents, _set_parents, doc="Self's
parents: the variables referred to in self's declaration.")

/usr/lib/python2.5/site-packages/pymc/PyMCObjects.pyc in
gen_lazy_function(self)
352 cache_depth =
self._cache_depth)
353
--> 354 self._value.force_compute()
355
356 def get_value(self):

/usr/lib/python2.5/site-packages/pymc/LazyFunction.so in
pymc.LazyFunction.LazyFunction.force_compute()

/home/sduni/Python/PyMC/examples/cf_rickerMod.py in <lambda>(N, ev,
rm, K)
32
33 # Per capita growth rate is density dependent - Ricker
---> 34 Er = Lambda("Eri", lambda N=N[:-1], ev=ev, rm=rm, K=K: exp(rm

* (1 - (N/K)) + ev))

35
36 # Pmean equal to number this year times per capita rate

TypeError: unsupported operand type(s) for /: 'list' and 'float'

AssertionError



# The original code was: