# Finding the Poisson rate parameter with PyMC3

I'm trying to compute the rate parameter of fake set of poisson data, where I set the parameter.

When I run PyMC the posterior distribution always peaks around the true rate parameter, but never seems to hit it.

For example, in this scenario, the true rate parameter is 10 but the posterior peaks at 10.05. What have I done wrong? Here's the code:

true_rate = 10
n_obs = 10,000
poisson_accidents = np.random.poisson(true_rate, 10000)

# fit a pymc3 model to get the rate
basic_model = pm.Model()

with basic_model:
lambda_prior = pm.Uniform('lambda_prior',0,20)

#likelihood
Y_obs = pm.Poisson('Y_obs',mu = lambda_prior, observed = poisson_accidents)

with basic_model:
trace = pm.sample(2000,cores = 1,chains = 10,tune = 5000)
pm.traceplot(trace)
plt.show()


EDIT

I have increased the number of iterations as suggested:

with basic_model:
trace = pm.sample(10000,cores = 1,chains = 10,tune = 5000)
pm.traceplot(trace)
plt.show()


This still results in a peak that is not centered at 10: EDIT 2, 15K samples, peaking at 9.93 Am I just being persnickety, or is there something wrong?

• Maybe the process didnt converge - try increasing number of samples(to more than 2000). Nov 13, 2018 at 12:08
• @DaFanat The hamiltonian Monte Carlo methods that are used my PyMC3 usually converge in 2000 interations. The problem is relatively simple, so I doubt this is the problem. I think OP has his standards too high. Nov 13, 2018 at 14:07
• OP, what is your expected result? Do you expect the posterior mode to be 10? Nov 13, 2018 at 14:16
• I suppose I just wanted to check that I hadn't committed some grave mcmc sin. As long as you can't see anything obvious. If you want to put your comment in an answer, I'll mark it as correct. Nov 13, 2018 at 14:19

First of all, you have far too many chains for a problem like this. Second of all, your tuning parameter is far too high. Something along the lines of

trace = pm.sample(2000,chains = 4, tune = 1000)


Should be sufficient for a problem like this.

Remember that there is sampling error in the data. The mean of the generated data is rarely ever exactly 10. Since you are using a uniform prior, then the posterior mode is the maximum likelihood estimator (which for the poisson is the mean of the observations, if I recall correctly). There might be a little bias because you are using a prior which is not supported on the parameter's domain, but with this many observations, I would expect the bias (if it exists) to be small.

I see nothing wrong here, the behaviour is exactly as to be expected. HMC is not omniscient, and so expecting the posterior mode to be exactly the true parameter value is holding a very high standard. Instead, I would plot HPD intervals (with pm.plot_posterior(trace)) and determine if the true parameter lies within the interval.

• Will having "too many chains" and too many burn in samples reduce the quality of results, or is it just unnecessary? Also, what do you mean by the prior not being in the parameter's domain? 10 is within normal(0,20) ? Nov 13, 2018 at 14:48
• Too many chains and too many warmup iterations is just not needed. Makes the computation take longer. The parameter for the poisson is supported on the positive real line, where as you have passed a prior supported on [0,20]. Nov 13, 2018 at 15:07