Probability Interval for Gamma Distribution

Question

There is something that I am doing wrong in the exercise below and I would appreciate some help figuring it out.

Let $X_1, X_2, \ldots, X_5$ be a random sample from a $\Gamma \left(3,3 \right)$ distribution . Define $W=\sum_{i=1}^5 X_i $.

1. Obtain the distribution of $W$ and of $2W/3$

2. Find $c_1$ and $c_2$ such that $P \left( c_1 < W <c_2 \right)=0.95 $

---

Since $W$ is the sum of iid gamma variables, $W \sim \Gamma \left(15,3 \right)$ and now a straightforward transformation argument shows that $Y=\frac{2}{3} W \sim \chi^2 \left(30 \right) $. This completes question 1.

Now for question 2, my thoughts were to transform these probabilities using the result from question 1, namely that $\frac{2}{3} W \sim \chi^2 \left(30 \right) $. Hence, unless I have made a horrible mistake here:

$$P \left( \frac{2}{3} c_1 < Y < \frac{2}{3} c_2 \right)=0.95 $$

We can find the constants from a chisquare table for $30$ degrees of freedom. In my case, they are $16.791$ for the lower bound and $46.979$ for the upper one.

Then if we equate $\frac{2}{3}c_1=16.791$ and solve for $c_1$ we should be able to find the constant for the above gamma confidence interval, right?

The problem is that according to $R$, $P \left( W<c_1 \right) =1 $, instead of the desired $0.025$ which means I have gone wrong somewhere. I do not think it is in question though $1$ as I have verified the result many times.

Thank you.

Answer 1

Your algebra sounds like it's all fine.

You don't give enough information to guess how you decided "according to $\text{R}$, $P(W<c_1)=1$, instead of the desired $0.025$".

This is a common issue with the gamma distribution - there are two different common parameterizations, both reasonably widespread, and if we're not careful, we can think we're dealing with one when we're actually doing the other. (Actually, there's a third paramaterization that comes up a lot when dealing with gamma GLMs, the shape-mean parameterization, but that one is usually more obvious when it occurs.)

Wikipedia (permalink to today's version) gives both forms, see the column on the right. Confusingly, it swaps the role of what I see as the more conventional parameter names (to my mind $\beta$ is more often the scale, $\theta$ is more often the rate).

[While it comes up a lot when dealing with software, it's not simply a software issue, because the issue often happens between humans as well.]

As it happens, R is a particular culprit with this kind of issue, the help for the collection of gamma distribution functions seemingly going out of its way to muddy the water (I'm using 3.0.2 at the time of writing, but the issue has been there for ages).

My guess is you might be calling the pgamma function in R with unnamed arguments, but supplying shape and scale arguments, like so ...

 pgamma(16.791*3/2,15,3)
[1] 1

... when R defaults to shape and rate arguments.

When in doubt, draw. Here's what you want to be calculating with:

As you see, the proportion of that to the left of $\frac{3}{2}16.791$ (around 25) looks roughly right.

Here's the actual density I suspect you're calculating with:

plot(x,dgamma(x,15,3),col=3,type="l")

Here the proportion of that to the left of $\frac{3}{2}16.791$ will be very close to 1.

The R help is, unfortunately, less than clear -- even actively misleading. The sentence under "Description" implies that the supplied parameters are shape and scale -- and the description of rate under "Arguments" confirms the impression(!) -- but the argument list in the function itself isn't ambiguous, if non-obvious to a novice user:

pgamma(q, shape, rate = 1, scale = 1/rate, lower.tail = TRUE, log.p = FALSE)

Do you see how the scale defaults to a function of whatever is specified for the rate? If you just supply two unnamed arguments after the q, they are the shape and the rate, and then scale is obtained by taking reciprocals.

Which is to say, you need to use a named argument to get what you want:

 pgamma(16.791*3/2,15,scale=3)
[1] 0.02500255

When there is any potential for doubt, you should probably name your arguments anyway, to make them explicit to human readers.