So I am trying to figure out this problem.

enter image description here

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$

  • 1
    $\begingroup$ Please read the self-study tag wiki which explains how self study questions work. In particular, it explains you should not be provided with step-by-step solutions. Indeed first you need to tell us what you do understand and where you first run into difficulty. [If you really have no idea, you aren't in a position to ask yet -- you should go back and reread your notes, and draw some pictures.] $\endgroup$
    – Glen_b
    Nov 1 '15 at 23:53
  • $\begingroup$ @Glen_b Updated with some partial work. Any help would be appreciated ! $\endgroup$
    – Jenna Maiz
    Nov 3 '15 at 4:47
  • $\begingroup$ You have a density for $S/\sigma^2$ ... can you identify a posterior density for $\sigma^2$? $\endgroup$
    – Glen_b
    Nov 3 '15 at 4:57
  • $\begingroup$ Not sure. Do we really need that ? $\endgroup$
    – Jenna Maiz
    Nov 3 '15 at 4:59
  • $\begingroup$ How would you propose identifying the region of the posterior for $\sigma^2$ that has the highest density without knowing what that density is? [Surely you must have done some Bayesian calculation before this, and have seen examples of others. What happened there? --- and if you have not, you are not in a position to ask for help with this question yet.] $\endgroup$
    – Glen_b
    Nov 3 '15 at 5:01

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:

    enter image description here $\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.


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.

  • $\begingroup$ There is even an R package (on CRAN) doing this for you. For an example, see stats.stackexchange.com/questions/171458/… $\endgroup$ Sep 21 '19 at 9:45
  • 1
    $\begingroup$ @kjetil: That could have saved me some programming time! $\endgroup$
    – Ben
    Sep 21 '19 at 10:59
  • $\begingroup$ Glen_b's method actually is an advance over this approach because it implicitly recognizes that the derivative of the objective can easily be expressed in terms of the densities at the endpoints and therefore all critical points occur for intervals where those densities are equal. This is a substantial simplification and, as he suggests, leads to effective and efficient numerical algorithms. $\endgroup$
    – whuber
    Aug 14 '20 at 14:05
  • 1
    $\begingroup$ Ben, I don't recognize @Glen's answer in your straw-man description of it. He's clear it's an "outline of the steps involved." I didn't use the word "method" in my comment as if it meant a capital-M Method. Your approach, although perfectly correct, requires more analysis and work, beginning with having to formulate the function $U.$ The diagrams in Glen's answer make the problem clear and also point to a much simpler way to compute the NR iterations. His approach also generalizes in an obvious manner to multimodal PDFs, whereas your formulation cannot: $U$ does not exist then. $\endgroup$
    – whuber
    Aug 15 '20 at 14:15
  • 1
    $\begingroup$ Just to clarify; my answer was deliberately given as "an outline" because the question was clearly a homework style question. Explicit algebraic steps are easy to identify given the outline of the reasoning but I left them for the student. Indeed the basic idea is easy for students to remember, even when specific calculations are forgotten. $\endgroup$
    – Glen_b
    Aug 17 '20 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.