pymc normal prior + normal likelihood does not converge correctly?

540 views Asked by At

I am new to pymc and bayesian statistics. Here I am trying to implement an extremely simple pymc model in order to compare with the theoretical result. In my test case, I assumed a normal prior as mu~N(20,20) and a likelihood assumed to be data~N(mu,10)

Let there are 10 observations with values as shown in the code below, in this simple model, the theoretical result of the posterior pdf should be N(6.84,6.67); however, I am confused why my pymc model produce the results far from the theoretical results. Here is my code. I am wondering if the problem from my code or my concept of bayesian statistics. Thanks for your help.

from __future__ import division
import pymc as pm
import numpy as np

data=[2.944,-13.361,7.143,16.235,-6.917, 8.580,12.540,
          -15.937,-14.409, 5.711]
mean=pm.Normal('mean',mu=20.,tau=1./20)
prec=1./10

obs=pm.Normal('obs',mu=mean,tau=prec,value=data,observed=True)
model=pm.Model([obs, mean])
mcmc=pm.MCMC(model)

mcmc.sample(500,100)

mcmc.stats()

The results shows the posterior mean is N(0.25,1) which is far away from the theoretical results.

'mean': {'95% HPD interval': array([-1.80515483,  2.06741224]),
  'mc error': 0.04850118114229577,
  'mean': 0.22188919237458093,
  'n': 400,
  'quantiles': {2.5: -1.7106911626432717,
   25: -0.39834886222214749,
   50: 0.24108945921296354,
   75: 0.85983578287420315,
   97.5: 2.282198749772455},
  'standard deviation': 0.99310330871482888}
1

There are 1 answers

0
Enson Portela On

Should mean=pm.Normal('mean',mu=20.,tau=1./20**2)?