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:
@pm.deterministic def coupling(phi=phi, a=a,b=b,c=c): """ Calculate the J-coupling from Karplus forumula """ ang=phi*3.1415/180. 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.