Confidence regions from a Gamma distribution

Question

Let a random variable with posterior distribution be given by $X \sim \textrm{Gamma}(10,12)$. This is the result of Jeffreys prior $\pi_J(\lambda)$ multiplied by the Likelihood of i.i.d. $X_1,...,X_n \sim \textrm{Exp}(\lambda)$.

Also, $n=10$ and $\sum_{i=1}^{10} X_i = 12$.

I'm trying to compute a one-sided Bayesian confidence region with level $0.05$ that takes the form $[\alpha,\infty)$, i.e. finding $\alpha$ such that $P(X≥a)=0.95$. This is what I have so far,

$$P(X≥a)=Q(10,\alpha/12) \Rightarrow 12Q^{-1}(10,19/12) \approx 65.1049$$

Where $Q$ is the regularized incomplete gamma function.

The value seems very high to me! Is this actually correct?

This is actually no the part I'm struggling with. I can't manage to compute a one-sided Bayesian confidence region with level $0.05$ that takes the form $(0,b)$, i.e. finding $\alpha$ such that $P(0<X<b)=0.95$. I don't know how to express this as a regularized incomplete gamma function nor inverse cumulative distribution function. If I run some python code, the program eventually fails since it's dealing with both $b$ and $\exp(b)$ which are conflicting and cannot be solved for simultaneously.

Can someone please help me compute $b$?

Answer 1

It appears that you're using the shape-rate parameterization of the gamma. Presumably also $\lambda$ is the rate parameter of the exponential. These will be taken as given for what follows.

Let $\underline{x} = (x_1,x_2,...,x_n)$ be the observed values of $X_1,...,X_n\stackrel{\text{iid}}{\sim}\text{Exp}(\lambda)$. Abusing notation a little, and noting that the Jeffrey's prior for the exponential parameter is $\propto 1/\lambda$ on the positive half line (an improper prior), we have

\begin{eqnarray*} f_{\lambda|\underline{X}=\underline{x}} &\propto& f_{\underline{X}|\lambda}(\underline{x}) \cdot f_\lambda\\ & \propto & \left[\prod_i \lambda e^{-\lambda x_i}\right]\cdot \lambda^{-1} \: I_{(0,\infty)}\\ & = & \lambda^{n-1} e^{-\lambda\sum_i x_i} \: I_{(0,\infty)} \end{eqnarray*}

which we can recognize as the kernel of a $\text{Gamma}(n,\sum_i x_i)$ in the shape-rate parameterization. This much you already knew (though there were multiple issues with the way you were writing it that were causing me some confusion).

Now, a 95% credible interval for $\lambda$ of the form $[a,\infty)$ (again, abusing notation) puts $a$ at the 5th percentile of the distribution. If you have the inverse gamma cdf function, this is simple; with $n=10$ and $\sum_i x_i=12$ you get the value $a\approx 0.4521$ (in R this is qgamma(.05,10,12)). This value makes sense.

(The value $65.1$ you had came from using the wrong parameterization)

With the usual (lower) regularized incomplete gamma, you have to find $a$ such that

$P(as,n)=0.05$

Assuming by $Q$ you are looking at the upper rather than the lower regularized incomplete gamma, as with the notation here:

https://en.wikipedia.org/wiki/Incomplete_gamma_function#Regularized_gamma_functions_and_Poisson_random_variables

and equivalently the upper regularized gamma here

https://mathworld.wolfram.com/RegularizedGammaFunction.html

you'd instead be solving $Q(as,n)=0.95$. That is, finding $a$ such that $Q(12a,10) = 0.95$. This should give the $0.452$-ish value above.

For the other tail, you'd solve $Q(12b,10) = 0.05$

From the values you give, it's not quite clear to me how your inverse function $Q^{-1}$ is set up.