This is my first post and I am a newby in pymc. I am trying to model a non-linear system (see below for a further explanation). I create my synthetic data with:

data = np.random.normal(8., data_std, 1000)
obs = pm.Normal("obs", coupling, value=data, observed=True)

My stochastic variable is:

phi = pm.Uniform("phi", -180, 180., value=-150)

and it is related to the data though the deterministic function:

def coupling(phi=phi, a=a,b=b,c=c):
     Calculate the J-coupling  from Karplus forumula
     return a*np.cos(ang)**2+b*np.cos(ang)+c

(a,b,c are scalars. In the future, I'll try to see what happens when they are also stochatic variables)

Then my model is:

model = pm.Model([obs, phi])
mcmc = pm.MCMC(model)
mcmc.sample(20000, burn=1000, burn_till_tuned=False)

What I don't understand is why the posterior distribution of phi does not depend on the dispersion of the data. That is, I get the same std of phi_trace whether data_std is 0.1 or 2. What I am doing wrong?

phi_trace = mcmc.trace("phi")[:]
print "\n Phi std:", phi_trace.std()

Thanks in advance, Ramon

PS. I know the nonlinearity of the system may cause the model to get trapped to one of the solutions, say close to -150 and not explore 0 or +150. I'll try to deal with that in the future...

PS2. Some background: phi is a dihedral angle of a molecule and the coupling is what can be measured experimentally. In fact a uniform distribution for phi is a bad prior, but I wanted to see the difference when using more informative priors.


I am answering my own question... The reason is that I am defining the observed data (obs) as a normal distribution, but I am not including the precision of the distribution as a stochastic variable. The correct model should be:

tau = pm.Normal('tau', 5., .05) #or other possible choices...
obs = pm.Normal("obs", coupling, tau, value=data, observed=True)
model = pm.Model([obs, phi, tau])

Now the posterior of phi depends on the dispersion of data, as expected...


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.