Credible interval for gamma prior

Question

Let's cosnider $X_1, X_2,...,X_n \sim \textrm{Poisson}(\lambda)$ i.i.d, and $\textrm{Gamma}(\alpha, \beta)$ as prior. Then posterior distribuion is $\textrm{Gamma}(\sum_i{X_i} + \alpha, n + \beta)$. We can deduce that bayesian estimator will be given as: $$\hat T = \frac{\sum X_i + \alpha}{n + \beta}$$

I want to ask you for help in finidng $(1-c)$ credible interval for this problem. I've read that credible interval is such interval $(l, r)$ that satisfies: $$\int_l^r p(\theta \mid x) = 1 - c$$

But I'm not sure how to derive this with it. Do you know how this $l, r$ can be derived?

Answer 1

If you find $r$ such that $$\int_0^r p(\theta | x) = 1-c/2$$ and $l$ such that $$\int_0^l p(\theta | x) = c/2$$ Then $(l,r)$ is a credible interval because $$\int_l^r p(\theta | x) = 1-c$$

The value of $r$ satisfying the above is, by definition, the quantile function (inverse of the CDF) evaluated at $1-c/2$. The quantile function of the gamma distribution does not have a convenient closed form but can be approximated. In R, you can use qgamma.

You can also obtain a credible interval by sampling from the posterior distribution. I.e., take $10000$ samples from a $\text{Gamma}(X_i + \alpha, n + \beta)$ and look at the empirical quantiles. This is oftentimes the only available option when the models get more complicated...

Answer 2

Hint, using a similar problem: Suppose the gamma prior distribution has shape parameter $4$ and rate parameter $1/3$ and that the sum of $50$ independent observations from $\mathsf{Pois}(\lambda)$ is $256.$ The sample mean is $5.12.$

Then, by Bayes' Theorem, the posterior distribution is gamma with shape parameter $260$ and rate parameter $50.333,$ so that a 95% Bayesian credible interval estimate for $\lambda$ is approximately $(4.557,\, 5.812).$

In R, where qgamma is a gamma quantile function (inverse CDF):

qgamma(c(.025,.975), 260,50.333) 
[1] 4.556734 5.812086

Notes: (1) R uses the shape/rate parameterization for gamma distributions.

(2) The sampling method, suggested by in the Answer by @philbo_baggins, gives one or two decimal places of accuracy with a sample of 10,000:

set.seed(2022)
quantile(rgamma(10^4, 260, 50.333), c(.025,.975))
    2.5%    97.5% 
4.567375 5.810150 
quantile(rgamma(10^4, 260, 50.333), c(.025,.975))
    2.5%    97.5% 
4.557472 5.809427 

(3) By contrast, one style of frequentist 95% confidence interval gives $(4.531,\ 5.788).$

(256+2 + qnorm(c(.025,.975))*sqrt(256+1))/50
[1] 4.531588 5.788412