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 .


$\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
    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

  • 1
    $\begingroup$ 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. $\endgroup$ Aug 24 '14 at 21:38
  • $\begingroup$ @StéphaneLaurent Thanks for your answer. Can you give some code, or some reference to an example with an analogous solution? $\endgroup$
    – EAguirre
    Aug 24 '14 at 22:00
  • $\begingroup$ I don't know Matlab. $\endgroup$ Aug 25 '14 at 5:53
  • 2
    $\begingroup$ Presumably "HPD" = "Highest Posterior Density". The interval you find should extend from $0.04713$ to $0.3946$. $\endgroup$
    – whuber
    Aug 25 '14 at 14:42
  • $\begingroup$ 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 $\endgroup$
    – EAguirre
    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[:] = sfx[::-1]
    lowestHeightIdx = np.min(np.where(np.cumsum(sfx) > (1. - alpha)))
    lowestHeight = sfx[lowestHeightIdx]
    return x[fx >= lowestHeight]

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.