1. La siguiente tabla muestra el número de accidentes mortales por año asociados con aerolíneas en todo el mundo durante un período de diez años (Fuente: Statistical Abstract of the United States). Suponga que el número de accidentes mortales en cada año son condicionalmente independientes y siguen una distribución Poisson con parámetro \(\theta\). Establezca una distribución previa para \(\theta\) y determine la distribución posterior con base en los datos de 1976 a 1985. Bajo este modelo, calcule un intervalo predictivo al 95% para el número de accidentes fatales en 1986.
Año 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985
Accidentes 24 25 31 31 22 21 26 20 16 22

  1. Un laboratorio está estimando la tasa de tumorigenesis en dos cepas de ratones, A y B. Los ratones tipo A han sido bien estudiados e información de otros laboratorios sugiere que los ratones tipo A tienen conteos de tumores que siguen una distribución de Poisson con media \(\theta_A = 12\). Se desconoce la tasa promedio de los tumores para los ratones tipo B, \(\theta_B\), pero existe suficiente evidencia empírica para asegurar que los ratones tipo B están relacionados con los ratones tipo A. Los conteos de tumores observados para las dos cepas de ratones son \(\boldsymbol{y}_A = (12, 9, 12, 14, 13, 13, 15, 8, 15, 6)\) y \(\boldsymbol{y}_B = (11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)\).

    1. Encuentre las distribuciones posteriores, las medias posteriores, las desviaciones estándar posteriores y los intervalos de credibilidad del 95% para \(\theta_A\) y \(\theta_B\), asumiendo modelos Gamma-Poisson independientes para cada grupo, con distribuciones previas \(\theta_A\sim\textsf{Gamma}(120,10)\) y \(\theta_B\sim\textsf{Gamma}(12,1)\). ¿Por qué estas distribuciones previas son razonables?

    2. Calcule y grafique la media posterior de \(\theta_B\) bajo la distribución previa \(\theta_B\sim \textsf{Gamma}(12m,m)\), para cada valor de \(m\in\{1, 2,\ldots,50\}\). Describa qué tipo de información previa sobre \(\theta_B\) sería necesaria para que la media posterior de \(\theta_B\) sea cercana a la de \(\theta_A\).

  2. Muestre que bajo el modelo Gamma-Poisson \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\), para \(i = 1,\ldots,n\), con \(\theta \sim \textsf{Gamma}(a,b)\), se tiene que la distribución predictiva posterior es Binomial Negativa con parámetros \(a+s\) y \(b+n\), donde \(s=\sum_{i=1}^n y_i\).

  3. Considere el modelo Gamma-Poisson \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\), para \(i = 1,\ldots,n\), con \(\theta \sim \textsf{Gamma}(\alpha,\beta)\), y la función de perdida cuadrática \(L(\theta,\theta^*)=(\theta-\theta^*)^2\).

    1. Muestre que el estimador que minimiza la perdida esperada posterior es de la forma \(\hat\theta = a + b\,\bar{y}\), donde \(a > 0\), \(b \in (0, 1)\) y \(\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\).
    2. Calcule el riesgo frecuentista tanto de \(\hat\theta\) como del estimador de máxima verosimilitud.
    3. Calcule el riesgo Bayesiano de \(\hat\theta\).
    4. Suponga que un investigador quiere recolectar una muestra lo suficientemente grande para que el riesgo Bayesiano después del experimento sea la mitad del riesgo Bayesiano antes del experimento. Encuentre ese tamaño de muestra.

  4. Considere el modelo Beta-Binomial \(x\mid\theta\sim \textsf{Bin}(n,\theta)\) con \(\theta\sim \textsf{Beta}(\sqrt{n}/2,\sqrt{n}/2)\), donde \(n\) es conocido.

    1. ¿Cuál es la distribución posterior?
    2. ¿Cuál es el estimador que minimiza la perdida esperada posterior si la función de perdida es \(L(\theta,\hat\theta)=(\theta-\hat\theta)^2\)? Llame tal estimador \(\hat{\theta}\) y muestre que el riesgo frecuentista correspondiente es constante.
    3. Sea \(\theta_0=x/n\). Encuentre el riesgo frecuentista de este estimador. Compare los riesgos de \(\hat{\theta}\) y \(\theta_0\), para \(n=10,50,100\). ¿Qué se puede concluir?

  5. Sea \(x\sim \textsf{Bin}(n,\theta)\) con \(n\) conocido, y \(\theta \sim \textsf{Beta}(\alpha,\beta)\). Considere la función de perdida \(L(\theta,\hat\theta)=(\theta-\hat\theta)/(\theta(1-\theta))\). Encuentre el estimador que minimiza la perdida esperada posterior para esta función de perdida.

  6. Sea \(L(\theta,\hat\theta) = \omega(\theta)(\theta - \hat\theta)^2\) la función de perdida cuadrática ponderada, donde \(\omega(\theta)\) es una función no negativa. Muestre que el estimador que minimiza la perdida esperada posterior tiene la forma \[ \hat\theta=\frac{\textsf{E}(\omega(\theta)\,\theta\mid \boldsymbol{y})}{\textsf{E}(\omega(\theta)\mid \boldsymbol{y})}\,. \]

  7. Un investigador del Departamento de Ingeniería Electrónica y Eléctrica de la Universidad de Bath (Inglaterra) necesitaba analizar unos datos sobre los tiempos de falla de un determinado tipo de alambre de metal (en este problema el tiempo de falla se define como el número de veces que una máquina podría tensionar el alambre antes de romperse). Los siguientes datos corresponden a \(n = 14\) tiempos de falla de una parte del experimento: \[ \boldsymbol{y} = (495,541,1461,1555,1603,2201,2750,3468,3516,4319,6622,7728,13159,21194)\,. \] El modelo más simple para los datos del tiempo de falla involucra la distribución Exponencial, \(y_i \mid \lambda \stackrel{\text{iid}}{\sim} \textsf{Exponencial} ( \lambda )\), con \(\textsf{E}(y_i \mid \lambda) = \lambda\), para \(i=1,\ldots,n\).

    1. Sea \(X \sim \textsf{Gamma} (\alpha, \beta)\), con \(\alpha > 0\) y \(\beta > 0\). Encuentre la función de densidad de probabilidad de \(Y=\frac{ 1 }{ X }\). La distribución de \(Y\) se denomina distribución Gamma Inversa con parámetros \(\alpha\) y \(\beta\), lo que se denota con \(Y\sim\textsf{GI} (\alpha,\beta)\).

    2. Encuentre la distribución posterior de \(\lambda\), asumiendo que \(\lambda\sim\textsf{GI}(a,b)\), con \(a > 0\) y \(b > 0\).

    3. Muestre que la media posterior de \(\lambda\) es un promedio ponderado de la media previa \(\textsf{E}(\lambda)\) y la media muestral \(\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\).

    4. Se tiene información externa de otro experimento de acuerdo con el cual la distribución previa de \(\lambda\) debería tener una media \(\mu_0 = 4500\) y una desviación estándar \(\sigma_0 = 1800\). Grafique las distribuciones previa y posterior de \(\lambda\) en el mismo gráfico con \(\lambda\) en el eje horizontal tomando valores de 1000 a 12000, identificando qué curva corresponde a qué densidad.

    5. Muestre que el estimador de máxima verosimilitud (MLE, por sus siglas en inglés) de \(\lambda\) y la información observada de Fisher son respectivamente \[ \hat\lambda_{\text{MLE}} = \bar{y} \qquad\text{y}\qquad \hat{I}(\hat\lambda_{\text{MLE}}) = \frac{n}{\bar{y}^2}\,. \]

      Recuerde que la información observada de Fisher se define como \[ \hat{I}(\hat\lambda_{\text{MLE}}) = \left[ -\frac{\partial^2}{\partial\lambda^2}\log p(\boldsymbol{y}\mid\lambda) \right]_{\lambda = \hat\lambda_{\text{MLE}}}\,, \] donde \(p(\boldsymbol{y}\mid\lambda) = \prod_{i=1}^n p(y_i\mid\lambda)\) es la distribución muestral conjunta de \(\boldsymbol{y}\).

    6. Encuentre la estimación, el coeficiente de variación y un intervalo de de credibilidad/confianza al 95% para \(\lambda\), bajo el paradigma Bayesiano como Frecuentista (usando la Normalidad asintóica del MLE, Bootstrap paramétrico y Bootstrap no paramétrico).

      Recuerde que asintóticamente se tiene que \(\hat\lambda_{\text{MLE}}\sim\textsf{N}(\lambda,\hat{I}^{-1}(\hat\lambda_{\text{MLE}}))\).

Solución

  1. La siguiente tabla muestra el número de accidentes mortales por año asociados con aerolíneas en todo el mundo durante un período de diez años (Fuente: Statistical Abstract of the United States). Suponga que el número de accidentes mortales en cada año son condicionalmente independientes y siguen una distribución Poisson con parámetro \(\theta\). Establezca una distribución previa para \(\theta\) y determine la distribución posterior con base en los datos de 1976 a 1985. Bajo este modelo, calcule un intervalo predictivo al 95% para el número de accidentes fatales en 1986.
Año 1976 1977 1978 1979 1980 1981 1982 1983 1984 1985
Accidentes 24 25 31 31 22 21 26 20 16 22

Sea \(y_i\) el número de accidentes mortales en el año \(i\) para \(i=1,\ldots,n\), donde el año \(i=1\) corresponde a 1976, el año \(i=2\) a 1977, etc. Asumiendo que el número de accidentes mortales en cada año son condicionalmente independientes y siguen una distribución Poisson con parámetro \(\theta\), se tiene que \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\). Así, se considera una distribución previa conjuga para \(\theta\) de la forma \(\theta\sim\textsf{Gamma}(a,b)\) con \(a=1\) y \(b=1\). Esta previa es una distribución previa razonable ya que es aproximadamente constante para \(\theta > 10\) (en este contexto es aceptable asumir que la tasa promedio del número de accidentes mortales será considerablemente superior a 10), aunque se encuentre concentrada débilmente al rededor de 1 pues \(\textsf{E}(\theta)=1\) y \(\textsf{CV}(\theta) = 1\) (otra manera de haber construido la previa consistiría en elegir \(a\) y \(b\) tales que \(\textsf{E}(\theta)=\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\) y \(\textsf{CV}(\theta) = 1\); aunque estrictamente esto no es válido porque se están usando los datos para estipular el estado de conocimiento anterior al proceso de observación, hay una rama de la estadística Bayesiana denominada Bayes Empírico [empirical Bayes] que construye las distribuciones previas de esta manera).

