# Confidence regions from a Gamma distribution

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. 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$$?

• By "Bayesian confidence region" do you instead intend a credible interval? ... en.wikipedia.org/wiki/Credible_interval Commented Dec 18, 2022 at 7:11
• Then (in relation to the use of 'confidence') it seems the error is not yours but someone else's. That now raises another question. Is this coursework (or other self-study)? If so, please see the self-study tag wiki info ... stats.stackexchange.com/tags/self-study/info Commented Dec 18, 2022 at 7:17
• If I am ready to give an answer, rest assured you will have it. It's not like I am waiting to be prompted in order to begin answering the question. I want to be sure I am clear what is being asked. It would help if you would clarify which parameterization you're using. Can you clarify what $\alpha$ is in $P(X≥a)=Q(10,\alpha/12)$ and where $a$ went? Commented Dec 18, 2022 at 7:24
• Please add any info provided in the comments in the question body. Also as Glen_b said, please add the self-study tag. Commented Dec 18, 2022 at 8:10
• Done @User1865345 :) Commented Dec 18, 2022 at 8:11

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.

• Thanks for your answer @Glen_b. I'm getting 0.468307637002674 when solving for $b$ in $𝑄(12𝑏,10)=0.05$. This is wrong. Any idea why? Commented Dec 18, 2022 at 13:39
• It looks like your Q function's arguments are not in the same order as the definition I used from the links that I gave. Swap them to match the reversed order in your function. You may also need to pay attention to the exact way the Q-function is defined to make sure it is definitely the upper regularized incomplete gamma. Commented Dec 18, 2022 at 13:47
• The answer should be roughly 1.3, I think Commented Dec 18, 2022 at 13:50
• Yes, that's the correct answer I've got mine accurate to 3dp. The only problem is that the regularized incomplete gamma function returns a univariate polynomial of degree $9$ in terms of $b$ multiplied by $e^{-12b}$ which only wolfram alpha can solve. My python program still breaks as it isn't able to solve such complex equations which, to be honest, really sucks. I'm considering switching to R instead. Commented Dec 18, 2022 at 14:02
• wouldn't scipy.stats.gamma.ppf do what you need in python? Commented Dec 18, 2022 at 14:28