Intervalo HPD para la media

Question

Supongamos que tenemos iid observación con el siguiente modelo de $ Y_t \sim \mathcal{N}(\mu,1/\mu) , t=1,2,..T$

La pregunta es " Suponiendo un plano antes de $(0 ,\infty )$ encontrar un 95% de HPD intervalo de $\mu"$

La información adicional que el problema que da es: $T=6; \sum{y_i}^{2}=76;\sum{y_i}=18$

---

este es mi primer intento, $p(\mu)=\mu^{-1}1_{y>0}$

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

---

Creo que este problema no tiene una solución analítica, por lo que la solución es numérico. La cuestión es que el algoritmo ,paso a paso, se puede utilizar para resolverlo .

---

Actualización

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

¿Cómo puedo encontrar (a,b) que me dan un 0.95 creíble región con el mínimo distiance entre a y b

Answer 1

En breve, dada una matriz de $x, p(x)$ pares de:

• Ordenar la matriz por $p(x)$, disminuyendo.
• Encontrar el primer elemento para el que la suma acumulativa de la densidad de la matriz ordenada es $\ge 1 - \alpha$.
• Encontrar los elementos de la matriz con $p(x)\ge$ este elemento de la densidad. Estos $x$ formulario valores de la más alta densidad en el intervalo(s).

He aquí un ejemplo de implementación en R, desde el código asociado con Hacer Bayesiano de Análisis de Datos:

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 ) )
}

Mi propia en python, una adaptación de la anterior:

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]