Multinomial, Lambda and reshape

[Hajas, Wayne]

It took me a while to figure this out so I may as well document it here.  The problem is solved in http://idisk.me.com/fonnesbeck-Public/Mt.py but I didn't immediatley figure out the details.

I want to use the multinomial distribution.

My first attempt was

Mort=Uniform('Mort',.01,.99,.1)
nAge=10
P=[]
for t in range(-0+nAge):
    @deterministic
    def ProbAge(Mort=Mort, nAge=nAge, Age=t):
        return (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-t))
    ProbAge.__name__="P["+str(t)+"]"
    P.append(ProbAge)
   AF=Multinomial('AF',n=nAnimal,p=P,value=AgeFreq,observed=True)

There were no error-messages.  When I ran the model, the parameter values remained constant.  Something was wrong.

Then I found the Lambda function and tried this.

Mort=Uniform('Mort',.01,.99,.1)
nAge=10
Age=range(nAge)
P = Lambda('P', lambda Mort=Mort: (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-Age)))

The error I got was

NameError: global name 'exp' is not defined

Then I used reshape for my numerical constants.

 Mort=Uniform('Mort',.01,.99,.1)
 nAge=10
 Age=reshape(range(nAge),nAge)
 P = Lambda('P', lambda Mort=Mort: (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-Age)))

And everything seems to work OK.

Wayne Hajas.

[anand.prabhakar.patil]

Hi Wayne,

Thanks for persevering, glad you got it working. Comments on why your first attempts may not have worked follow.

On Wed, Jan 21, 2009 at 11:13 PM, Hajas, Wayne <Wayne...@dfo-mpo.gc.ca> wrote:

My first attempt was

Mort=Uniform('Mort',.01,.99,.1)
nAge=10
P=[]
for t in range(-0+nAge):
    @deterministic
    def ProbAge(Mort=Mort, nAge=nAge, Age=t):
        return (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-t))
    ProbAge.__name__="P["+str(t)+"]"
    P.append(ProbAge)
   AF=Multinomial('AF',n=nAnimal,p=P,value=AgeFreq,observed=True)

There were no error-messages.  When I ran the model, the parameter values remained constant.  Something was wrong.

I'm assuming the multinomial was originally outside the loop and tried the following:

from pymc import *
from numpy import *

Mort=Uniform('Mort',.01,.99,.1)
nAge=10
nAnimal = 15

P=[]
for t in range(-0+nAge):
    @deterministic
    def ProbAge(Mort=Mort, nAge=nAge, Age=t):
        return (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-t))
    ProbAge.__name__="P["+str(t)+"]"
    P.append(ProbAge)

AgeFreq = rmultinomial(nAnimal, Container(P).value)
AF=Multinomial('AF',n=nAnimal,p=P,value=AgeFreq,observed=True)

As you say, Mort 'flatlines', meaning the MCMC is never able to accept a jump.

You've been bitten by one of the most subtle, confusing behaviors of Python, which unfortunately causes bugs in PyMC programs that use the decorator syntax relatively frequently. It has to do with Python's 'lexical scoping'. Rather than giving a confusing explanation, I'll point out where the problem is and how to fix it, and leave you to experiment and get used to how it works.

The problem is in the declaration of the deterministic:

    @deterministic
    def ProbAge(Mort=Mort, nAge=nAge, Age=t):

        return (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge- *** t ***))

I've put stars around the problem, which is that you've used 't' in the function body rather than Age. Replace 't' with 'Age' and everything works fine. As a general rule, it's critical to be disciplined about declaring every variable in the body of a deterministic as a parent to avoid such issues.

Then I found the Lambda function and tried this.

Mort=Uniform('Mort',.01,.99,.1)
nAge=10
Age=range(nAge)
P = Lambda('P', lambda Mort=Mort: (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-Age)))

The error I got was

NameError: global name 'exp' is not defined

That probably means a 'from numpy import *' line got commented out when you were trying this iteration. The 'exp' function isn't part of core Python by default, meaning if you just start up a shell and type 'exp' you'll see the same thing. Exp has to be imported from numpy before you can use it.

Incidentally, you'll still get an error here because your deterministic P, when it does  (-1+nAge-Age), is trying to add the list Age to the number nAge. Lists can't be added to numbers in python, so Age has to be converted to a proper numpy array in order for this to work...

Then I used reshape for my numerical constants.

 Mort=Uniform('Mort',.01,.99,.1)
 nAge=10
 Age=reshape(range(nAge),nAge)
 P = Lambda('P', lambda Mort=Mort: (1-exp(Mort)) / (1-exp(nAge*Mort)) * exp(Mort*(-1+nAge-Age)))

which is exactly what this does. Alternative ways to get the same effect would be:

Age = arange(nAge)
Age = asarray(range(nAge))
Age = array(range(nAge))

Anand