# hiperparametros
a <- 1
b <- 1
# previa
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
curve(expr = dgamma(x, shape = a, rate = b), from = 0, to = 40, n = 1000, lwd = 2, 
      xlab = expression(theta), ylab = "Densidad", main = "Distr. Previa")

Así, la distribución posterior es \(\theta\mid\boldsymbol{y}\sim\textsf{Gamma}(a+s,b+n)\), donde \(\boldsymbol{y}=(y_1,\ldots, y_n)\) es el conjunto de datos observado y \(s=\sum_{i=1}^n y_i\) es el estadístico suficiente correspondiente. En este caso, \(n=10\) y \(s=238\) y en consecuencia, \(\theta\mid\boldsymbol{y}\sim\textsf{Gamma}(239,11)\). Luego, la estimación puntual de \(\theta\) es \(\hat\theta=\textsf{E}(\theta\mid\boldsymbol{y}) = 21.727\) y el intervalo de credibilidad correspondiente al 95% es \((19.060\,;\,24.567)\).

# datos
y <- c(24, 25, 31, 31, 22, 21, 26, 20, 16, 22)
n <- length(y)
print(n)
## [1] 10
# estadistico suficiente
s <- sum(y)
print(s)
## [1] 238
# posterior
ap <- a + s
print(ap)
## [1] 239
bp <- b + n
print(bp)
## [1] 11
# grafico
par(mfrow = c(1,2), mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
curve(expr = dgamma(x, shape = ap, rate = bp), from = 0, to = 40, n = 1000, lwd = 2, col = 4, 
      xlab = expression(theta), ylab = "Densidad", main = "Distr. posterior")
curve(expr = dgamma(x, shape = a, rate = b), n = 1000, add = T)
legend("topright", legend = c("Posterior","Previa"), col = c(1,4), lwd = 2, bty = "n")
# 
curve(expr = dgamma(x, shape = ap, rate = bp), from = 16, to = 28, n = 1000, lwd = 2, col = 4, 
      xlab = expression(theta), ylab = "Densidad", main = "Distr. posterior")
curve(expr = dgamma(x, shape = a, rate = b), n = 1000, add = T)
abline(v = ap/bp, lwd = 2, lty = 2, col = 2)
abline(v = qgamma(p = c(0.025, 0.975), shape = ap, rate = bp), lwd = 2, lty = 3, col = 3)
legend("topright", legend = c("Media","IC 95%"), col = c(2,3), lwd = 2, bty = "n")

# estimacion puntual
round(ap/bp, 3)
## [1] 21.727
# intervalo de creidibilidad al 95%
round(qgamma(p = c(0.025, 0.975), shape = ap, rate = bp), 3)
## [1] 19.060 24.567

Finalmente, un análisis de sensibilidad muestra que los resultados no cambian radicalmente al usar la previa difusa \(a=1\) y \(b=1\) (i.e., \(\textsf{E}(\theta)=1\) y \(\textsf{CV}(\theta) = 1\) con densidad aproximadamente uniforme para \(\theta > 10\)) en comparación la previa empírica \(a=1\) y \(b=1/\bar{y}=0.042\) (i.e., \(\textsf{E}(\theta)=\bar{y}=23.8\) y \(\textsf{CV}(\theta) = 1\)). La siguiente tabla muestra la estimación puntual de \(\theta\) junto con el intervalo de credibilidad correspondiente al 95%, utilizando estas dos distribuciones previas (la previa de Jeffreys sería otra alternativa para llevar acabo el análisis de sensibilidad).

# previa difusa (previa 1)
a1 <- 1
b1 <- 1
# previa empírica (previa 2)
a2 <- 1
b2 <- 1/mean(y)
# posterior usando previa difusa
ap1 <- a1 + s
bp1 <- b1 + n
# posterior usando previa empírica
ap2 <- a2 + s
bp2 <- b2 + n
# inferencia
tab <- rbind(c(ap1/bp1, qgamma(p = c(0.025, 0.975), shape = ap1, rate = bp1)), 
             c(ap2/bp2, qgamma(p = c(0.025, 0.975), shape = ap2, rate = bp2)))
rownames(tab) <- c("Previa difusa","Previa empírica")
colnames(tab) <- c("Estimación","Q2.5%","Q97.5%")
knitr::kable(x = tab, digits = 3)
Estimación Q2.5% Q97.5%
Previa difusa 21.727 19.060 24.567
Previa empírica 23.800 20.878 26.911
  1. Un laboratorio está estimando la tasa de tumorigenesis en dos cepas de ratones, A y B. Los ratones tipo A han sido bien estudiados e información de otros laboratorios sugiere que los ratones tipo A tienen conteos de tumores que siguen una distribución de Poisson con media \(\theta_A = 12\). Se desconoce la tasa promedio de los tumores para los ratones tipo B, \(\theta_B\), pero existe suficiente evidencia empírica para asegurar que los ratones tipo B están relacionados con los ratones tipo A. Los conteos de tumores observados para las dos cepas de ratones son \(\boldsymbol{y}_A = (12, 9, 12, 14, 13, 13, 15, 8, 15, 6)\) y \(\boldsymbol{y}_B = (11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)\).
  1. Encuentre las distribuciones posteriores, las medias posteriores, las varianzas posteriores y los intervalos de credibilidad del 95% para \(\theta_A\) y \(\theta_B\), asumiendo modelos Gamma-Poisson independientes para cada grupo, con distribuciones previas \(\theta_A\sim\textsf{Gamma}(120,10)\) y \(\theta_B\sim\textsf{Gamma}(12,1)\). ¿Por qué estas distribuciones previas son razonables?

Estas distribuciones previas son razonables porque están centradas al rededor de 12 unidades dado que \(\textsf{E}(\theta_A) = 120/10 = 12\) y \(\textsf{E}(\theta_B) = 12/1 = 12\) (la evidencia empírica indica que los ratones tipo A han sido bien estudiados y tienen media 12 y los ratones tipo B están relacionados con los ratones tipo A). La diferencia radica en el grado de concentración de las densidades a priori al rededor de 12, pues \(\textsf{CV}(\theta_A)=1/\sqrt{120}\approx 9\%\) mientras que \(\textsf{CV}(\theta_B)=1/\sqrt{12}\approx 29\%\). Esto último concuerda con el estado de información externa al conjunto de datos de acuerdo con el contexto del problema, porque estamos más seguros acerca de la evidencia previa acerca de los ratones tipo A en comparación con los ratones tipo B.

# grafico
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
curve(expr = dgamma(x, shape = 120, rate = 10), from = 0, to = 25, n = 1000, col = 2, lwd = 1, ylim = c(0,0.5), xlab = expression(theta), ylab = "Densidad", main = "Distr. Previa")
curve(expr = dgamma(x, shape = 12, rate = 1), n = 1000, col = 4, lwd = 1, add = TRUE)
abline(v = 12, lty = 2)
legend("topright", legend = c("Tipo A", "Tipo B"), col = c(2, 4), lty = 1, lwd = 2, bty = "n")

Así, la distribución posterior en ambos casos es tipo Gamma. En particular, se tiene que las distribuciones posteriores correspondientes son \(\theta_A\mid\boldsymbol{y}_A\sim\textsf{Gamma}(237,20)\) y \(\theta_B\mid\boldsymbol{y}_B\sim\textsf{Gamma}(125,14)\).

# datos
yA <- c(12, 9, 12, 14, 13, 13, 15, 8, 15, 6)
yB <- c(11, 11, 10, 9, 9, 8, 7, 10, 6, 8, 8, 9, 7)
# tamaños
nA <- length(yA)
nB <- length(yB)
# estadístico suficiente
sA <- sum(yA)
sB <- sum(yB)
# previa
aA <- 120
bA <- 10
aB <- 12
bB <- 1
# posterior
apA <- aA + sA
bpA <- bA + nA
print(c(apA, bpA))
## [1] 237  20
apB <- aB + sB
bpB <- bB + nB
print(c(apB, bpB))
## [1] 125  14

Por tal motivo las medias posteriores, las varianzas posteriores, y los intervalos de credibilidad correspondientes al 95% son los que se presentan en la tabla a continuación.

# inferencia
tab <- rbind(c(apA/bpA, apA/bpA^2, qgamma(p = c(0.025, 0.975), shape = apA, rate = bpA)), 
             c(apB/bpB, apB/bpB^2, qgamma(p = c(0.025, 0.975), shape = apB, rate = bpB)))
rownames(tab) <- c("Tipo A","Tipo B")
colnames(tab) <- c("Media","Varianza","Q2.5%","Q97.5%")
knitr::kable(x = tab, digits = 3)
Media Varianza Q2.5% Q97.5%
Tipo A 11.850 0.593 10.389 13.405
Tipo B 8.929 0.638 7.432 10.560 
# gráfico
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
# tipo A
curve(expr = dgamma(x, shape = apA, rate = bpA), from = 0, to = 25, n = 1000, col = 2, lwd = 2, ylim = c(0,0.5), xlab = expression(theta), ylab = "Densidad", main = "")
curve(expr = dgamma(x, shape = aA,  rate = bA),  n = 1000, col = 2, lwd = 1, lty = 2, add = TRUE)
curve(expr = dgamma(x, shape = apB, rate = bpB), n = 1000, col = 4, lwd = 2, lty = 1, add = TRUE)
curve(expr = dgamma(x, shape = aB,  rate = bB),  n = 1000, col = 4, lwd = 1, lty = 2, add = TRUE)
legend("topright", legend = c("Posterior A","Posterior B","Previa A","Previa B"), col = c(2,4,2,4), lty = c(1,1,2,2), lwd = c(2,3,1,1), bty = "n")

  1. Calcule y grafique la media posterior de \(\theta_B\) bajo la distribución previa \(\theta_B\sim \textsf{Gamma}(12m,m)\), para cada valor de \(m\in\{1, 2,\ldots,50\}\). Describa qué tipo de información previa sobre \(\theta_B\) sería necesaria para que la media posterior de \(\theta_B\) sea cercana a la de \(\theta_A\).

Como se evidencia en las Figuras, la información previa sobre \(\theta_B\) para que \(\textsf{E}(\theta_B\mid\boldsymbol{y}_B)\) sea cercana a \(\textsf{E}(\theta_A\mid\boldsymbol{y}_A)=11.850\) (línea en color negro) debe ser muy informativa al rededor de \(\textsf{E}(\theta_B)=12m/m=12\), lo cual se logra incrementando el valor de \(m\), de forma que \(\textsf{Var}(\theta_B)=12m/m^2=12/m\) decrezca haciendo que la información previa sea mas fuerte.

# gráfico
mgrid <- c(1,8,15,25,50)
col <- RColorBrewer::brewer.pal(n = 9, name = "Set1")
par(mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
plot(NA, NA, xlim = c(6,14), ylim = c(0,1), xlab = expression(theta), ylab = "Densidad posterior", main = "Tipo B")
for (j in 1:length(mgrid)) {
  curve(expr = dgamma(x, shape = 12*mgrid[j] + sB, rate = 1*mgrid[j] + nB),  n = 1000, col = col[j], lwd = 2, lty = 1, add = TRUE)
  abline(v = (12*mgrid[j] + sB)/(1*mgrid[j] + nB), lty = 2,  col = col[j])
}
abline(v = (12 + sA)/(1 + nA), lty = 1)
legend("topright", legend = c(paste0("m = ", mgrid), "Media post. A"), col = c(col[1:length(mgrid)], "black"), lty = 1, lwd = 2, bty = "n")

mgrid <- 1:50
mmean <- (12*mgrid + sB)/(1*mgrid + nB)
mvar  <- (12*mgrid + sB)/(1*mgrid + nB)^2
par(mfrow = c(1,2), mar = c(3,3,1.4,1.4), mgp = c(1.75,.75,0))
# media 
plot(x = mgrid, y = mmean, type = "p", pch = 18, cex = 0.8, ylim = c(9, 12), xlab = "m", ylab = "Media posterior", main = "Tipo B")
abline(h = (12 + sA)/(1 + nA), lty = 1)
# cv
plot(x = mgrid, y = sqrt(mvar)/mmean, type = "p", pch = 18, cex = 0.8, xlab = "m", ylab = "CV posterior", main = "Tipo B")

  1. Muestre que bajo el modelo Gamma-Poisson \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\), para \(i = 1,\ldots,n\), con \(\theta \sim \textsf{Gamma}(a,b)\), se tiene que la distribución predictiva posterior es Binomial Negativa con parámetros \(a+s\) y \(b+n\), donde \(s=\sum_{i=1}^n y_i\).

En este caso, se tiene que la distribución predictiva posterior está dada por: \[ \begin{align*} p(y^*\mid\boldsymbol{y}) &= \int_{\Theta} p(y^*\mid\theta)\,p(\theta\mid\boldsymbol{y})\,\text{d}\theta\\ &= \int_0^\infty \frac{e^{-\theta}\theta^{y*}}{y^*!}\,\frac{(b+n)^{a+s}}{\Gamma(a+s)}\,\theta^{a+s-1}e^{-(b+n)\theta}\,\textsf{d}\theta\\ &= \frac{(b+n)^{a+s}}{\Gamma(a+s)}\,\frac{1}{y^*!}\int_0^\infty \theta^{y^*+a+s-1}e^{-(b+n+1)\theta}\,\textsf{d}\theta\\ &= \frac{(b+n)^{a+s}}{\Gamma(a+s)}\,\frac{1}{y^*!}\,\frac{\Gamma(y^*+a+s)}{(b+n+1)^{y^*+a+s}}\\ &= \frac{\Gamma(y^*+a+s)}{\Gamma(a+s)\,\Gamma(y^*+1)}\,\frac{(b+n)^{a+s}}{(b+n+1)^{y^*+a+s}}\\ &= \frac{\Gamma(y^*+a+s)}{\Gamma(a+s)\,\Gamma(y^*+1)}\,\left(\frac{b+n}{b+n+1}\right)^{a+s}\,\left(\frac{1}{b+n+1}\right)^{y^*} \end{align*} \] donde \(\Theta=(0,\infty)\) es el espacio de parámetros y \(\boldsymbol{y}=(y_1,\ldots,y_n)\) es el vector de observaciones, \(s=\sum_{i=1}^n y_i\) un estadístico suficiente para \(\theta\), y \(n\) es el tamaño de la muestra. Por lo tanto, la distribución predictiva posterior resulta ser Binomial Negativa con parámetros \(a+s\) y \(b+n\).

  1. Considere el modelo Gamma-Poisson \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\), para \(i = 1,\ldots,n\), con \(\theta \sim \textsf{Gamma}(\alpha,\beta)\), y la función de perdida cuadrática \(L(\theta,\theta^*)=(\theta-\theta^*)^2\).

  2. Considere el modelo Gamma-Poisson \(y_i\mid\theta\stackrel{\text{iid}}{\sim}\textsf{Poisson}(\theta)\), para \(i = 1,\ldots,n\), con \(\theta \sim \textsf{Gamma}(\alpha,\beta)\), y la función de perdida cuadrática \(L(\theta,\theta^*)=(\theta-\theta^*)^2\).

  1. Muestre que el estimador que minimiza la perdida esperada posterior es de la forma \(\hat\theta = a + b\,\bar{y}\), donde \(a > 0\), \(b \in (0, 1)\) y \(\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\).

