Credible interval for gamma prior 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?
 A: 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...
A: 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

