HPD interval for the mean

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

• You need a numerical approximation of the inverse cumulative function, otherwise you can simulate from the posterior and use a procedure to get the HPD based on a sample. Aug 24 '14 at 21:38
• @StéphaneLaurent Thanks for your answer. Can you give some code, or some reference to an example with an analogous solution? Aug 24 '14 at 22:00
• I don't know Matlab. Aug 25 '14 at 5:53
• Presumably "HPD" = "Highest Posterior Density". The interval you find should extend from $0.04713$ to $0.3946$.
– whuber
Aug 25 '14 at 14:42
• Thanks @whuber for your answer. I don´t know how you find such interval ; May be I was not clear; my question is :how to find a pseudocode (or a code in any language) that give me the exact length of the interval Aug 25 '14 at 17:09

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]