La distribución posterior es \[ \begin{align*} p(\theta\mid\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\theta)\,p(\theta)\\ &\propto \prod_{i=1}^n \frac{e^{-\theta}\theta^{y_i}}{y_i!} \times \frac{\beta^\alpha}{\Gamma(\alpha)}\,\theta^{\alpha-1}\, e^{-\beta\theta}\\ &\propto \theta^{\alpha+s -1}\,e^{-(\beta+n)\theta} \end{align*} \] donde \(s=\sum_{i=1}^n y_i\). Entonces, \(\theta\mid \boldsymbol{y}\sim \textsf{Gamma}(\alpha+s,\beta+n)\).

Si la función de pérdida es la función de pérdida cuadrática, entonces el estimador que minimiza la pérdida posterior esperada es la media posterior. Así, se tiene que \[ \hat{\theta}=\hat{\theta}(\boldsymbol{y})=\textsf{E}(\theta\mid\boldsymbol{y})=\frac{\alpha+s}{\beta+n} =\frac{\alpha}{\beta+n}+\frac{n}{\beta+n}\frac{s}{n} \] y por lo tanto, \(\hat{\theta}=a+b\bar{y}\), donde \(a=\alpha/(\beta+n)\), \(b=n/(\beta+n)\), y \(\bar{y}=\frac{s}{n}\). Como \(\alpha,\beta,n>0\), entonces \(a>0\) y \(b\in (0,1)\).

  1. Calcule el riesgo frecuentista tanto de \(\hat\theta\) como del estimador de máxima verosimilitud.

El riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) es \[ R_{\textsf{F}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}\mid\theta}(L(\theta,\hat\theta)) = \textsf{E}_{\boldsymbol{y}\mid\theta}((\theta-\hat\theta)^2) \] el cual corresponde al error cuadrático medio (MSE, por sus siglas en inglés) evaluado en \(\hat\theta\), i.e, \[ \begin{align*} R_{\textsf{F}}(\theta,\hat\theta) &= \textsf{MSE}(\hat\theta)\\ &= \textsf{Var}_{\boldsymbol{y}\mid\theta}(\hat\theta) + \left(\textsf{E}_{\boldsymbol{y}\mid\theta}(\hat\theta) - \theta\right)^2\\ &= \textsf{Var}_{\boldsymbol{y}\mid\theta}(a+b\bar{y}) + \left(\textsf{E}_{\boldsymbol{y}\mid\theta}(a+b\bar{y}) - \theta\right)^2\\ &= b^2\textsf{Var}_{\boldsymbol{y}\mid\theta}(\bar{y})+ \left(a+b\textsf{E}_{\boldsymbol{y}\mid\theta}(\bar{y}) - \theta\right)^2\\ &= b^2\frac{n\theta}{n^2}+ \left(a+b\frac{n\theta}{n} - \theta\right)^2\\ &= b^2\frac{\theta}{n}+ \left(a+(b-1)\theta \right)^2\\ &= (b-1)^2\theta^2 + \left(\frac{b^2}{n}+2a(b-1)\right)\theta + a^2\\ &=c_1\theta^2 + c_2\theta + c_3 \end{align*} \] donde \(c_1=(b-1)^2\), \(c_2= \frac{b^2}{n}+2a(b-1)\), y \(c_3=a^2\). Después de algo de álgebra, se obtiene que \[ \begin{align*} &c_1=(b-1)^2=\left(\frac{n}{\beta + n}-1\right)^2\quad\Rightarrow\quad c_1=\frac{\beta^2}{(\beta+n)^2}\\ &c_2=\frac{b^2}{n}+2a(b-1)=\frac{\left(\frac{n}{\beta+n}\right)^2}{n}+2\frac{\alpha}{\beta+n}\left(\frac{n}{\beta+n}-1\right) \quad\Rightarrow\quad c_2=\frac{n-2\alpha\beta}{(\beta+n)^2}\\ &c_3=\left(\frac{\alpha}{\beta+n}\right)^2\quad\Rightarrow\quad c_3=\frac{\alpha^2}{(\beta+n)^2} \end{align*} \]

De otra parte, se demuestra que el estimador de máxima verosimilitud (MLE, por sus siglas en inglés) es \(\hat{ \theta }_{ \text{MLE} } = \bar{y}\) (¡ejercicio!).

