FInding the high density region for a $\chi^2$ (Chi-Squared) distribution?

Question

So I am trying to figure out this problem.

My approach so far has been to consider $\frac{S_o}{\sigma^2} \sim \chi^{2}_k$ as the prior, thus making the posterior $\frac{( S + S_o)}{\sigma^2}\sim\chi^{2}_{n+k}$. Since not much else is given I choose the improper prior: $S_o = 0 $ and $k = 0$ giving a posterior distribution as $\frac{S}{\sigma^2} \sim \chi^{2}_n$ . At this stage I am not sure how to proceed. Any solutions with the steps would be highly appreciated. Thanks !

Update: I calculated the $S= \Sigma_i (x_i-28)^2 = 2588 $ and $\chi^2_{48} \approx 37.689$

Answer 1

Outline of the steps involved.

1. find the posterior density for $(\sigma^2|x)$.

2. see that it's unimodal; note therefore that if you pick some specific value of density $h$, the region with density $\geq h$ will be a contiguous interval. Any such interval will include a proportion of the distribution that you can easily calculate:

   $\qquad\qquad$(A pair of highest density regions for some density $f(x)$)

   A high value of $h$ will lead to a region that includes a small amount of density concentrated close to the mode. A smaller value of $h$ will lead to a region that includes a larger amount of density. With a continuous unimodal ddensity, for any given amount of shaded area (probability) you can find some $h$ that defines the bounds of an interval which includes that much.

3. You need to identify an $\sigma^2_l$ and $\sigma^2_u$ such that $f(\sigma^2_l)=f(\sigma^2_u)$ (i.e. $h$ above) and where the area between them is the desired probability. That interval between those bounds is a HPD interval for $\sigma^2$. There's several (easy) calculations involved in getting the lower and upper bound of the interval.

Answer 2

The answer by Glen_b frames the HPD as a set of simultaneous equations that can be solved via numerical methods. This is one possible way to compute the HPD. An alternative method is to frame the HPD as an optimisation problem, and solve this via numerical methods. The best way to do this depends on the shape of the density function, but in the case of the chi-squared distribution or inverse-gamma distribution (the distribution of its inverse random variable) you have a unimodal continuous density.

Suppose you have a random variable $X \sim f$ where the density is a continuous unimodal (i.e., quasi-concave) function. Denote the density by distribution function by $F$ and the quantile function by $Q$. In this case the HPD will be a contiguous interval. There are two broad ways to derive the interval, both of which have advantages and disadvantages. One optimisation method is to use the density cut-off as your argument variable, which I will call "bottom-to-top" optimisation; the other method is to use one of the tail areas as your argument variable, which I will call "left-to-right" optimisation. I will show how to do the latter optimisation below.

---

Left-to-right optimisation: For any value $0 < \alpha < 1$ the upper bound $U$ can be written as a function of the lower bound $L$ as the function:

$$U(L) = Q(1-\alpha+F(L)) \quad \quad \quad \text{for all } 0 \leqslant L \leqslant Q(\alpha).$$

The highest density region (HDR) can be obtained by solving the optimisation problem:

$$\underset{0 \leqslant L \leqslant Q(\alpha)}{\text{Maximise}} \ \ U(L)-L.$$

This is a non-linear optimisation problem, but it can be solved numerically via standard iterative methods (e.g., by Newton-Raphson iteration). Solving the optimisation problem yields the lower bound for the HDR interval, and the corresponding upper bound can then be obtained by substitution into the above function. These values will satisfy the equations:

$$f(\hat{L}) = f(\hat{U}) \quad \quad \quad \quad \quad \mathbb{P}(\hat{L} \leqslant X \leqslant \hat{U}) = 1-\alpha.$$

If you would like to approach this optimisation problem analytically, you will need to derive the appropriate equation for the iterations for the optimisation (e.g., Newton-Raphson iterations). Alternatively, it is quite simple to program this optimisation problem using the nlm function in R to automate the HDR interval. (I have actually programmed this function myself, but I won't add the code here, since this is a self-study question. I might add it later once some time has elapsed.)

---

Implementation in R: This optimisation method is implemented in the stat.extend package in R, which includes functions for HDRs for a range of univariate probability distributions. In this case your pivotal quantity leads to a confidence interval based on the inverse gamma distribution, which is implemented in the HDR.invgamma function.