5
$\begingroup$

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.

enter image description here

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:

enter image description here

EDIT 2, 15K samples, peaking at 9.93

enter image description here

Am I just being persnickety, or is there something wrong?

Improve this question
$\endgroup$
  • 1
    $\begingroup$ Maybe the process didnt converge - try increasing number of samples(to more than 2000). $\endgroup$ – yoav_aaa Nov 13 '18 at 12:08
  • 1
    $\begingroup$ @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. $\endgroup$ – Demetri Pananos Nov 13 '18 at 14:07
  • $\begingroup$ OP, what is your expected result? Do you expect the posterior mode to be 10? $\endgroup$ – Demetri Pananos Nov 13 '18 at 14:16
  • $\begingroup$ 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. $\endgroup$ – RNs_Ghost Nov 13 '18 at 14:19
2
$\begingroup$

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 = 2, 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.

Improve this answer
$\endgroup$
  • $\begingroup$ 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) ? $\endgroup$ – RNs_Ghost Nov 13 '18 at 14:48
  • 1
    $\begingroup$ 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]. $\endgroup$ – Demetri Pananos Nov 13 '18 at 15:07

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.