Así, se tiene que le riesgo frecuentista de \(\hat{ \theta }_{ \text{MLE} }\) es \[ \begin{align*} R_{\textsf{F}}(\theta,\hat{ \theta }_{ \text{MLE} }) &= \textsf{MSE}(\hat{ \theta }_{ \text{MLE} })\\ &= \textsf{Var}_{\boldsymbol{y}\mid\theta}(\hat{ \theta }_{ \text{MLE} }) + \left(\textsf{E}_{\boldsymbol{y}\mid\theta}(\hat{ \theta }_{ \text{MLE} }) - \theta\right)^2\\ &= \textsf{Var}_{\boldsymbol{y}\mid\theta}(\bar{y}) + \left(\textsf{E}_{\boldsymbol{y}\mid\theta}(\bar{y}) - \theta\right)^2\\ &= \frac{n\theta}{n^2}+ \left(\frac{n\theta}{n} - \theta\right)^2\\ &= \frac{\theta}{n} \end{align*} \]

  1. Calcule el riesgo Bayesiano de \(\hat\theta\).

El riesgo Bayesiano de \(\hat\theta\) es \[ \begin{align*} R_{\textsf{B}}(\theta,\hat\theta) &= \textsf{E}_{\theta}(R_{\textsf{F}}(\theta,\hat\theta))\\ &= \int_{\Theta} R_{\textsf{F}}(\theta,\hat\theta)\, p(\theta)\,\textsf{d}\theta\,,\\ &= \int_{\Theta} (c_1\theta^2 + c_2\theta + c_3)\,p(\theta) \, \textsf{d}\theta\\ &= c_1 \int_{\Theta} \theta^2\, p(\theta) \, \textsf{d}\theta + c_2 \int_{\Theta} \theta \, p(\theta) \, \textsf{d}\theta + c_3 \int_{\Theta} p(\theta) \, \textsf{d}\theta\\ &=c_1\textsf{E}_\theta(\theta^2) + c_2\textsf{E}_\theta(\theta)+c_3\\ &=c_1\left(\textsf{Var}_\theta(\theta)+(\textsf{E}_\theta(\theta))^2\right) + c_2\textsf{E}_\theta(\theta)+c_3\\ &=c_1\left(\frac{\alpha}{\beta^2}+\frac{\alpha^2}{\beta^2}\right) + c_2\frac{\alpha}{\beta}+c_3\\ &=c_1\frac{\alpha(1-\alpha)}{\beta^2}+c_2\frac{\alpha}{\beta}+c_3 \end{align*} \] donde \[ \begin{align*} c_1=\frac{\beta^2}{(\beta+n)^2},\quad c_2=\frac{n-2\alpha\beta}{(\beta+n)^2},\quad\text{y}\quad c_3=\frac{\alpha^2}{(\beta+n)^2}\,. \end{align*} \]

  1. Suponga que un investigador quiere recolectar una muestra lo suficientemente grande para que el riesgo Bayesiano después del experimento sea la mitad del riesgo Bayesiano antes del experimento. Encuentre ese tamaño de muestra.

Se necesita un tamaño de muestra tal que \[ \textsf{E}_{\theta\mid\boldsymbol{y}}(R_{\textsf{F}}(\theta,\hat\theta))=\frac12\,\textsf{E}_{\theta}(R_{\textsf{F}}(\theta,\hat\theta)). \] Ya se sabe que el riesgo Bayesiano antes del experimento es de la forma \(\textsf{E}_{\theta}(R_{\textsf{F}}(\theta,\hat\theta))=c_1\frac{\alpha(1-\alpha)}{\beta^2}+c_2\frac{\alpha}{\beta}+c_3\).

Siguiendo la misma idea, se tiene que el riesgo Bayesiano de \(\hat\theta\) después del experimento es \[ \begin{align*} \textsf{E}_{\theta\mid\boldsymbol{y}}(R_{\textsf{F}}(\theta,\hat\theta)) &=\int_\Theta R_{\textsf{F}}(\theta,\hat\theta)\,p(\theta\mid\boldsymbol{y})\, \textsf{d}\theta\\ &= \int_\Theta (c_1\theta^2 + c_2\theta + c_3) p(\theta\mid\boldsymbol{y}) \, \textsf{d}\theta\\ &= c_1 \int_\Theta \theta^2\, p(\theta\mid\boldsymbol{y}) \, \textsf{d}\theta + c_2 \int_\Theta \theta \, p(\theta\mid\boldsymbol{y}) \, \textsf{d}\theta + c_3 \int_\Theta \pi(\theta\mid\boldsymbol{y}) \, d\theta\\ &=c_1\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta^2) + c_2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)+c_3\\ &=c_1\left(\textsf{Var}_{\theta\mid\boldsymbol{y}}(\theta) + (\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta))^2\right) + c_2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)+c_3\\ &=c_1\left(\frac{\alpha+s }{(\beta + n)^2}+\frac{\left(\alpha+s\right)^2}{(\beta+n)^2}\right) + c_2\frac{\alpha+s}{\beta+n}+c_3\\ &=c_1\frac{\left(\alpha+s\right)\left(1-\left(\alpha+s\right)\right)}{(\beta + n)^2}+c_2\frac{\alpha}{\beta + n}+c_3 \end{align*} \]

donde \[ \begin{align*} c_1=\frac{\beta^2}{(\beta+n)^2},\quad c_2=\frac{n-2\alpha\beta}{(\beta+n)^2},\quad\text{y}\quad c_3=\frac{\alpha^2}{(\beta+n)^2}\,. \end{align*} \]

Entonces, se requiere el tamaño de la muestra \(n\) que satisfaga la ecuación \[ 2c_1\frac{\left(\alpha+s\right)\left(1-\left(\alpha+s\right)\right)}{(\beta + n)^2}+2c_2\frac{\alpha}{\beta + n}+2c_3 = c^*_1\frac{\alpha(1-\alpha)}{\beta^2}+c^*_2\frac{\alpha}{\beta}+c^*_3 \]

donde \[ \begin{align*} c^*_1=\frac{\beta^2}{(\beta+0)^2}=1,\quad c^*_2=\frac{0-2\alpha\beta}{(\beta+0)^2}=-2\frac{\alpha}{\beta},\quad\text{y}\quad c^*_3=\frac{\alpha^2}{(\beta+0)^2}=\frac{\alpha^2}{\beta^2}\,. \end{align*} \]

Se observa que los hiperparámetros \(\alpha\) y \(\beta\) son cantidades conocidas, y por lo tanto, las constantes \(c_k\) y \(c_k^*\), para \(k=1,2,3\), también lo son. Así mismo, el estadístico suficiente \(s=\sum_{i=1}^n y_i\) también es una cantidad conocida después de que el experimento tenga lugar. En consecuencia, en la ecuación anterior la única cantidad desconocida es el tamaño de muestra \(n\).

  1. Considere el modelo Beta-Binomial \(x\mid\theta\sim \textsf{Bin}(n,\theta)\) con \(n\) conocido, y \(\theta\sim \textsf{Beta}(\sqrt{n}/2,\sqrt{n}/2)\).
  1. ¿Cuál es la distribución posterior?

La distribución posterior es: \[ \begin{align*} p(\theta\mid x) &= \frac{p(x\mid\theta)\,p(\theta)}{\int_\Theta p(x\mid\theta)\,p(\theta)\, \textsf{d}\theta}\\ &= \frac{\binom{n}{x} \theta^x (1-\theta)^{n-x}\frac{1}{B(\sqrt{n}/2,\sqrt{n}/2)} \theta^{\sqrt{n}/2-1} (1-\theta)^{\sqrt{n}/2-1}}{\int_\Theta \binom{n}{x} \theta^x (1-\theta)^{n-x}\frac{1}{B(\sqrt{n}/2,\sqrt{n}/2)} \theta^{\sqrt{n}/2-1} (1-\theta)^{\sqrt{n}/2-1} \, \textsf{d}\theta}\\ &= \frac{\theta^{\sqrt{n}/2 + x-1} (1-\theta)^{\sqrt{n}/2 + n-x-1}}{\int_\Theta \theta^{\sqrt{n}/2 + x-1} (1-\theta)^{\sqrt{n}/2 + n-x-1} \, \textsf{d}\theta}\\ &= \frac{\theta^{\sqrt{n}/2 + x-1} (1-\theta)^{n-x+\sqrt{n}/2-1}}{\textsf{B}(\sqrt{n}/2 + x,\sqrt{n}/2 + n-x)}\\ \end{align*} \] y por lo tanto, \(\theta\mid x \sim \textsf{Beta}(\sqrt{n}/2 + x,\sqrt{n}/2 + n-x)\).

  1. ¿Cuál es el estimador que minimiza la perdida esperada posterior si la función de perdida es \(L(\theta,\hat\theta)=(\theta-\hat\theta)^2\)? Llame tal estimador \(\hat{\theta}\) y muestre que el riesgo frecuentista correspondiente es constante.

