[pymc] translating simple model

[bblais]

Hello,

I just started with PyMC, and am running into a few problems. Let me
describe the model, and then show the code I thought should work, :),
and then tell you the error. I'm trying to predict the win/loss for a
particular game, and I'm given measurable variable (which I call beta)
which is a ratio of quantities from the two teams playing. My data,
therefore, is a series of values of beta, and whether the home team
won, so for example I'd have:

Data:
beta=[1.18, 0.8, 0.47]
xi= [ 1 , 1 , 0]

(obviously I'll have more data, but I wanted to make sure it worked
first). In the model, the probability for a win depends on beta, and
a parameter r which I want to eventually fit or look at the posterior
for. Thus my model is:

beta ~ Uniform(0,10) (observed)
r ~ Uniform(0,2) (not observed)
lambda_s = beta **r (definition)
ps=lambda_s/(1+lambda_s) (also a definition)
xi ~ Bernoulli(ps) (observed)

I defined my model for PyMC as model1.py:

# simple model, with only beta
from pymc import Uniform,deterministic,Bernoulli

beta=Uniform('beta',0,10,observed=True,value=[1.18,.8,.47])
r=Uniform('r',0,2)

@deterministic
def lambda_s(beta=beta,r=r):
return beta**r

@deterministic
def ps(lambda_s=lambda_s):
return lambda_s/(1+lambda_s)

xi=Bernoulli('xi',p=ps,observed=True,value=[1,1,0])

I tried to run my model with:

from pylab import *
import pymc
import model1


model=pymc.MCMC(model1)
model.sample(iter=10000, burn=1000,thin=2)

r=model.trace('r')[:]

hist(r,30)

draw()
show()

Even at the import model1 I get the error:

error: (len(p)==1 || len(p)==len(x)) failed for hidden np


I must be overlooking something basic here, I think, but does anyone
have an idea what?

thanks!

Brian Blais

--
You received this message because you are subscribed to the Google Groups "PyMC" group.
To post to this group, send email to py...@googlegroups.com.
To unsubscribe from this group, send email to pymc+uns...@googlegroups.com.
For more options, visit this group at http://groups.google.com/group/pymc?hl=en.

[Chris Fonnesbeck]

Hi Brian,

Your model runs fine as written for me (saved to a file called
blais.py):


In [1]: import blais

In [2]: from pymc import *

In [3]: M = MCMC(blais)

In [4]: M.isample(4000)
==============
PyMC console
==============

PyMC is now sampling. Use the following commands to query or
pause the sampler.


Commands:
i -- index: print current iteration index
p -- pause: interrupt sampling and return to the main
console.
Sampling can be resumed later with icontinue().
h -- halt: stop sampling and truncate trace. Sampling
cannot be
resumed for this chain.
b -- bg: return to the main console. The sampling will
still
run in a background thread. There is a
possibility of
malfunction if you interfere with the Sampler's
state or the database during sampling. Use this
at your
own risk.

pymc > Sampling terminated successfully.

In [6]: M.stats()
Out[6]:
{'lambda_s': {'95% HPD interval': array([[ 1.01501719, 1.38370687],
[ 0.64007366, 0.97028501],
[ 0.22098603, 0.90296912]]),
'mc error': array([ 0.00428021, 0.00384332,
0.00788886]),
'mean': array([ 1.19913781, 0.7938924 , 0.48870203]),
'n': 4000,
'quantiles': {2.5: array([ 1.01178093, 0.6465784 ,
0.22867721]),
25: array([ 1.10456662, 0.70447092,
0.30565643]),
50: array([ 1.19809621, 0.78375248,
0.43847681]),
75: array([ 1.29685278, 0.87475432,
0.63586855]),
97.5: array([ 1.38187827, 0.98483738,
0.94961679])},
'standard deviation': array([ 0.11205486, 0.10146593,
0.21235839])},
'ps': {'95% HPD interval': array([[ 0.50559349, 0.58198893],
[ 0.39027129, 0.49245922],
[ 0.18098981, 0.4745054 ]]),
'mc error': array([ 0.00088822, 0.00118887, 0.00345084]),
'mean': array([ 0.54408889, 0.44078865, 0.31536188]),
'n': 4000,
'quantiles': {2.5: array([ 0.50292798, 0.39267999,
0.18611659]),
25: array([ 0.52484279, 0.41330768,
0.23410173]),
50: array([ 0.54506086, 0.4393841 ,
0.30482022]),
75: array([ 0.56462164, 0.46659678,
0.38870394]),
97.5: array([ 0.58016326, 0.49618039,
0.48707869])},
'standard deviation': array([ 0.02332925, 0.03123384,
0.09125281])},
'r': {'95% HPD interval': array([ 0.13518393, 1.99948427]),
'mc error': 0.021657817560558582,
'mean': 1.0705769064596478,
'n': 4000,
'quantiles': {2.5: 0.070761638926591944,
25: 0.60087240984440038,
50: 1.0919518991428174,
75: 1.5704998073209293,
97.5: 1.9541717610812714},
'standard deviation': 0.56866845960872192}}

What version of pymc are you running?

[Brian Blais]

On Apr 19, 2010, at 18:44 , Chris Fonnesbeck wrote:

Hi Brian,

Your model runs fine as written for me (saved to a file called

blais.py):

What version of pymc are you running?

Thanks!  I was running 2.0, which was installed via easy_install on my OSX machine.  Now I compiled the svn version (2.1.alpha), and even though it gives weird linking errors (mostly about arch ppc flags) it does seem to work anyway. 

Glad I wasn't missing anything obvious!

bb

--

Brian Blais

bbl...@gmail.com

http://web.bryant.edu/~bblais

http://bblais.blogspot.com/