HPD interval for the mean

Question

Suppose we have iid observation with the following model $ Y_t \sim \mathcal{N}(\mu,1/\mu) , t=1,2,..T$

The question is " Assuming a flat prior on $(0 ,\infty )$ find a 95% HPD interval for $\mu"$

The aditional information that the problem gives is: $T=6; \sum{y_i}^{2}=76;\sum{y_i}=18$

---

this is my first try, $p(\mu)=\mu^{-1}1_{y>0}$ then

$ \begin{eqnarray*} p(\mu\vert Y) & \propto & p(Y \vert \mu) p(\mu) & \propto \mu^{T}\exp[\sum{(-\mu/2){}(y_i-\mu)}^{2} ] \mu^{-1}1_{y>0} \end{eqnarray*} $

$Pr(\mu ϵ C│y)=\int_{a}^{b}p(\mu\vert Y)d\mu=\int_{a}^{b} \mu^{T-1}\exp[\sum{(-\mu/2){}(y_i-\mu)}^{2} ] 1_{y>0}d\mu=0.95$

---

I think this problem doesn´t have an analytical solution, so the solution is numeric. The question is which algorithm ,step by step, can I use to solve it .

---

Update

$\sum{(y_i-\mu)}^{2}= \sum{y_i}^{2}-2\mu\sum{y_i}+T\mu^{2}=76-36\mu+6\mu^{2}$ $\ (-\mu/2) \sum{(y_i-\mu)}^{2}= -38\mu+18\mu^{2}-3\mu^{3}$

$Pr(\mu ϵ C│y)=\int_{a}^{b}p(\mu\vert Y)d\mu=\int_{a}^{b} \mu^{5}\exp[-38\mu+18\mu^{2}-3\mu^{3} ] 1_{y>0}d\mu=0.95$

    %%Matlab Code
    a=0.0001;b=0.0001;step=0.0001;
    x=[a:step:1];
    fx= x.^(5).*exp(-38.*x+18.*(x.^2)-3.*(x.^3)); plot(x,fx)

How can I find (a,b) that give me a 0.95 credible region with the minimum distiance between a and b

Answer 1

In brief, given an array of $x, p(x)$ pairs:

• Sort the array by $p(x)$, decreasing.
• Find the first element for which the cumulative sum of density of the sorted array is $\ge 1 - \alpha$.
• Find the array elements with $p(x)\ge$ this element's density. These $x$ values form the highest density interval(s).

Here's an example implementation in R, from the code associated with Doing Bayesian Data Analysis:

HDIofGrid = function( probMassVec , credMass=0.95 ) {
   # Arguments:
   #   probMassVec is a vector of probability masses at each grid point.
   #   credMass is the desired mass of the HDI region.
   # Return value:
   #   A list with components:
   #   indices is a vector of indices that are in the HDI
   #   mass is the total mass of the included indices
   #   height is the smallest component probability mass in the HDI
   # Example of use: For determining HDI of a beta(30,12) distribution
   #   approximated on a grid:
   #   > probDensityVec = dbeta( seq(0,1,length=201) , 30 , 12 )
   #   > probMassVec = probDensityVec / sum( probDensityVec )
   #   > HDIinfo = HDIofGrid( probMassVec )
   #   > show( HDIinfo )
   sortedProbMass = sort( probMassVec , decreasing=T )
   HDIheightIdx = min( which( cumsum( sortedProbMass ) >= credMass ) )
   HDIheight = sortedProbMass[ HDIheightIdx ]
   HDImass = sum( probMassVec[ probMassVec >= HDIheight ] )
   return( list( indices = which( probMassVec >= HDIheight ) ,
                 mass = HDImass , height = HDIheight ) )
}

My own in python, adapted from the above:

def hdi(x, fx, alpha):
    """
    Given x and fx, returns hdi
    """
    sfx = fx.copy()
    sfx.sort()
    sfx[:] = sfx[::-1]
    lowestHeightIdx = np.min(np.where(np.cumsum(sfx) > (1. - alpha)))
    lowestHeight = sfx[lowestHeightIdx]
    return x[fx >= lowestHeight]