Si la función de perdida es \(L(\theta,\hat\theta)=(\theta-\hat\theta)^2\), i.e., la función de perdida cuadrática, entonces el estimador que minimiza la función de perdida esperada a posteriori es la media posterior. Así, se tiene que \[ \hat{\theta}=\textsf{E}(\theta\mid x)=\frac{x+\sqrt{n}/2}{(x+\sqrt{n}/2)+(n-x+\sqrt{n}/2)} = \frac{x+\sqrt{n}/2}{\sqrt{n}+ n} = \frac{\frac{x}{\sqrt{n}}+\frac12}{1 + \sqrt{n}}\,. \] Ahora, si la función de perdida es la función de perdida cuadrática, entonces el riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) es \[ R_{\textsf{F}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}\mid\theta}(L(\theta,\hat\theta)) = \textsf{E}_{\boldsymbol{y}\mid\theta}((\theta-\hat\theta)^2) \] el cual corresponde al error cuadrático medio (MSE, por sus siglas en inglés) evaluado en \(\hat\theta\), i.e, \[ \begin{align*} R_{\textsf{F}}(\theta,\hat\theta) &= \textsf{MSE}(\hat\theta)\\ &= \textsf{Var}_{x\mid\theta}(\hat\theta) + \left(\textsf{E}_{x\mid\theta}(\hat\theta) - \theta\right)^2\\ &= \textsf{Var}_{x\mid\theta}\left(\frac{\frac{x}{\sqrt{n}}+\frac12}{1 + \sqrt{n}}\right) + \left(\textsf{E}_{x\mid\theta}\left(\frac{\frac{x}{\sqrt{n}}+\frac12}{1 + \sqrt{n}}\right) - \theta\right)^2\\ &= \left(\frac{\frac{1}{\sqrt{n}}}{1 + \sqrt{n}}\right)^2\textsf{Var}_{x\mid\theta}(x)+ \left(\frac{\frac{\textsf{E}_{x\mid\theta}(x)}{\sqrt{n}}+\frac12}{1 + \sqrt{n}} - \theta\right)^2\\ &= \frac{1}{n(1 + \sqrt{n})^2}\,n\theta(1-\theta) + \left(\frac{\frac{n\theta}{\sqrt{n}}+\frac12}{1 + \sqrt{n}}-\theta\right)^2\\ &=\frac{\theta(1-\theta)}{(1 + \sqrt{n})^2} + \left(\frac{\sqrt{n}\theta+\frac12 -\theta - \sqrt{n}\theta}{1 + \sqrt{n}}\right)^2\\ &=\frac{\theta(1-\theta)}{(1 + \sqrt{n})^2} + \frac{(\frac12 -\theta)^2}{(1 + \sqrt{n})^2}\\ &=\frac{\theta -\theta^2+\frac14 -\theta + \theta^2}{(1 + \sqrt{n})^2}\\ &=\frac{1}{4(1 + \sqrt{n})^2}\\ \end{align*} \] lo cual resulta ser constante.

  1. Sea \(\theta_0=x/n\). Encuentre el riesgo frecuentista de este estimador. Compare los riesgos de \(\hat{\theta}\) y \(\theta_0\), para \(n=10,50,100\). ¿Qué se puede concluir?

El riesgo de \(\theta_0=x/n\) es \[ \begin{align*} R_{\textsf{F}}(\theta,\theta_0)&= \textsf{MSE}(\theta_0)\\ &= \textsf{Var}_{x\mid\theta}(\theta_0) + \left(\textsf{E}_{x\mid\theta}(\theta_0) - \theta\right)^2\\ &= \textsf{Var}_{x\mid\theta}(x/n) + \left(\textsf{E}_{x\mid\theta}(x/n) - \theta\right)^2\\ &= \frac{1}{n^2}\textsf{Var}_{x\mid\theta}(x)+ \left(\frac{\textsf{E}_{x\mid\theta}(x)}{n}-\theta\right)^2\\ &= \frac{n\theta(1-\theta)}{n^2}+ \left[\frac{n\theta}{n}-\theta\right]^2\\ &= \frac{\theta(1-\theta)}{n}\\ \end{align*} \]

Los gráficos de los riesgos de \(\hat\theta\) y \(\theta_0\) para \(n=10,50,100\) muestran que ninguno de los estimadores domina al otro para todos los valores de \(\theta\in\Theta\).

Se observa que \(\hat\theta\) domina \(\theta_0\) cuando:

Un patrón claro consiste en que a medida que le tamaño de muestra se incrementa:

# r.1 = riesgo de theta^
# r.12= riesgo de theta0
r.1  <- function(theta, n) theta^0/(4*(1+sqrt(n))^2)
r.2  <- function(theta, n) theta*(1-theta)/n
grid <- seq(from = 0.0, to = 1.0, length = 1000)
# n = 10
fun <- function(x) r.1(x, n=10) - r.2(x, n=10)
uni1 <- uniroot(fun, c(0.0, 0.5))$root
uni2 <- uniroot(fun, c(0.5, 1.0))$root
plot(grid, r.2(theta=grid, n=10), type = 'l', cex.axis=0.7, xlab = expression(theta), ylab = 'Riesgo', lwd = 1, las=1, col="red", main="n=10", ylim=c(0,0.025))
lines(grid, r.1(theta=grid, n=10), lwd = 1, col="blue")
legend("bottom", c(expression(hat(theta)),expression(theta[0])), col = c("blue", "red"), lty=1, lwd=2, bty="n") 
abline(v = uni1, lty = 3)
abline(v = uni2, lty = 3)
points(uni1, r.1(uni1, n=10), pch = 16, cex = 0.8)
points(uni2, r.1(uni2, n=10), pch = 16, cex = 0.8)
text( 0.3, 0.013, '(0.174;0.014)', cex = 1.0 )
text( 0.7, 0.013, '(0.825;0.014)', cex = 1.0 )
segments(uni1, y0=0.01, x1 = uni2, lwd=3)
text( 0.5, 0.008, '0.650', cex = 1.0)

# n = 50
fun <- function(x) r.1(x, n=50) - r.2(x, n=50)
uni1 <- uniroot(fun, c(0.0, 0.5))$root
uni2 <- uniroot(fun, c(0.5, 1.0))$root
plot(grid, r.2(theta=grid, n=50), type = 'l', cex.axis=0.7, xlab = expression(theta), ylab = 'Riesgo', lwd = 1, las=1, col="red", main="n=50", ylim=c(0,0.025))
lines(grid, r.1(theta=grid, n=50), lwd = 1, col="blue")
legend("top", c(expression(hat(theta)),expression(theta[0])), col = c("blue", "red"), lty=1, lwd=2, bty="n")    
abline(v = uni1, lty = 3)
abline(v = uni2, lty = 3)
points(uni1, r.1(uni1, n=50), pch = 16, cex = 0.8)
points(uni2, r.1(uni2, n=50), pch = 16, cex = 0.8)
text( 0.1, 0.006, '(0.258;0.003)', cex = 1.0 )
text( 0.9, 0.006, '(0.741;0.003)', cex = 1.0 )
segments(uni1, y0=0.01, x1 = uni2, lwd=3)
text( 0.5, 0.008, '0.482', cex = 1.0)

# n = 100
fun <- function(x) r.1(x, n=100) - r.2(x, n=100)
uni1 <- uniroot(fun, c(0.0, 0.5))$root
uni2 <- uniroot(fun, c(0.5, 1.0))$root
plot(grid, r.2(theta=grid, n=100), type = 'l', cex.axis=0.7, xlab = expression(theta), ylab = 'Riesgo', lwd = 1, las=1, col="red", main="n=100", ylim=c(0,0.025))
lines(grid, r.1(theta=grid, n=100), lwd = 1, col="blue")
legend("top",  c(expression(hat(theta)),expression(theta[0])), col = c("blue", "red"), lty=1, lwd=2, bty="n")   
abline(v = uni1, lty = 3)
abline(v = uni2, lty = 3)
points(uni1, r.1(uni1, n=100), pch = 16, cex = 0.8)
points(uni2, r.1(uni2, n=100), pch = 16, cex = 0.8)
text( 0.1, 0.004, '(0.291;0.002)', cex = 1.0 )
text( 0.9, 0.004, '(0.708;0.002)', cex = 1.0 )
segments(uni1, y0=0.01, x1 = uni2, lwd=3)
text( 0.5, 0.008, '0.416', cex = 1.0)

  1. Sea \(x\sim \textsf{Bin}(n,\theta)\) con \(n\) conocido, y \(\theta \sim \textsf{Beta}(\alpha,\beta)\). Considere la función de perdida \(L(\theta,\hat\theta)=(\theta-\hat\theta)/(\theta(1-\theta))\). Encuentre el estimador que minimiza la perdida esperada posterior para esta función de perdida.

En este caso la distribución posterior es \(\theta\mid x\sim \textsf{Beta}(\alpha+x,\beta+n-x)\). Entonces, se tiene que la perdida esperada posterior es \[ \begin{align*} \textsf{E}_{\theta\mid x}(L(\theta,\hat\theta)) &= \int_{\Theta} L(\theta,\hat\theta)\, p(\theta\mid x)\, \textsf{d}\theta\\\ &= \int_\Theta \frac{(\theta-\hat\theta)^2}{\theta(1-\theta)} \frac{\theta^{\alpha+x-1}(1-\theta)^{\beta+n-x-1}}{B(\alpha+x,\beta+n-x)} \,\textsf{d}\theta\\ &= \frac{1}{\textsf{B}(\alpha+x,\beta+n-x)}\int_\Theta (\theta-\hat\theta)^2 \theta^{\alpha+x-2}(1-\theta)^{\beta+n-x-2} \,\textsf{d}\theta\\ &= \frac{\textsf{B}(\alpha+x-1,\beta+n-x-1)}{\textsf{B}(\alpha+x,\beta+n-x)}\int_\Theta (\theta-\hat\theta)^2 \frac{\theta^{(\alpha+x-1)-1}(1-\theta)^{(\beta+n-x-1)-1}}{B(\alpha+x-1,\beta+n-x-1)} \,\textsf{d}\theta\\ &= \frac{\textsf{B}(\alpha+x-1,\beta+n-x-1)}{\textsf{B}(\alpha+x,\beta+n-x)}\int_\Theta L^*(\theta,\hat\theta)\, p^*(\theta\mid x) \,\textsf{d}\theta \end{align*} \] donde \[ p^*(\theta\mid x)=\frac{\theta^{(\alpha+x-1)-1}(1-\theta)^{(\beta+n-x-1)-1}}{\textsf{B}(\alpha+x-1,\beta+n-x-1)}\qquad\text{y}\qquad L^*(\theta,\hat\theta)=(\theta-\hat\theta)^2. \] Dado que \[ \int_\Theta L(\theta,\hat\theta)\, p(\theta\mid x)\, \textsf{d}\theta = c\,\int_\Theta L^*(\theta,\hat\theta)\, p^*(\theta\mid x) \,\textsf{d}\theta \] donde \(c\) es una constante que no depende de \(\theta\) dada por \[ c=\frac{\textsf{B}(\alpha+x-1,\beta+n-x-1)}{\textsf{B}(\alpha+x,\beta+n-x)}\,, \] \(L^*(\theta,\hat\theta)\) es la función de perdida cuadrática, entonces el estimador que minimiza la perdida esperada posterior de acuerdo con \(p^*(\theta\mid x)\), y en consecuencia con \(p(\theta\mid x)\), es la media posterior de \(p^*(\theta\mid x)\). Así se tiene que \[ \hat\theta=\textsf{E}(\theta\mid x)=\frac{\alpha+x-1}{(\alpha+x-1)+(\beta+n-x-1)}=\frac{\alpha+x-1}{\alpha+\beta+n-2}\,. \]

  1. Sea \(L(\theta,\hat\theta) = \omega(\theta)(\theta - \hat\theta)^2\) la función de perdida cuadrática ponderada, donde \(\omega(\theta)\) es una función no negativa. Muestre que el estimador que minimiza la perdida esperada posterior tiene la forma \[ \hat\theta=\frac{\textsf{E}(\omega(\theta)\,\theta\mid \boldsymbol{y})}{\textsf{E}(\omega(\theta)\mid \boldsymbol{y})}\,. \]

Se tiene que \[ \begin{align*} \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)(\theta - \hat\theta)^2) &= \int_\Theta \omega(\theta)(\theta-\hat\theta)^2\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta\\ &=\int_\Theta\omega(\theta)(\theta-\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)+\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)^2\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta\\ &=\int_\Theta\omega(\theta)(\theta-\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta))^2\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta +2\int_\Theta\omega(\theta)(\theta-\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta))(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta +\int_\Theta\omega(\theta)(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)^2\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta \end{align*} \]

La primera integral es constante respecto a \(\hat\theta\) y no afecta el proceso de minimización de respecto a \(\hat\theta\). De otra parte, la segunda integral es \[ \begin{align*} 2(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)\int_\Theta\omega(\theta)(\theta-\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta))\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta &=2(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)\left( \int_\Theta\omega(\theta)\,\theta\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta -\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\int_\Theta\omega(\theta)\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta \right)\\ &=2(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)\left( \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) -\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\,\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) \right) \end{align*} \] De esta forma, tomando la primera derivada de \(\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)(\theta - \hat\theta)^2)\) respecto a \(\hat\theta\), se tiene que \[ \begin{align*} \frac{\textsf{d}}{\textsf{d}\hat\theta}\,\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)(\theta - \hat\theta)^2) &=-2\left( \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) - \textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) \right) +\frac{\textsf{d}}{\textsf{d}\hat\theta}\int_\Theta\omega(\theta)(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)^2\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta\\ &=-2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) +\frac{\textsf{d}}{\textsf{d}\hat\theta}(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)^2\int_\Theta\omega(\theta)\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta\\ &=-2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) -2(\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)-\hat\theta)\int_\Theta\omega(\theta)\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta\\ &=-2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) -2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\int_\Theta\omega(\theta)\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta+2\hat\theta\int_\Theta\omega(\theta)\,p(\theta\mid \boldsymbol{y})\,\textsf{d}\theta\\ &=-2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) -2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta))+2\hat\theta\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta))\\ &=-2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\hat\theta\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) \end{align*} \]

Igualando esta derivada a cero y despejando para \(\hat\theta\) con el fin de hallar el punto crítico correspondiente, se tiene que \[ -2 \textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta) +2\hat\theta\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) = 0 \] y en consecuencia, de \(\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)(\theta - \hat\theta)^2)\) respecto a \(\hat\theta\) es \[ \hat\theta = \frac{\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)\,\theta)}{\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta))} \]

Además, la segunda derivada \[ \frac{\textsf{d}^2}{\textsf{d}\hat\theta^2}\,\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)(\theta - \hat\theta)^2)=2\textsf{E}_{\theta\mid\boldsymbol{y}}(\omega(\theta)) >0 \] para todo \(\hat\theta\) (pues \(\omega(\theta)\) es una función no negativa), lo que significa que el punto crítico correspondiente en efecto minimiza la función objetivo.

  1. Un investigador del Departamento de Ingeniería Electrónica y Eléctrica de la Universidad de Bath (Inglaterra) necesitaba analizar unos datos sobre los tiempos de falla de un determinado tipo de alambre de metal (en este problema el tiempo de falla se define como el número de veces que una máquina podría tensionar el alambre antes de romperse). Los siguientes datos corresponden a \(n = 14\) tiempos de falla de una parte del experimento: \[ \boldsymbol{y} = (495,541,1461,1555,1603,2201,2750,3468,3516,4319,6622,7728,13159,21194)\,. \] El modelo más simple para los datos del tiempo de falla involucra la distribución Exponencial, \(y_i \mid \lambda \stackrel{\text{iid}}{\sim} \textsf{Exponencial} ( \lambda )\), con \(\textsf{E}(y_i \mid \lambda) = \lambda\), para \(i=1,\ldots,n\).
  1. Sea \(X \sim \textsf{Gamma} (\alpha, \beta)\), con \(\alpha > 0\) y \(\beta > 0\). Encuentre la función de densidad de probabilidad de \(Y=\frac{ 1 }{ X }\). La distribución de \(Y\) se denomina distribución Gamma Inversa con parámetros \(\alpha\) y \(\beta\), lo que se denota con \(Y\sim\textsf{GI} (\alpha,\beta)\).

Para encontrar la densidad de \(Y = \frac{1}{X}\), donde \(X \sim \textsf{Gamma}(\alpha, \beta)\), seguimos los pasos del cambio de variable para distribuciones de probabilidad.

La función de densidad de \(X \sim \textsf{Gamma}(\alpha, \beta)\) es: \[ f_X(x) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}, \quad x > 0. \]

Sea \(Y = \frac{1}{X}\). Entonces, \(X = \frac{1}{Y}\) y la transformación es válida para \(Y > 0\). Ahora, derivamos \(X\) con respecto a \(Y\) para encontrar el jacobiano del cambio de variable: \[ \frac{\text{d}X}{\text{d}Y} = -\frac{1}{Y^2}. \] El valor absoluto del jacobiano es: \[ \left| \frac{\text{d}X}{\text{d}Y} \right| = \frac{1}{Y^2}. \]

La densidad de \(Y\) se obtiene usando la fórmula del cambio de variable: \[ f_Y(y) = f_X(x) \left| \frac{\text{d}X}{\text{d}Y} \right|, \] donde \(x = \frac{1}{y}\). Sustituyendo: \[ f_Y(y) = f_X\left(\frac{1}{y}\right) \cdot \frac{1}{y^2}. \]

Sustituyendo la densidad de \(X\): \[ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} \left(\frac{1}{y}\right)^{\alpha - 1} e^{-\beta \frac{1}{y}} \cdot \frac{1}{y^2}. \]

Simplificamos los términos: \[ \left(\frac{1}{y}\right)^{\alpha - 1} = y^{-(\alpha - 1)}, \quad e^{-\beta \frac{1}{y}} \text{ permanece igual}, \quad \frac{1}{y^2} = y^{-2}. \]

Por tanto: \[ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{-(\alpha - 1)} e^{-\beta \frac{1}{y}} y^{-2}. \]

Agrupando los exponentes de \(y\): \[ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{-(\alpha + 1)} e^{-\beta \frac{1}{y}}, \quad y > 0. \]

Esta es la densidad de una distribución Gamma Inversa (\(\textsf{GI}(\alpha, \beta)\)), cuya forma general es: \[ f_Y(y) = \frac{\beta^\alpha}{\Gamma(\alpha)} y^{-(\alpha + 1)} e^{-\frac{\beta}{y}}, \quad y > 0. \]

Por lo tanto, hemos mostrado que si \(X \sim \textsf{Gamma}(\alpha, \beta)\), entonces \(Y = \frac{1}{X} \sim \textsf{GI}(\alpha, \beta)\).

  1. Encuentre la distribución posterior de \(\lambda\), asumiendo que \(\lambda\sim\textsf{GI}(a,b)\), con \(a > 0\) y \(b > 0\).

En este problema, la distribución de los datos está modelada como \(y_i \mid \lambda \stackrel{\text{iid}}{\sim} \textsf{Exponencial}(\lambda)\) para \(i = 1, \ldots, n\), mientras que el parámetro \(\lambda\) sigue una distribución a priori \(\textsf{GI}(a, b)\), con \(a > 0\) y \(b > 0\). La función de densidad de \(\lambda \sim \textsf{GI}(a, b)\) está dada por: \[ p(\lambda) = \frac{b^a}{\Gamma(a)} \lambda^{-(a+1)} e^{-\frac{b}{\lambda}}, \quad \lambda > 0. \] La verosimilitud de los datos \(\boldsymbol{y} = (y_1, \ldots, y_n)\) bajo la distribución exponencial es: \[ p(\boldsymbol{y} \mid \lambda) = \prod_{i=1}^n \lambda e^{-\lambda y_i} = \lambda^n e^{-\lambda \sum_{i=1}^n y_i}. \] Denotando por \(s = \sum_{i=1}^n y_i\), la verosimilitud toma la forma compacta \(p(\boldsymbol{y} \mid \lambda) = \lambda^n e^{-\lambda s}\).

La distribución posterior de \(\lambda\) es proporcional al producto de la verosimilitud y la densidad a priori, es decir: \[ p(\lambda \mid \boldsymbol{y}) \propto p(\boldsymbol{y} \mid \lambda) p(\lambda). \] Sustituyendo las expresiones de la verosimilitud y la densidad a priori: \[ p(\lambda \mid \boldsymbol{y}) \propto \lambda^n e^{-\lambda s} \cdot \lambda^{-(a+1)} e^{-\frac{b}{\lambda}} = \lambda^{n - a - 1} e^{-\lambda s} e^{-\frac{b}{\lambda}}. \] Reconociendo que esta expresión coincide con la densidad de una distribución Gamma Inversa, concluimos que la distribución posterior de \(\lambda\) es \(\textsf{GI}(a^*, b^*)\), donde los parámetros actualizados son \(a^* = a + n\) y \(b^* = b + s\).

Por lo tanto, la distribución posterior de \(\lambda\) es: \[ \lambda \mid \boldsymbol{y} \sim \textsf{GI}(a + n, b + s), \] donde \(a + n\) refleja la actualización de la información sobre el parámetro \(\lambda\) basada en el tamaño de la muestra \(n\), y \(b + s\) incorpora la suma total de los tiempos de falla \(s = \sum_{i=1}^n y_i\). Esta forma de la distribución posterior muestra cómo los datos se combinan con la información previa para actualizar nuestro conocimiento sobre \(\lambda\).

  1. Muestre que la media posterior de \(\lambda\) es un promedio ponderado de la media previa \(\textsf{E}(\lambda)\) y la media muestral \(\bar{y}=\frac{1}{n}\sum_{i=1}^n y_i\).

Para demostrar que la media posterior de \(\lambda\) es un promedio ponderado de la media previa \(\textsf{E}(\lambda)\) y la media muestral \(\bar{y} = \frac{1}{n}\sum_{i=1}^n y_i\), partimos del hecho de que, dada la distribución posterior \(\lambda \mid \boldsymbol{y} \sim \textsf{GI}(a^*, b^*)\), la media de una distribución Gamma Inversa está dada por: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \frac{b^*}{a^* - 1}, \quad \text{para } a^* > 1. \] En este caso, los parámetros de la posterior son \(a^* = a + n\) y \(b^* = b + s\), donde \(s = \sum_{i=1}^n y_i\). Sustituyendo estos valores, la media posterior de \(\lambda\) se convierte en: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \frac{b + s}{a + n - 1}. \]

Ahora, analizamos la relación entre la media posterior y las medias previa y muestral. La media previa de \(\lambda\), dada su distribución \(\textsf{GI}(a, b)\), es: \[ \textsf{E}(\lambda) = \frac{b}{a - 1}, \quad \text{para } a > 1. \] La media muestral de los datos es \(\bar{y} = \frac{s}{n}\). Sustituyendo \(s = n\bar{y}\) en la expresión de la media posterior, obtenemos: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \frac{b + n\bar{y}}{a + n - 1}. \]

Podemos reescribir esta expresión como un promedio ponderado de \(\textsf{E}(\lambda)\) y \(\bar{y}\). Multiplicamos y dividimos cada término en el numerador por su denominador respectivo para separar los términos: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \frac{b}{a - 1} \cdot \frac{a - 1}{a + n - 1} + \bar{y} \cdot \frac{n}{a + n - 1}. \] Aquí, los coeficientes \(\omega = \frac{a - 1}{a + n - 1}\) y \(1 - \omega = \frac{n}{a + n - 1}\) representan los pesos de la media previa y la media muestral, respectivamente, y satisfacen \(\omega + (1 - \omega) = 1\).

Por lo tanto, la media posterior puede escribirse como: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \omega \, \textsf{E}(\lambda) + (1 - \omega) \, \bar{y}, \] donde los pesos \(\omega\) y \(1 - \omega\) dependen de los parámetros \(a\), \(n\), y el tamaño de la muestra. Esto demuestra que la media posterior de \(\lambda\) es un promedio ponderado de la media previa y la media muestral.

  1. Se tiene información externa de otro experimento de acuerdo con el cual la distribución previa de \(\lambda\) debería tener una media \(\mu_0 = 4500\) y una desviación estándar \(\sigma_0 = 1800\). Grafique las distribuciones previa y posterior de \(\lambda\) en el mismo gráfico con \(\lambda\) en el eje horizontal tomando valores de 1000 a 12000, identificando qué curva corresponde a qué densidad.

La información externa sobre la distribución previa de \(\lambda\) indica que su media es \(\mu_0 = 4500\) y su desviación estándar es \(\sigma_0 = 1800\). Usando la relación entre la media y la desviación estándar de la distribución Gamma Inversa, podemos determinar los parámetros \(a\) y \(b\) de la distribución previa.

Para una distribución \(\lambda \sim \textsf{GI}(a, b)\), la media y la desviación estándar están dadas por: \[ \textsf{E}(\lambda) = \frac{b}{a - 1}, \quad \textsf{SD}(\lambda) = \sqrt{\frac{b^2}{(a-1)^2(a-2)}} \quad \text{para } a > 2. \] Sustituyendo \(\textsf{E}(\lambda) = \mu_0 = 4500\) y \(\textsf{SD}(\lambda) = \sigma_0 = 1800\), obtenemos: \[ \frac{b}{a - 1} = 4500, \quad \sqrt{\frac{b^2}{(a-1)^2(a-2)}} = 1800. \]

Resolviendo el sistema, primero reescribimos \(b\) en términos de \(a\) usando la ecuación de la media: \[ b = 4500(a - 1). \] Sustituyendo en la ecuación de la desviación estándar: \[ \sqrt{\frac{[4500(a - 1)]^2}{(a-1)^2(a-2)}} = 1800. \] Simplificando: \[ \sqrt{\frac{4500^2}{a-2}} = 1800. \] Elevando al cuadrado y resolviendo para \(a\): \[ \frac{4500^2}{a-2} = 1800^2 \implies a - 2 = \frac{4500^2}{1800^2} = 6.25 \implies a = 8.25. \] Usando \(b = 4500(a - 1)\), calculamos: \[ b = 4500(8.25 - 1) = 4500 \cdot 7.25 = 32625. \]

Por lo tanto, la distribución previa de \(\lambda\) es \(\textsf{GI}(8.25, 32625)\).

La distribución posterior se calculó previamente como \(\textsf{GI}(a^*, b^*)\), donde: \[ a^* = a + n, \quad b^* = b + s. \] Sustituyendo \(a = 8.25\), \(b = 32625\), \(n = 14\), y \(s = \sum_{i=1}^{14} y_i = 71208\), obtenemos: \[ a^* = 8.25 + 14 = 22.25, \quad b^* = 32625 + 71208 = 103833. \] Por lo tanto, la distribución posterior de \(\lambda\) es \(\textsf{GI}(22.25, 103833)\).

Para graficar las distribuciones previa y posterior, usamos sus densidades: \[ p(\lambda) = \frac{b^a}{\Gamma(a)} \lambda^{-(a+1)} e^{-\frac{b}{\lambda}}, \quad \lambda > 0. \]

# Parámetros de la distribución previa
a_prior <- 8.25
b_prior <- 32625
# Parámetros de la distribución posterior
a_post <- 22.25
b_post <- 103833
# Rango de lambda
lambda <- seq(1000, 12000, length.out = 1000)
# Densidad de la distribución previa
prior_density <- (b_prior^a_prior / gamma(a_prior)) * lambda^(-(a_prior + 1)) * exp(-b_prior / lambda)
# Densidad de la distribución posterior
post_density <- (b_post^a_post / gamma(a_post)) * lambda^(-(a_post + 1)) * exp(-b_post / lambda)
# Gráfico
plot(lambda, prior_density, type = "l", col = "blue", lwd = 2, ylim = c(0,0.00045),
     ylab = "Densidad", xlab = expression(lambda),
     main = "Distribuciones Previa y Posterior de Lambda")
lines(lambda, post_density, col = "red", lwd = 2, lty = 2)
legend("topright", legend = c("Previa", "Posterior"),
       col = c("blue", "red"), lwd = 2, lty = c(1, 2), bty = "n")

En el gráfico, la distribución previa, representada en azul, refleja la información externa inicial sobre \(\lambda\), con una media de 4500 y una desviación estándar de 1800. La distribución posterior, en rojo, actualiza esta información a partir de los datos observados, reflejando tanto la información previa como los tiempos de falla del experimento. La posterior es más concentrada debido a la mayor cantidad de información aportada por los datos.

  1. Muestre que el estimador de máxima verosimilitud (MLE, por sus siglas en inglés) de \(\lambda\) y la información observada de Fisher son respectivamente \[ \hat\lambda_{\text{MLE}} = \bar{y} \qquad\text{y}\qquad \hat{I}(\hat\lambda_{\text{MLE}}) = \frac{n}{\bar{y}^2}\,. \]

Recuerde que la información observada de Fisher se define como \[ \hat{I}(\hat\lambda_{\text{MLE}}) = \left[ -\frac{\partial^2}{\partial\lambda^2}\log p(\boldsymbol{y}\mid\lambda) \right]_{\lambda = \hat\lambda_{\text{MLE}}}\,, \] donde \(p(\boldsymbol{y}\mid\lambda) = \prod_{i=1}^n p(y_i\mid\lambda)\) es la distribución muestral conjunta de \(\boldsymbol{y}\).

La función de densidad conjunta de los datos \(\boldsymbol{y} = (y_1, \ldots, y_n)\) bajo el modelo \(\textsf{Exponencial}(\lambda)\) es: \[ p(\boldsymbol{y} \mid \lambda) = \prod_{i=1}^n \frac{1}{\lambda} e^{-y_i / \lambda} = \lambda^{-n} e^{-\frac{1}{\lambda} \sum_{i=1}^n y_i}. \] Definiendo \(s = \sum_{i=1}^n y_i\), esta verosimilitud conjunta se escribe como \(p(\boldsymbol{y} \mid \lambda) = \lambda^{-n} e^{-s / \lambda}\). Para encontrar el estimador de máxima verosimilitud (MLE) de \(\lambda\), consideramos el logaritmo de la función de verosimilitud: \[ \log p(\boldsymbol{y} \mid \lambda) = -n \log \lambda - \frac{s}{\lambda}. \] Derivamos esta función con respecto a \(\lambda\) para encontrar los puntos críticos: \[ \frac{\partial}{\partial \lambda} \log p(\boldsymbol{y} \mid \lambda) = -\frac{n}{\lambda} + \frac{s}{\lambda^2}. \] Igualando esta derivada a cero: \[ -\frac{n}{\lambda} + \frac{s}{\lambda^2} = 0 \quad \implies \quad s = n\lambda \quad \implies \quad \lambda = \frac{s}{n} = \bar{y}. \] Aquí, \(\bar{y}\) es la media muestral, definida como \(\bar{y} = \frac{1}{n} \sum_{i=1}^n y_i\). Ahora verificamos que \(\lambda = \bar{y}\) corresponde a un máximo. La segunda derivada del logaritmo de la verosimilitud es: \[ \frac{\partial^2}{\partial \lambda^2} \log p(\boldsymbol{y} \mid \lambda) = \frac{n}{\lambda^2} - \frac{2s}{\lambda^3}. \] Sustituyendo \(s = n\lambda\) en esta expresión: \[ \frac{\partial^2}{\partial \lambda^2} \log p(\boldsymbol{y} \mid \lambda) = \frac{n}{\lambda^2} - \frac{2n\lambda}{\lambda^3} = \frac{n}{\lambda^2} - \frac{2n}{\lambda^2} = -\frac{n}{\lambda^2}. \] Dado que \(-\frac{n}{\lambda^2} < 0\) para \(\lambda > 0\), concluimos que \(\lambda = \bar{y}\) es un máximo local. Como la verosimilitud es unimodal, este es también un máximo global.

La información observada de Fisher se define como: \[ \hat{I}(\hat{\lambda}_{\text{MLE}}) = \left[ -\frac{\partial^2}{\partial \lambda^2} \log p(\boldsymbol{y} \mid \lambda) \right]_{\lambda = \hat{\lambda}_{\text{MLE}}}. \] Sustituyendo la segunda derivada calculada y evaluándola en \(\lambda = \hat{\lambda}_{\text{MLE}} = \bar{y}\), obtenemos: \[ \hat{I}(\hat{\lambda}_{\text{MLE}}) = -\left(-\frac{n}{\bar{y}^2}\right) = \frac{n}{\bar{y}^2}. \]

En conclusión, el estimador de máxima verosimilitud de \(\lambda\) es \(\hat{\lambda}_{\text{MLE}} = \bar{y}\), y la información observada de Fisher evaluada en este estimador es \(\hat{I}(\hat{\lambda}_{\text{MLE}}) = \frac{n}{\bar{y}^2}\).

  1. Encuentre la estimación, el coeficiente de variación y un intervalo de de credibilidad/confianza al 95% para \(\lambda\), bajo el paradigma Bayesiano como Frecuentista (usando la Normalidad asintóica del MLE, Bootstrap paramétrico y Bootstrap no paramétrico).

Recuerde que asintóticamente se tiene que \(\hat\lambda_{\text{MLE}}\sim\textsf{N}(\lambda,\hat{I}^{-1}(\hat\lambda_{\text{MLE}}))\).

Estimación de \(\lambda\) bajo el paradigma Bayesiano

En el paradigma Bayesiano, asumimos que \(\lambda \mid \boldsymbol{y} \sim \textsf{GI}(a^*, b^*)\), con parámetros actualizados: \[ a^* = a + n, \quad b^* = b + \sum_{i=1}^n y_i = b + s. \] La media posterior de \(\lambda\) es: \[ \textsf{E}(\lambda \mid \boldsymbol{y}) = \frac{b^*}{a^* - 1}. \] El coeficiente de variación (CV) posterior se calcula como: \[ \textsf{CV}(\lambda \mid \boldsymbol{y}) = \frac{\sqrt{\textsf{Var}(\lambda \mid \boldsymbol{y})}}{\textsf{E}(\lambda \mid \boldsymbol{y})}, \] donde: \[ \textsf{Var}(\lambda \mid \boldsymbol{y}) = \frac{b^{*2}}{(a^* - 1)^2(a^* - 2)} \quad \text{para } a^* > 2. \] El intervalo de credibilidad al 95% para \(\lambda\) es el intervalo central basado en la distribución \(\textsf{GI}(a^*, b^*)\), que puede calcularse numéricamente a partir de la densidad de la distribución Gamma Inversa.

Estimación de \(\lambda\) bajo el paradigma Frecuentista

El estimador de máxima verosimilitud (MLE) de \(\lambda\) es \(\hat{\lambda}_{\text{MLE}} = \bar{y}\). Bajo la normalidad asintótica del MLE, sabemos que: \[ \hat{\lambda}_{\text{MLE}} \sim \textsf{N}(\lambda, \hat{I}^{-1}(\hat{\lambda}_{\text{MLE}})), \] donde la información observada de Fisher es: \[ \hat{I}(\hat{\lambda}_{\text{MLE}}) = \frac{n}{\hat{\lambda}_{\text{MLE}}^2}. \] La varianza asintótica del MLE es \(\textsf{Var}(\hat{\lambda}_{\text{MLE}}) = \hat{I}^{-1}(\hat{\lambda}_{\text{MLE}}) = \frac{\hat{\lambda}_{\text{MLE}}^2}{n}\). Por tanto, un intervalo de confianza asintótico al 95% para \(\lambda\) es: \[ \left(\hat{\lambda}_{\text{MLE}} - 1.96 \cdot \sqrt{\frac{\hat{\lambda}_{\text{MLE}}^2}{n}}, \hat{\lambda}_{\text{MLE}} + 1.96 \cdot \sqrt{\frac{\hat{\lambda}_{\text{MLE}}^2}{n}}\right). \]

El coeficiente de variación (CV) se calcula como: \[ \textsf{CV}(\hat{\lambda}_{\text{MLE}}) = \frac{\sqrt{\textsf{Var}(\hat{\lambda}_{\text{MLE}})}}{\hat{\lambda}_{\text{MLE}}} = \frac{\sqrt{\frac{\hat{\lambda}_{\text{MLE}}^2}{n}}}{\hat{\lambda}_{\text{MLE}}} = \frac{1}{\sqrt{n}}. \]

Bootstrap paramétrico

En el enfoque bootstrap paramétrico, generamos múltiples muestras replicadas \(\boldsymbol{y}^*\) de tamaño \(n\), simuladas a partir de la distribución \(\textsf{Exponencial}(\hat{\lambda}_{\text{MLE}})\). Para cada muestra replicada, calculamos el estimador \(\hat{\lambda}_{\text{MLE}}^*\). Las replicaciones nos permiten calcular un intervalo de confianza percentil al 95%, definido como el intervalo central que contiene el 95% de los valores replicados.

Bootstrap no paramétrico

En el bootstrap no paramétrico, se generan replicaciones de los datos \(\boldsymbol{y}^*\) mediante muestreo con reemplazo de los datos originales \(\boldsymbol{y}\). Para cada replicación, calculamos \(\hat{\lambda}_{\text{MLE}}^*\), obteniendo una distribución empírica de \(\hat{\lambda}_{\text{MLE}}\). Un intervalo de confianza percentil al 95% se calcula directamente como el intervalo central que contiene el 95% de los valores replicados.

# Cargar librerías necesarias
library(invgamma)  # Para la densidad de la Gamma Inversa
library(knitr)     # Para generar tablas

# Datos
y <- c(495, 541, 1461, 1555, 1603, 2201, 2750, 3468, 3516, 4319, 6622, 7728, 13159, 21194)
n <- length(y)
s <- sum(y)
y_bar <- mean(y)

# Parámetros previos
a_prior <- 8.25  # Asumidos
b_prior <- 32625  # Asumidos

# Bayesiano: Posterior
a_post <- a_prior + n
b_post <- b_prior + s
lambda_mean_post <- b_post / (a_post - 1)
lambda_sd_post <- sqrt(b_post^2 / ((a_post - 1)^2 * (a_post - 2)))
cv_post <- lambda_sd_post / lambda_mean_post

# Intervalo de credibilidad al 95% (percentiles de la Gamma Inversa)
cred_interval <- qinvgamma(c(0.025, 0.975), shape = a_post, rate = b_post)

# Frecuentista: Normalidad asintótica del MLE
lambda_mle <- y_bar
lambda_var <- lambda_mle^2 / n
lambda_sd <- sqrt(lambda_var)
cv_mle <- lambda_sd / lambda_mle
conf_interval <- c(lambda_mle - 1.96 * lambda_sd, lambda_mle + 1.96 * lambda_sd)

# Bootstrap paramétrico
set.seed(123)
B <- 1000
parametric_bootstrap <- replicate(B, mean(rexp(n, rate = 1 / lambda_mle)))
boot_param_interval <- quantile(parametric_bootstrap, c(0.025, 0.975))
cv_boot_param <- sd(parametric_bootstrap) / mean(parametric_bootstrap)

# Bootstrap no paramétrico
nonparam_bootstrap <- replicate(B, mean(sample(y, size = n, replace = TRUE)))
boot_nonparam_interval <- quantile(nonparam_bootstrap, c(0.025, 0.975))
cv_boot_nonparam <- sd(nonparam_bootstrap) / mean(nonparam_bootstrap)

# Crear una tabla con los resultados
resultados <- data.frame(
  Método = c("Posterior Bayesiano", 
             "MLE Frecuentista", 
             "Bootstrap Paramétrico", 
             "Bootstrap No Paramétrico"),
  Estimación = c(lambda_mean_post, lambda_mle, lambda_mle, lambda_mle),
  `CV` = c(cv_post, cv_mle, cv_boot_param, cv_boot_nonparam),
  `Límite Inferior 95%` = c(cred_interval[1], conf_interval[1], boot_param_interval[1], boot_nonparam_interval[1]),
  `Límite Superior 95%` = c(cred_interval[2], conf_interval[2], boot_param_interval[2], boot_nonparam_interval[2])
)

# Mostrar la tabla con knitr
kable(resultados, digits = 2, col.names = c("Método", "Estimación", "CV", "Límite Inferior 95%", "Límite Superior 95%"),
      caption = "Estimaciones, coeficiente de variación (CV) e intervalos de credibilidad/confianza al 95% para λ.")
Estimaciones, coeficiente de variación (CV) e intervalos de credibilidad/confianza al 95% para λ.
Método Estimación CV Límite Inferior 95% Límite Superior 95%
Posterior Bayesiano 4858.21 0.22 3186.03 7381.96
MLE Frecuentista 5043.71 0.27 2401.66 7685.77
Bootstrap Paramétrico 5043.71 0.26 2844.19 7995.15
Bootstrap No Paramétrico 5043.71 0.28 2587.97 8156.07