1. Sean \(x\), \(y\) y \(z\) variables aleatorias con función de densidad conjunta (discreta o continua) dada por \(p(x,y,z) \propto p(x,z)p(y,z)p(z)\). Muestre que:

    1. \(p(x\mid y,z)\propto p(x,z)\), i.e., \(p(x\mid y,z)\) es función de \(x\) y \(z\).
    2. \(p(y\mid x,z)\propto p(y,z)\), i.e., \(p(y\mid x,z)\) es función de \(y\) y \(z\).
    3. \(x\) y \(y\) son condicionalmente independientes dado \(z\).

  2. Sean \(A\), \(B\), y \(C\) proposiciones de falso-verdadero. Suponga que \(A\) y \(B\) son condicionalmente independientes, dado \(C\). Muestre que:

    1. \(A\) y \(B^C\) son condicionalmente independientes, dado \(C\).
    2. \(A^C\) y \(B^C\) son condicionalmente independientes, dado \(C\).

  3. Muestre que si \(y\mid\theta\) tiene distribución exponencial con parámetro \(\theta\), entonces la distribución \(\textsf{Gamma}\) sirve como distribución previa conjugada para hacer inferencias sobre \(\theta\), dada una muestra aleatoria de valores de \(y\).

  4. Suponga que su estado de información previo para \(\theta\), la proporción de individuos que apoyan la pena de muerte en California, es \(\textsf{Beta}\) con media \(\textsf{E}(\theta) = 0.6\) y desviación estándar \(\textsf{DE}(\theta) = 0.3\).

    1. Determine los hiperparámetros de la distribución previa y dibuje la función de densidad correspondiente.
    2. Se toma una muestra aleatoria de 1,000 californianos y el 65% apoya la pena de muerte. Calcule tanto la media como la desviación estándar posterior para \(\theta\). Dibuje la función de densidad posterior correspondiente.
    3. Examine la sensibilidad de la distribución posterior a diferentes valores de la media y de la desviación estándar a priori, incluyendo una distribución previa no informativa.

  5. Un ingeniero está inspeccionando un gran envío de piezas con fines de control de calidad y decide probar diez elementos seleccionados al azar. Históricamente, la proporción de artículos defectuosos \(\theta\) ha sido de alrededor del 1% y muy rara vez ha estado por encima del 2%.

    1. Determine una distribución previa conjugada para \(\theta\) de acuerdo con la información histórica, y además, usando esta distribución previa, encuentre la distribución posterior de \(\theta\) dada una muestra aleatoria de tamaño diez.
    2. Suponga que el ingeniero no encuentra componentes defectuosos en su proceso de observación. ¿Cuál es la distribución posterior de \(\theta\)? ¿Cuál es la media posterior de \(\theta\)?
    3. Calcule el estimador de máxima verosimilitud para \(\theta\). Como estimador puntual, ¿en este caso es preferible el estimador máximo verosímil o la media posterior? ¿Por qué?

  6. Considere el modelo Beta-Binomial \[ y\mid\theta \sim \textsf{Binomial}(n,\theta) \qquad\text{y}\qquad \theta \sim \textsf{Beta}(a,b) \] con \(y\in\mathcal{Y}=\{0,\ldots,n\}\) y \(\theta\in\Theta=[0,1]\).

    1. Muestre que la distribución marginal de \(y\) es \[ p(y) = \frac{\Gamma(n+1)}{\Gamma(y+1)\,\Gamma(n-y+1)}\,\frac{\Gamma(a+b)}{\Gamma(a+b+n)}\,\frac{\Gamma(a+y)\,\Gamma(b+n-y)}{\Gamma(a)\,\Gamma(b)}\,. \]
    2. Muestre que la media marginal y la varianza marginal de \(y\) son respectivamente \[ \textsf{E}(y) = \frac{na}{a+b} \qquad\text{y}\qquad \textsf{Var}(y) = \frac{nab(a+b+n)}{(a+b)^2(a+b+1)}\,. \] Sugerencia: \[ \textsf{E}(X) = \textsf{E}(\textsf{E}(X\mid Y)) \qquad\text{y}\qquad \textsf{Var}(X) = \textsf{E}(\textsf{Var}(X\mid Y)) + \textsf{Var}(\textsf{E}(X\mid Y))\,. \]
    3. Muestre que el promedio posterior de \(\theta\) es un promedio ponderado entre la media previa de \(\theta\) y el número de éxitos promedio, es decir, \[ \textsf{E}(\theta\mid y) = \omega\,\textsf{E}(\theta) + (1-\omega)\,\bar{y} \] donde \(\omega = \frac{b}{b+n}\), \(1 - \omega = \frac{n}{b+n}\) y \(\bar{y} = y/n\).

  7. Muestre que las distribuciones Bernoulli, Binomial, Multinomial, Poisson, Exponencial, Beta, Gamma, y Normal hacen parte de la familia exponencial.

  8. Jeffreys (H. Jeffreys (1961). Theory of Probability. Oxford Univ. Press.) sugirió una regla para generar una distribución previa de un parámetro \(\theta\) asociado con la distribución muestral \(p(y\mid\theta)\). La distribución previa de Jeffreys es de la forma \(p_J(\theta)\propto\sqrt{I(\theta)}\), donde \[ I(\theta) = -\textsf{E}_{y\mid\theta}\left( \frac{\text{d}^2}{\text{d}\theta^2}\log p(y\mid\theta) \right) \] es la información esperada de Fisher.

    1. Sea \(y\mid\theta\sim\textsf{Bin}(n,\theta)\). Muestre que la distribución previa de Jeffreys para esta distribución muestral es \[ p_J(\theta)\propto \theta^{-\frac12} (1-\theta)^{-\frac12}\,. \]
    2. Reparametrice la distribución Binomial con \(\psi = \textsf{logit}(\theta)\), de forma que \[p(y\mid\psi) \propto e^{\psi y}(1+e^\psi)^{-n}\,.\] Obtenga la distribución previa de Jeffreys para esta distribución muestral.
    3. Tome la distribución previa de la parte a. y aplique la fórmula del cambio de variables para obtener la densidad previa de \(\psi\). Esta densidad debe coincidir con la obtenida en el inciso b.. Esta propiedad de invarianza bajo reparametrización es la característica fundamental de la previa de Jeffreys.

  9. El estimador óptimo del parámetro \(\theta\in \Theta\subset\mathbb{R}\) de acuerdo con la regla de Bayes se define como el estimador \(\hat\theta=\hat\theta(\boldsymbol{y})\) que minimiza la perdida esperada posterior dada por \[ \textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta)) = \int_{\Theta} L(\theta,\hat\theta)\, p(\theta\mid\boldsymbol{y})\, \textsf{d}\theta\,, \] donde \(L(\theta,\hat\theta)\) es una función de perdida (costo que conlleva estimar \(\theta\) por medio de \(\hat\theta\)) y \(\boldsymbol{y}=(y_1,\ldots, y_n)\) es un conjunto de datos observado.

    1. Muestre que si \(L(\theta,\hat\theta) = (\theta - \hat\theta)^2\), entonces el estimador óptimo de acuerdo con la regla de Bayes es la media posterior \(\hat\theta=\textsf{E}(\theta\mid\boldsymbol{y})\).
    2. Muestre que si \(L(\theta,\hat\theta) = |\theta - \hat\theta|\), entonces el estimador óptimo de acuerdo con la regla de Bayes es la mediana posterior \(\hat\theta=(\theta\mid\boldsymbol{y})_{0.5}\).
    3. El riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) se define como \[ R_{\textsf{F}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}\mid\theta}(L(\theta,\hat\theta)) = \int_{\mathcal{Y}} L(\theta,\hat\theta)\, p(\boldsymbol{y}\mid\theta)\,\textsf{d}\boldsymbol{y}\,, \] i.e., el valor medio de la función de perdida \(L(\theta,\hat\theta)\) a través de todos los valores de \(\boldsymbol{y} \in \mathcal{Y}\). De otra parte, el riesgo Bayesiano \(R_{\textsf{B}}(\theta,\hat\theta)\) se define como \[ 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\,, \] i.e., el valor medio del riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) a priori a través de todos los valores de \(\theta\in\Theta\). Muestre que \[ R_{\textsf{B}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}}(\textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta)))\,, \] donde \(\textsf{E}_{\boldsymbol{y}}(\cdot)\) denota el valor esperado respecto a la distribución marginal \(p(\boldsymbol{y})\).

  10. Suponga que \(y_i\mid\boldsymbol{\theta}\stackrel{\text{iid}}{\sim} p(y_i\mid\boldsymbol{\theta})\), para \(i=1,\ldots,n\), con \(\boldsymbol{\theta}\sim p(\boldsymbol{\theta})\), donde \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_k)\) y \(\boldsymbol{y} = (y_1\ldots,y_n)\).

    Sea \(\boldsymbol{\theta}_\text{MAP}\) el máximo a posteriori (MAP, por sus siglas en inglés), el cual se define como el valor de \(\boldsymbol{\theta}\) que maximiza la distribución posterior posterior, a saber, \[ \boldsymbol{\theta}_\text{MAP} = \textsf{arg max}_{\boldsymbol{\theta}}\,\log p(\boldsymbol{\theta}\mid\boldsymbol{y}) = \textsf{arg max}_{\boldsymbol{\theta}}\,\left(\log p(\boldsymbol{y}\mid\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})\right)\,, \] y además, sea \(\mathbf\Sigma_\text{MAP} = -\mathbf{H}^{-1}\), donde \(\mathbf{H}\) es la matriz Hessiana cuya \((i,j)\)-ésima entrada está dada por \[ H_{i,j} = \frac{\partial^2}{\partial\theta_i\,\partial\theta_j}\,\log p (\boldsymbol{\theta}\mid\boldsymbol{y}) \Bigg|_{\boldsymbol{\theta}=\boldsymbol{\theta}_\text{MAP}} = \frac{\partial^2}{\partial\theta_i\,\partial\theta_j}\,\left(\log p(\boldsymbol{y}\mid\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})\right)\Bigg|_{\boldsymbol{\theta}=\boldsymbol{\theta}_\text{MAP}}\,. \] Se observa que la constante de normalización \(p(\boldsymbol{y})\) no depende de \(\boldsymbol{\theta}\), y por lo tanto, no interfiere en la maximización de \(\log p(\boldsymbol{\theta}\mid\boldsymbol{y})\). El Teorema del Límite Central Bayesiano (BCLT, por sus siglas en inglés) indica que \(\boldsymbol{\theta}\mid \boldsymbol{y} \approx \textsf{N}_k(\boldsymbol{\theta}_\text{MAP},\mathbf\Sigma_\text{MAP})\), cuando \(n\rightarrow\infty\).

    Usando los datos del caso de victimas de violencia sexual, aproxime la distribución posterior de \(\theta\) por medio del BCLT. Dibuje la distribución exacta y la aproximada en un mismo gráfico.

  11. (Teorema de De Finetti para variables binarias) Sea \((x_1, x_2, \dots)\) una secuencia infinita de variables aleatorias binarias intercambiable (i.e., \(x_1,\ldots,x_n\) es intercambiable para todo \(n\in\mathbb{N}\)). Demuestre que existe una variable aleatoria \(\theta\) con soporte en \((0, 1)\) tal que \[ p(x_1,\ldots,x_n) = \int_0^1 \prod_{i=1}^n \theta^{x_i} (1 - \theta)^{1 - x_i} \, p(\theta) \, \textsf{d}(\theta), \] donde \(p(\theta)\) es una distribución de probabilidad sobre \((0, 1)\).

Solución

  1. Sean \(x\), \(y\) y \(z\) variables aleatorias con función de densidad conjunta (discreta o continua) dada por \(p(x,y,z) \propto p(x,z)p(y,z)p(z)\). Muestre que:
  1. \(p(x\mid y,z)\propto p(x,z)\), i.e., \(p(x\mid y,z)\) es función de \(x\) y \(z\).

Se tiene que \[ p(x\mid y,z) = \frac{p(x,y,z)}{p(y,z)} = \frac{p(x,z)p(y,z)p(z)}{p(y,z)} \propto p(x,z)\,, \] y por lo tanto \(p(x\mid y,z)\) es función de \(x\) y \(z\).

  1. \(p(y\mid x,z)\propto p(y,z)\), i.e., \(p(y\mid x,z)\) es función de \(y\) y \(z\).

Se tiene que \[ p(y\mid x,z) = \frac{p(x,y,z)}{p(x,z)} = \frac{p(x,z)p(y,z)p(z)}{p(x,z)} \propto p(y,z)\,, \] y por lo tanto p(yx,z)$ es función de \(y\) y \(z\).

  1. \(x\) y \(y\) son condicionalmente independientes dado \(z\).

Se tiene que \[ p(x,y\mid z) = \frac{p(x,y,z)}{p(z)} = \frac{p(x,z)p(y,z)p(z)}{p(z)} = \frac{p(x,z)p(z)}{p(z)}\,\frac{p(y,z)p(z)}{p(z)} \propto p(x\mid z)\,p(y\mid z)\,, \] y por lo tanto \(x\) y \(y\) son condicionalmente independientes dado \(z\).

  1. Sean \(A\), \(B\), y \(C\) proposiciones de falso-verdadero. Suponga que \(A\) y \(B\) son condicionalmente independientes, dado \(C\). Muestre que:
  1. \(A\) y \(B^C\) son condicionalmente independientes, dado \(C\).

Si \(A\) y \(B\) son condicionalmente independientes, dado \(C\), entonces \(\textsf{Pr}(B\mid A,C) = \textsf{Pr}(B\mid C)\), y por lo tanto \(1-\textsf{Pr}(B\mid A,C) = 1-\textsf{Pr}(B\mid C)\), y así, \(\textsf{Pr}(B^C\mid A,C) = \textsf{Pr}(B^C\mid C)\). Luego, \(A\) y \(B^C\) son condicionalmente independientes, dado \(C\).

  1. \(A^C\) y \(B^C\) son condicionalmente independientes, dado \(C\).

\(A^C\) y \(B^C\) son condicionalmente independientes, dado \(C\), dado que \[ \begin{align*} \textsf{Pr}(A^C\cap B^C\mid C) &= \textsf{Pr}((A\cup B)^C\mid C)\\ &= 1-(\textsf{Pr}(A\mid C)+\textsf{Pr}(B\mid C)-\textsf{Pr}(A\cap B\mid C))\\ &= 1-\textsf{Pr}(A\mid C)-\textsf{Pr}(B\mid C))+\textsf{Pr}(A\mid C)\textsf{Pr}(B\mid C)\,. \end{align*} \] Luego, \(\textsf{Pr}(A^C\cap B^C\mid C) = 1-\textsf{Pr}(A\mid C)-\textsf{Pr}(B\mid C)(1 - \textsf{Pr}(A\mid C))= (1 - \textsf{Pr}(A\mid C))(1 - \textsf{Pr}(B\mid C))\), de donde se sigue que \(A^C\) y \(B^C\) son condicionalmente independientes, dado \(C\).

  1. Muestre que si \(y\mid\theta\) tiene distribución exponencial con parámetro \(\theta\), entonces la distribución \(\textsf{Gamma}\) sirve como distribución previa conjugada para hacer inferencias sobre \(\theta\), dada una muestra aleatoria de valores de \(y\).

Si \(y_i\stackrel{\text{iid}}{\sim}\theta\sim \textsf{Exp}(\theta)\), para \(i=1,\ldots,n\), y \(\theta\sim\textsf{Gamma}(a,b)\), entonces la distribución posterior de \(\theta\) es: \[ \begin{align*} p(\theta\mid\boldsymbol{y}) &\propto p(\boldsymbol{y}\mid\theta)\,p(\theta)\\ &= \prod_{i=1}^n \textsf{Exp}(y_i\mid\theta) \,\textsf{Gamma}(\theta)\\ &= \prod_{i=1}^n \theta e^{-\theta\,y_i} \,\frac{b^a}{\Gamma(a)}\,\theta^{a-1}\,e^{-b\theta}\\ &\propto \theta^n e^{-\theta s}\theta^{a-1}\,e^{-b\theta}\\ &= \theta^{(a+n)-1} e^{-(b+s)\theta}\,, \end{align*} \] donde \(\boldsymbol{y}=(y_1,\ldots,y_n)\) y \(s=\sum_{i=1}^n y_i\). Por lo tanto \(\theta\mid\boldsymbol{y} \sim \textsf{Gamma}(a + n, b +s)\) dado que \(\theta^{(a+n)-1} e^{-(b+s)\theta}\) corresponde al núcleo de una distribución Gamma con parámetros \(\alpha=a+n\) y \(\beta = b + s\). Así, la distribución Gamma sirve como distribución previa conjugada para hacer inferencias sobre \(\theta\) porque la distribución posterior de \(\theta\) también pertenece a la familia de distribuciones Gamma.

  1. Suponga que su estado de información previo para \(\theta\), la proporción de individuos que apoyan la pena de muerte en California, es \(\textsf{Beta}\) con media \(\textsf{E}(\theta) = 0.6\) y desviación estándar \(\textsf{DE}(\theta) = 0.3\).
  1. Determine los hiperparámetros de Su distribución previa y dibuje la función de densidad correspondiente.

La información externa indica que \(\textsf{E}(\theta) = \frac{a}{a+b} = 0.6\) y \(\textsf{Var}(\theta) = \frac{ab}{(a+b)^2(a+b+1)} = 0.3^2\), el cual corresponde a un sistema de ecuaciones para \(a\) y \(b\). Así, se tiene que \(a+b = 5a/3\) de donde \(b = 2a/3\) y \(\frac{2a/3}{(25a/9)(5a/3+1)} = 0.3^2\). Despejando se tiene que \(a = 1\) y por lo tanto \(b = 2/3\).

# hiperparametro a
a <- (3/5)*((2/3)/((25/9)*0.3^2) - 1)
a
## [1] 1
# hiperparametro b
b <- 2*a/3
b
## [1] 0.6666667
# media
a/(a+b)
## [1] 0.6
# varianza
a*b/((a+b)^2*(a+b+1))
## [1] 0.09
# grafico de la previa
curve(expr = dbeta(x, shape1 = a, shape2 = b), from = 0.001, to = 0.999, lwd = 2, xlab = expression(theta), ylab = expression(paste("p","(",theta,")",sep="")), main = "Distr. previa")

  1. Se toma una muestra aleatoria de 1,000 californianos y el 65% apoya la pena de muerte. Calcule tanto la media como la desviación estándar posterior para \(\theta\). Dibuje la función de densidad posterior correspondiente.

La distribución posterior de \(\theta\) es \(\theta\mid s \sim\textsf{Beta}(a + s, b + n -s)\). Dado que \(a = 1\), \(b = 2/3\), \(n=1000\) y \(s=\sum_{i=1}^n y_i=650\), entonces \(\textsf{E}(\theta\mid s) = 0.6499\) y \(\textsf{DE}(\theta\mid s) = 0.0151\).

# distribucion posterior
a <- 1
b <- 2/3
n <- 1000
s <- 0.65*n
ap <- a + s
bp <- b + n - s
# media posterior
round(ap/(ap + bp), 4)
## [1] 0.6499
# DE posterior
round(sqrt(ap*bp/((ap+bp)^2*(ap+bp+1))), 4)
## [1] 0.0151
# grafico de la posterior
par(mfrow = c(1,2))
curve(expr = dbeta(x, shape1 = ap, shape2 = bp), from = 0.001, to = 0.999, lwd = 2, xlab = expression(theta), ylab = expression(paste("p","(",theta," | ",y,")",sep="")), main = "Distr. posterior")
curve(expr = dbeta(x, shape1 = ap, shape2 = bp), from = 0.6, to = 0.7, lwd = 2, xlab = expression(theta), ylab = expression(paste("p","(",theta," | ",y,")",sep="")), main = "Distr. posterior")

  1. Examine la sensibilidad de la distribución posterior a diferentes valores de la media y de la desviación estándar a priori, incluyendo una distribución previa no informativa.

Se observa que utilizando la previa uniforme o la previa de Jeffreys los resultados a posteriori son casi idénticos. Esto ocurre porque el tamaño de muestra es considerablemente grande, lo que opaca la información proveniente de la distribución previa.

# previa original
a  <- 1
b  <- 2/3
# previa uniforme
a1 <- 1
b1 <- 1
# previa Jeffreys
a2 <- 1/2
b2 <- 1/2
# datos
n  <- 1000
s  <- 0.65*n
# grafico de la posterior
curve(expr = dbeta(x, shape1 = a + s, shape2 = b + n - s), from = 0.6, to = 0.7, lwd = 2, lty = 1, xlab = expression(theta), ylab = expression(paste("p","(",theta," | ",y,")",sep="")), main = "Distr. posterior")
curve(expr = dbeta(x, shape1 = a1 + s, shape2 = b1 + n - s), from = 0.6, to = 0.7, lwd = 2, col = 2, lty = 2, add = T)
curve(expr = dbeta(x, shape1 = a2 + s, shape2 = b2 + n - s), from = 0.6, to = 0.7, lwd = 2, col = 3, lty = 3, add = T)
legend("topright", legend = c("a = 1, b = 2/3","a = 1, b = 1","a = 1/2, b = 1/2"), col = c(1,2,3), lwd = 2, lty = c(1,2,3), bty = "n")

  1. Un ingeniero está inspeccionando un gran envío de piezas con fines de control de calidad y decide probar diez elementos seleccionados al azar. Históricamente, la proporción de artículos defectuosos \(\theta\) ha sido de alrededor del 1% y muy rara vez ha estado por encima del 2%.
  1. Determine una distribución previa conjugada para \(\theta\) de acuerdo con la información histórica, y además, usando esta distribución previa, encuentre la distribución posterior de \(\theta\) dada una muestra aleatoria de tamaño diez.

A priori se sabe que \(\textsf{E}(\theta)=0.01\) y \(\textsf{Pr}(\theta < 0.02) \approx 0.97\). Este corresponde a un sistema de ecuaciones no lineales de dos incógnitas con infinitas soluciones. Por ejemplo, tomando \(a=5\) y \(b=495\) se obtienen las restricciones requeridas a priori (por su puesto otras soluciones también son posibles).

Usando la distribución previa \(\theta\sim\textsf{Beta}(5,495)\) se tiene que la distribución posterior de \(\theta\) para una muestra de tamaño \(n=10\) es \(\theta\mid s \sim\textsf{Beta}(5 + s, 505 - s)\).

# hiperparámetros
a <- 5
b <- 0.99*a/0.01
b
## [1] 495
# media
a/(a+b)
## [1] 0.01
# acumulada a 0.02 
pbeta(q = 0.02, shape1 = a, shape2 = b, lower.tail = T)
## [1] 0.9715067
  1. Suponga que el ingeniero no encuentra componentes defectuosos en su proceso de observación. ¿Cuál es la distribución posterior de \(\theta\)? ¿Cuál es la media posterior de \(\theta\)?

Si \(s=0\) entonces la distribución posterior de \(\theta\) es \(\theta\mid s \sim\textsf{Beta}(5, 505)\) y por lo tanto la media posterior de \(\theta\) es \(\textsf{E}(\theta\mid s) = 0.0098\).

# hiperparámetros
a <- 5
b <- 0.99*a/0.01
# posterior
n <- 10
s <- 0
ap <- a + s
bp <- b + n - s
# media posterior
ap/(ap+bp)
## [1] 0.009803922
  1. Calcule el estimador de máxima verosimilitud para \(\theta\). Como estimador puntual , ¿en este caso es preferible el estimador máximo verosímil o la media posterior? ¿Por qué?

En este caso la estimación máximo verosímil de \(\theta\) es \(\bar{y}= 0\). La media posterior es un mejor estimador porque, 1. tiene en consideración la información histórica acerca de las componentes defectuosas, y 2. no tiene sentido estimar la proporción poblacional de componentes defectuosos como cero dado que aunque las componentes defectuosas son pocas claramente a nivel poblacional no son exactamente iguales a cero.

  1. Considere el modelo Beta-Binomial \[ y\mid\theta \sim \textsf{Binomial}(n,\theta) \qquad\text{y}\qquad \theta \sim \textsf{Beta}(a,b) \] con \(y\in\mathcal{Y}=\{0,\ldots,n\}\) y \(\theta\in\Theta=[0,1]\).
  1. Muestre que la distribución marginal de \(y\) es \[ p(y) = \frac{\Gamma(n+1)}{\Gamma(y+1)\,\Gamma(n-y+1)}\,\frac{\Gamma(a+b)}{\Gamma(a+b+n)}\,\frac{\Gamma(a+y)\,\Gamma(b+n-y)}{\Gamma(a)\,\Gamma(b)}\,. \]

La distribución marginal de \(y\) se obtiene integrando la probabilidad conjunta \(p(y, \theta)\) sobre el espacio del parámetro \(\theta\). En el modelo dado, \(y \mid \theta \sim \textsf{Binomial}(n, \theta)\), cuya función de probabilidad es \(p(y \mid \theta) = \binom{n}{y} \theta^y (1-\theta)^{n-y}\), y \(\theta \sim \textsf{Beta}(a, b)\), cuya densidad es \(p(\theta) = \frac{\theta^{a-1} (1-\theta)^{b-1}}{B(a, b)}\), donde \(B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\). Por tanto, la marginal de \(y\) se expresa como \[ p(y) = \int_0^1 p(y \mid \theta) p(\theta) \, \text{d}\theta. \]

Sustituyendo las expresiones de \(p(y \mid \theta)\) y \(p(\theta)\), tenemos: \[ p(y) = \int_0^1 \binom{n}{y} \theta^y (1-\theta)^{n-y} \cdot \frac{\theta^{a-1} (1-\theta)^{b-1}}{B(a, b)} \, \text{d}\theta. \] Agrupando términos comunes dentro de la integral, la expresión se convierte en: \[ p(y) = \binom{n}{y} \frac{1}{B(a, b)} \int_0^1 \theta^{y+a-1} (1-\theta)^{n-y+b-1} \, \text{d}\theta. \]

La integral tiene la forma de una función beta, definida como: \[ B(x, y) = \int_0^1 t^{x-1} (1-t)^{y-1} \, \text{d}t = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)}. \] Comparando con la integral presente, identificamos los parámetros \(x = y+a\) y \(y = n-y+b\). Por tanto, la integral se evalúa como: \[ \int_0^1 \theta^{y+a-1} (1-\theta)^{n-y+b-1} \, \text{d}\theta = B(y+a, n-y+b) = \frac{\Gamma(y+a)\Gamma(n-y+b)}{\Gamma(a+b+n)}. \]

Sustituyendo este resultado en la expresión de \(p(y)\), obtenemos: \[ p(y) = \binom{n}{y} \frac{1}{B(a, b)} \cdot \frac{\Gamma(y+a)\Gamma(n-y+b)}{\Gamma(a+b+n)}. \] Usando que \(B(a, b) = \frac{\Gamma(a)\Gamma(b)}{\Gamma(a+b)}\), el término \(1/B(a, b)\) se convierte en \(\frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)}\). Así, la expresión se transforma en: \[ p(y) = \binom{n}{y} \cdot \frac{\Gamma(a+b)}{\Gamma(a)\Gamma(b)} \cdot \frac{\Gamma(y+a)\Gamma(n-y+b)}{\Gamma(a+b+n)}. \]

Finalmente, sustituimos \(\binom{n}{y} = \frac{\Gamma(n+1)}{\Gamma(y+1)\Gamma(n-y+1)}\) en la ecuación, obteniendo: \[ p(y) = \frac{\Gamma(n+1)}{\Gamma(y+1)\Gamma(n-y+1)} \cdot \frac{\Gamma(a+b)}{\Gamma(a+b+n)} \cdot \frac{\Gamma(y+a)\Gamma(n-y+b)}{\Gamma(a)\Gamma(b)}. \] Esto muestra que la distribución marginal de \(y\) es: \[ p(y) = \frac{\Gamma(n+1)}{\Gamma(y+1)\Gamma(n-y+1)} \cdot \frac{\Gamma(a+b)}{\Gamma(a+b+n)} \cdot \frac{\Gamma(a+y)\Gamma(b+n-y)}{\Gamma(a)\Gamma(b)}. \]

  1. Muestre que la media marginal y la varianza marginal de \(y\) son respectivamente \[ \textsf{E}(y) = \frac{na}{a+b} \qquad\text{y}\qquad \textsf{Var}(y) = \frac{nab(a+b+n)}{(a+b)^2(a+b+1)}\,. \] Sugerencia: \[ \textsf{E}(X) = \textsf{E}(\textsf{E}(X\mid Y)) \qquad\text{y}\qquad \textsf{Var}(X) = \textsf{E}(\textsf{Var}(X\mid Y)) + \textsf{Var}(\textsf{E}(X\mid Y))\,. \]

Para determinar la media y la varianza marginal de \(y\), utilizamos las propiedades de la linealidad de la esperanza y la descomposición de la varianza. En el modelo dado, \(y \mid \theta \sim \textsf{Binomial}(n, \theta)\) y \(\theta \sim \textsf{Beta}(a, b)\). La esperanza marginal de \(y\) se calcula utilizando la propiedad \(\textsf{E}(y) = \textsf{E}[\textsf{E}(y \mid \theta)]\). Dado que \(y \mid \theta\) tiene esperanza condicional \(\textsf{E}(y \mid \theta) = n\theta\), podemos escribir \(\textsf{E}(y) = \textsf{E}[n\theta] = n \textsf{E}(\theta)\). La esperanza de \(\theta\), cuando \(\theta \sim \textsf{Beta}(a, b)\), es \(\textsf{E}(\theta) = \frac{a}{a+b}\). Por lo tanto, la media marginal de \(y\) es: \[ \textsf{E}(y) = n \cdot \frac{a}{a+b} = \frac{na}{a+b}. \]

La varianza marginal de \(y\) se calcula utilizando la propiedad \(\textsf{Var}(y) = \textsf{E}[\textsf{Var}(y \mid \theta)] + \textsf{Var}[\textsf{E}(y \mid \theta)]\). La varianza condicional de \(y \mid \theta\) es \(\textsf{Var}(y \mid \theta) = n\theta(1-\theta)\). Tomando la esperanza de esta varianza condicional, tenemos \(\textsf{E}[\textsf{Var}(y \mid \theta)] = n \cdot \textsf{E}[\theta(1-\theta)]\). Para \(\theta \sim \textsf{Beta}(a, b)\), sabemos que \(\textsf{E}[\theta(1-\theta)] = \textsf{Var}(\theta) = \frac{ab}{(a+b)^2(a+b+1)}\). Sustituyendo, obtenemos: \[ \textsf{E}[\textsf{Var}(y \mid \theta)] = n \cdot \frac{ab}{(a+b)^2(a+b+1)}. \]

Por otro lado, \(\textsf{E}(y \mid \theta) = n\theta\), lo que implica que \(\textsf{Var}[\textsf{E}(y \mid \theta)] = \textsf{Var}(n\theta) = n^2 \textsf{Var}(\theta)\). Usando que \(\textsf{Var}(\theta) = \frac{ab}{(a+b)^2(a+b+1)}\), tenemos: \[ \textsf{Var}[\textsf{E}(y \mid \theta)] = n^2 \cdot \frac{ab}{(a+b)^2(a+b+1)}. \]

Sumando ambos términos, obtenemos la varianza marginal de \(y\): \[ \textsf{Var}(y) = \textsf{E}[\textsf{Var}(y \mid \theta)] + \textsf{Var}[\textsf{E}(y \mid \theta)] = n \cdot \frac{ab}{(a+b)^2(a+b+1)} + n^2 \cdot \frac{ab}{(a+b)^2(a+b+1)}. \] Factorizando el término común, la expresión resulta en: \[ \textsf{Var}(y) = \frac{ab}{(a+b)^2(a+b+1)} \cdot (n + n^2). \] Simplificando, obtenemos: \[ \textsf{Var}(y) = \frac{nab(a+b+n)}{(a+b)^2(a+b+1)}. \]

En conclusión, la media marginal de \(y\) es \(\textsf{E}(y) = \frac{na}{a+b}\), y la varianza marginal de \(y\) es \(\textsf{Var}(y) = \frac{nab(a+b+n)}{(a+b)^2(a+b+1)}\). Estos resultados reflejan la influencia de los parámetros \(a\) y \(b\) de la distribución beta y el tamaño \(n\) de la binomial en las características de \(y\).

  1. Muestre que el promedio posterior de \(\theta\) es un promedio ponderado entre la media previa de \(\theta\) y el número de éxitos promedio, es decir, \[ \textsf{E}(\theta\mid y) = \omega\,\textsf{E}(\theta) + (1-\omega)\,\bar{y} \] donde \(\omega = \frac{b}{b+n}\), \(1 - \omega = \frac{n}{b+n}\) y \(\bar{y} = y/n\).

Para demostrar que el promedio posterior de \(\theta\) es un promedio ponderado entre la media previa de \(\theta\) y el número de éxitos promedio \(\bar{y}\), comenzamos analizando el modelo jerárquico y las propiedades de la distribución posterior. Recordemos que el modelo se define como \(y \mid \theta \sim \textsf{Binomial}(n, \theta)\) y \(\theta \sim \textsf{Beta}(a, b)\).

La distribución posterior de \(\theta\) está dada por: \[ p(\theta \mid y) \propto p(y \mid \theta)p(\theta). \] Sustituyendo las expresiones de \(p(y \mid \theta)\) y \(p(\theta)\), tenemos: \[ p(\theta \mid y) \propto \binom{n}{y} \theta^y (1-\theta)^{n-y} \cdot \frac{\theta^{a-1}(1-\theta)^{b-1}}{B(a, b)}. \] Agrupando los términos que dependen de \(\theta\), la posterior de \(\theta\) es proporcional a: \[ p(\theta \mid y) \propto \theta^{y + a - 1} (1-\theta)^{n - y + b - 1}. \] Reconocemos que esto corresponde a la forma de una distribución Beta con parámetros \(a^* = a + y\) y \(b^* = b + n - y\), es decir: \[ \theta \mid y \sim \textsf{Beta}(a + y, b + n - y). \]

La media posterior de \(\theta\) se calcula como la media de una distribución Beta, que está dada por: \[ \textsf{E}(\theta \mid y) = \frac{a^*}{a^* + b^*}. \] Sustituyendo \(a^* = a + y\) y \(b^* = b + n - y\), obtenemos: \[ \textsf{E}(\theta \mid y) = \frac{a + y}{a + y + b + n - y}. \] Simplificando el denominador: \[ \textsf{E}(\theta \mid y) = \frac{a + y}{a + b + n}. \]

Para reescribir esta expresión como un promedio ponderado, recordamos que la media previa de \(\theta\) es: \[ \textsf{E}(\theta) = \frac{a}{a + b}, \] y definimos el promedio de éxitos como \(\bar{y} = y/n\). Reescribimos la media posterior en términos de \(\textsf{E}(\theta)\) y \(\bar{y}\). Multiplicando y dividiendo \(\textsf{E}(\theta \mid y)\) por \(a + b\), obtenemos: \[ \textsf{E}(\theta \mid y) = \frac{a}{a + b} \cdot \frac{a + b}{a + b + n} + \frac{y}{n} \cdot \frac{n}{a + b + n}. \]

Definiendo \(\omega = \frac{b}{b + n}\) y \(1-\omega = \frac{n}{b + n}\), podemos expresar: \[ \frac{a}{a + b} \cdot \frac{a + b}{a + b + n} = \omega \cdot \textsf{E}(\theta), \] y \[ \frac{y}{n} \cdot \frac{n}{a + b + n} = (1-\omega) \cdot \bar{y}. \]

Por lo tanto, la media posterior puede escribirse como: \[ \textsf{E}(\theta \mid y) = \omega \cdot \textsf{E}(\theta) + (1-\omega) \cdot \bar{y}. \]

Esto demuestra que la media posterior de \(\theta\) es un promedio ponderado entre la media previa de \(\theta\) y el promedio de éxitos \(\bar{y}\), donde los pesos están determinados por \(\omega = \frac{b}{b+n}\) y \(1-\omega = \frac{n}{b+n}\).

  1. Muestre que las distribuciones Bernoulli, Binomial, Multinomial, Poisson, Exponencial, Beta, Gamma, y Normal hacen parte de la familia exponencial.

La familia exponencial incluye distribuciones cuya función de probabilidad o densidad puede escribirse en la forma general: \[ p(x \mid \eta) = h(x) \exp\left(\eta^\top t(x) - a(\eta)\right), \] donde \(x\) es la variable aleatoria, \(\eta\) es el parámetro natural (o canónico), \(t(x)\) es el estadístico suficiente de los datos, \(a(\eta)\) es la función de normalización que asegura que \(p(x \mid \eta)\) sea válida, y \(h(x)\) es una función base que no depende de \(\eta\). A continuación, verificamos cómo las distribuciones mencionadas pertenecen a esta familia.

Distribución Bernoulli

La función de probabilidad para la distribución Bernoulli es: \[ p(x \mid \theta) = \theta^x (1-\theta)^{1-x}, \quad x \in \{0, 1\}. \] Reescribimos utilizando el parámetro natural \(\eta = \log\left(\frac{\theta}{1-\theta}\right)\), de donde: \[ \theta = \frac{e^\eta}{1+e^\eta}, \quad 1-\theta = \frac{1}{1+e^\eta}. \] Sustituyendo, tenemos: \[ p(x \mid \eta) = \exp\left(x\eta - \log(1+e^\eta)\right). \] Aquí, \(t(x) = x\), \(a(\eta) = \log(1+e^\eta)\), y \(h(x) = 1\), lo que confirma que Bernoulli pertenece a la familia exponencial.

Distribución Binomial

La función de probabilidad de la Binomial es: \[ p(x \mid \theta) = \binom{n}{x} \theta^x (1-\theta)^{n-x}, \quad x \in \{0, 1, \ldots, n\}. \] Con el parámetro natural \(\eta = \log\left(\frac{\theta}{1-\theta}\)), obtenemos: \[ p(x \mid \eta) = \binom{n}{x} \exp\left(x\eta - n\log(1+e^\eta)\right). \] Aquí, \(t(x) = x\), \(a(\eta) = n\log(1+e^\eta)\), y \(h(x) = \binom{n}{x}\), mostrando que Binomial pertenece a la familia exponencial.

Distribución Multinomial

La función de probabilidad para la Multinomial es: \[ p(x_1, \ldots, x_k \mid \theta_1, \ldots, \theta_k) = \frac{n!}{x_1! \cdots x_k!} \prod_{i=1}^k \theta_i^{x_i}, \quad \sum_{i=1}^k \theta_i = 1. \] Utilizando parámetros naturales \(\eta_i = \log(\theta_i)\), con la restricción \(\sum_{i=1}^k e^{\eta_i} = 1\), podemos escribir: \[ p(x \mid \eta) = \frac{n!}{x_1! \cdots x_k!} \exp\left(\sum_{i=1}^k x_i \eta_i - n\log\left(\sum_{i=1}^k e^{\eta_i}\right)\right). \] Aquí, \(t(x) = (x_1, \ldots, x_k)\), \(a(\eta) = n\log\left(\sum_{i=1}^k e^{\eta_i}\right)\), y \(h(x) = \frac{n!}{x_1! \cdots x_k!}\), confirmando que pertenece a la familia exponencial.

Distribución Poisson

La función de probabilidad de la Poisson es: \[ p(x \mid \theta) = \frac{\theta^x e^{-\theta}}{x!}, \quad x \in \{0, 1, 2, \ldots\}. \] Usando el parámetro natural \(\eta = \log(\theta)\), tenemos: \[ p(x \mid \eta) = \frac{1}{x!} \exp\left(x\eta - e^\eta\right). \] Aquí, \(t(x) = x\), \(a(\eta) = e^\eta\), y \(h(x) = \frac{1}{x!}\), mostrando que pertenece a la familia exponencial.

Distribución Exponencial

La función de densidad de la Exponencial es: \[ p(x \mid \theta) = \theta e^{-\theta x}, \quad x \geq 0. \] Usando el parámetro natural \(\eta = -\theta\), tenemos: \[ p(x \mid \eta) = \exp\left(\eta x - \log(-\eta)\right), \quad \eta < 0. \] Aquí, \(t(x) = x\), \(a(\eta) = -\log(-\eta)\), y \(h(x) = 1\), mostrando que pertenece a la familia exponencial.

Distribución Beta

La función de densidad de la Beta es: \[ p(x \mid \alpha, \beta) = \frac{x^{\alpha-1}(1-x)^{\beta-1}}{B(\alpha, \beta)}, \quad x \in (0, 1). \] Usando los parámetros naturales \(\eta_1 = \alpha-1\) y \(\eta_2 = \beta-1\), tenemos: \[ p(x \mid \eta) = \exp\left(\eta_1 \log(x) + \eta_2 \log(1-x) - \log(B(\alpha, \beta))\right). \] Aquí, \(t(x) = (\log(x), \log(1-x))\), \(a(\eta) = \log(B(\alpha, \beta))\), y \(h(x) = 1\), mostrando que pertenece a la familia exponencial.

Distribución Gamma

La función de densidad de la Gamma es: \[ p(x \mid \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha-1} e^{-\beta x}, \quad x > 0. \] Usando los parámetros naturales \(\eta_1 = \alpha - 1\) y \(\eta_2 = -\beta\), tenemos: \[ p(x \mid \eta) = \exp\left(\eta_1 \log(x) + \eta_2 x - \log\Gamma(\alpha) - \alpha \log(\beta)\right). \] Aquí, \(t(x) = (\log(x), x)\), \(a(\eta) = \log\Gamma(\alpha) + \alpha \log(\beta)\), y \(h(x) = 1\), confirmando que pertenece a la familia exponencial.

Distribución Normal

La función de densidad de la Normal es: \[ p(x \mid \theta, \sigma^2) = \frac{1}{\sqrt{2\pi \sigma^2}} \exp\left(-\frac{(x-\theta)^2}{2\sigma^2}\right). \] Usando los parámetros naturales \(\eta_1 = \frac{\theta}{\sigma^2}\) y \(\eta_2 = -\frac{1}{2\sigma^2}\), podemos escribir: \[ p(x \mid \eta) = \exp\left(\eta_1 x + \eta_2 x^2 - \frac{\theta^2}{2\sigma^2} - \frac{1}{2}\log(2\pi\sigma^2)\right). \] Aquí, \(t(x) = (x, x^2)\), \(a(\eta) = \frac{\theta^2}{2\sigma^2} + \frac{1}{2}\log(2\pi\sigma^2)\), y \(h(x) = 1\), mostrando que pertenece a la familia exponencial.

  1. Jeffreys (H. Jeffreys (1961). Theory of Probability. Oxford Univ. Press.) sugirió una regla para generar una distribución previa de un parámetro \(\theta\) asociado con la distribución muestral \(p(y\mid\theta)\). La distribución previa de Jeffreys es de la forma \(p_J(\theta)\propto\sqrt{I(\theta)}\), donde \[ I(\theta) = -\textsf{E}_{y\mid\theta}\left( \frac{\text{d}^2}{\text{d}\theta^2}\log p(y\mid\theta) \right) \] es la información esperada de Fisher.
  1. Sea \(y\mid\theta\sim\textsf{Bin}(n,\theta)\). Muestre que la distribución previa de Jeffreys para esta distribución muestral es \[ p_J(\theta)\propto \theta^{-\frac12} (1-\theta)^{-\frac12}\,. \]

Si \(y\mid\theta\sim\textsf{Bin}(n,\theta)\), entonces \(p(y\mid\theta) = \binom{n}{y}\theta^y(1-\theta)^{n-y}\), y por lo tanto \(\log p(y\mid\theta) = \log\binom{n}{y} + y\log\theta + (n-y)\log(1-\theta)\), de donde \[ \frac{\text{d}}{\text{d}\theta}\log p(y\mid\theta) = \frac{y}{\theta} - \frac{n-y}{1-\theta} \qquad\text{y}\qquad \frac{\text{d}^2}{\text{d}\theta^2}\log p(y\mid\theta)=-\left(\frac{y}{\theta^2}+\frac{n-y}{(1-\theta)^2}\right)\,. \] Así, \[ I(\theta) = -\textsf{E}_{y\mid\theta}\left( \frac{\text{d}^2}{\text{d}\theta^2}\log p(y\mid\theta) \right) = -\textsf{E}_{y\mid\theta}\left( -\left(\frac{y}{\theta^2}+\frac{n-y}{(1-\theta)^2}\right) \right) = \frac{\textsf{E}_{y\mid\theta}(y)}{\theta^2}+\frac{n-\textsf{E}_{y\mid\theta}(y)}{(1-\theta)^2} = \frac{n\theta}{\theta^2}+\frac{n-n\theta}{(1-\theta)^2} = \frac{n}{\theta(1-\theta)}\,. \] Entonces, \[ p_J(\theta) \propto \sqrt{I(\theta)} = \sqrt{\frac{n}{\theta(1-\theta)}} \propto \frac{1}{\sqrt{\theta(1-\theta)}} = \theta^{-1/2}(1-\theta)^{-1/2}\,. \]

  1. Reparametrice la distribución Binomial con \(\psi = \textsf{logit}(\theta)\), de forma que \[p(y\mid\psi) \propto e^{\psi y}(1+e^\psi)^{-n}\,.\] Obtenga la distribución previa de Jeffreys para esta distribución muestral.

Si \(\psi = \textsf{logit}(\theta)=\log\frac{\theta}{1-\theta}\), entonces \(\theta = \frac{e^\psi}{1+e^\psi}\). Por lo tanto, \[ p(y\mid\psi) = \binom{n}{y}\left( \frac{e^\psi}{1+e^\psi} \right)^y\left(1-\frac{e^\psi}{1+e^\psi}\right)^{n-y} = \binom{n}{y}\, \frac{e^{\psi y}}{(1+e^\psi)^y}\,\frac{1}{(1+e^\psi)^{n-y}} = \binom{n}{y} e^{\psi y}(1+e^\psi)^{-n}\,. \] Luego, \(\log p(y\mid\psi) = \log\binom{n}{y} + \psi y -n\log(1+e^\psi)\), de donde \[ \frac{\text{d}}{\text{d}\psi}\log p(y\mid\psi) = y - n\,\frac{e^\psi}{1+e^\psi} \qquad\text{y}\qquad \frac{\text{d}^2}{\text{d}\psi^2}\log p(y\mid\psi)=-n\,\frac{e^\psi}{(1+e^\psi)^2}\,. \] Así, \[ I(\psi) = -\textsf{E}_{y\mid\psi}\left( \frac{\text{d}^2}{\text{d}\psi^2}\log p(y\mid\psi) \right) = -\textsf{E}_{y\mid\psi}\left( -n\,\frac{e^\psi}{(1+e^\psi)^2} \right) = n\,\frac{e^\psi}{(1+e^\psi)^2}\,. \] Entonces, \[ p_J(\psi) \propto \sqrt{I(\psi)} = \sqrt{n\,\frac{e^\psi}{(1+e^\psi)^2}} \propto \frac{e^{\psi/2}}{1+e^\psi}\,. \]

  1. Tome la distribución previa de la parte a. y aplique la fórmula del cambio de variables para obtener la densidad previa de \(\psi\). Esta densidad debe coincidir con la obtenida en el inciso b.. Esta propiedad de invarianza bajo reparametrización es la característica fundamental de la previa de Jeffreys.

Si \(p_J(\theta) \propto \theta^{-1/2}(1-\theta)^{1-/2}\) y \(\psi = \textsf{logit}(\theta)=\log\frac{\theta}{1-\theta}\), entonces aplicando la fórmula del cambio de variables (e.g., https://online.stat.psu.edu/stat414/lesson/22/22.2) se tiene que \(\theta = \frac{e^\psi}{1+e^\psi}\) y \(\frac{\textsf{d}\theta}{\textsf{d}\psi}=e^{-\psi}/(1+e^{-\psi})^2\), de donde \[ p_J(\psi) = p_J(\theta)\,\Big|\frac{\textsf{d}\theta}{\textsf{d}\psi}\Big|= \left( \frac{e^\psi}{1+e^\psi} \right)^{-1/2}\left(1- \frac{e^\psi}{1+e^\psi}\right)^{-1/2}\,\frac{e^{-\psi}}{(1+e^{-\psi})^2} = \frac{e^{\psi/2}}{1+e^\psi}\,. \]

  1. El estimador óptimo del parámetro \(\theta\in \Theta\subset\mathbb{R}\) de acuerdo con la regla de Bayes se define como el estimador \(\hat\theta=\hat\theta(\boldsymbol{y})\) que minimiza la perdida esperada posterior dada por \[ \textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta)) = \int_{\Theta} L(\theta,\hat\theta)\, p(\theta\mid\boldsymbol{y})\, \textsf{d}\theta\,, \] donde \(L(\theta,\hat\theta)\) es una función de perdida (costo que conlleva estimar \(\theta\) por medio de \(\hat\theta\)) y \(\boldsymbol{y}=(y_1,\ldots, y_n)\) es un conjunto de datos observado.
  1. Muestre que si \(L(\theta,\hat\theta) = (\theta - \hat\theta)^2\), entonces el estimador óptimo de acuerdo con la regla de Bayes es la media posterior \(\hat\theta=\textsf{E}(\theta\mid\boldsymbol{y})\).

Diferenciado \(\textsf{E}_{\theta\mid\boldsymbol{y}}((\theta-\hat\theta)^2)\) respecto a \(\hat\theta\), se tiene que \[ \frac{\textsf{d}}{\textsf{d}\hat\theta}\textsf{E}_{\theta\mid\boldsymbol{y}}\left((\theta-\hat\theta)^2\right)=\frac{\textsf{d}}{\textsf{d}\hat\theta}\textsf{E}_{\theta\mid\boldsymbol{y}}\left(\theta^2-2\theta\hat\theta+\hat{\theta}^2\right) = -2\textsf{E}_{\theta\mid\boldsymbol{y}}(\theta) + 2\hat\theta\, \] y por lo tanto el punto crítico correspondiente es \(\hat\theta = \textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\). Así, por el criterio de la segunda derivada, como \[ \frac{\textsf{d}^2}{\textsf{d}\hat\theta^2}\textsf{E}_{\theta\mid\boldsymbol{y}}\left((\theta-\hat\theta)^2\right)= 2 > 0 \] para todo \(\hat\theta\), entonces en efecto la media posterior \(\hat\theta = \textsf{E}_{\theta\mid\boldsymbol{y}}(\theta)\) minimiza \(\textsf{E}_{\theta\mid\boldsymbol{y}}((\theta-\hat\theta)^2)\).

  1. Muestre que si \(L(\theta,\hat\theta) = |\theta - \hat\theta|\), entonces el estimador óptimo de acuerdo con la regla de Bayes es la mediana posterior \(\hat\theta=(\theta\mid\boldsymbol{y})_{0.5}\).

Diferenciado \(\textsf{E}_{\theta\mid\boldsymbol{y}}(|\theta-\hat\theta|)\) respecto a \(\hat\theta\) aplicando la regla de Leibniz (e.g., https://mathworld.wolfram.com/LeibnizIntegralRule.html), se tiene que \[ \frac{\textsf{d}}{\textsf{d}\hat\theta}\textsf{E}_{\theta\mid\boldsymbol{y}}\left(|\theta-\hat\theta|\right)=\frac{\textsf{d}}{\textsf{d}\hat\theta}\int_\Theta |\theta-\hat\theta|\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta = \frac{\textsf{d}}{\textsf{d}\hat\theta}\left( \int_{-\infty}^{\hat\theta} (\hat\theta-\theta)\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta + \int_{\hat\theta}^\infty (\theta-\hat\theta)\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta \right) = \int_{-\infty}^{\hat\theta}p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta - \int_{\hat\theta}^\infty \,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta\,, \] y por lo tanto, el punto crítico satisface que \[ \int_{-\infty}^{\hat\theta}p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta = \int_{\hat\theta}^\infty \hat\theta\,p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta \quad\text{de donde}\quad 2\int_{-\infty}^{\hat\theta}p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta = \int_{-\infty}^{\hat\theta}p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta+\int_{\hat\theta}^\infty p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta = 1\,, \] esto es, \(\int_{-\infty}^{\hat\theta}p(\theta\mid\boldsymbol{y})\,\textsf{d}\theta =1/2\). Luego, el punto crítico correspondiente es \(\hat\theta = (\theta\mid\boldsymbol{y})_{0.5}\). Así, por el criterio de la segunda derivada, como \[ \frac{\textsf{d}^2}{\textsf{d}\hat\theta^2}\textsf{E}_{\theta\mid\boldsymbol{y}}\left(|\theta-\hat\theta|\right)= 2 > 0 \] para todo \(\hat\theta\), entonces en efecto la mediana posterior \(\hat\theta = (\theta\mid\boldsymbol{y})_{0.5}\) minimiza \(\textsf{E}_{\theta\mid\boldsymbol{y}}(|\theta-\hat\theta|)\).

  1. El riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) se define como \[ R_{\textsf{F}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}\mid\theta}(L(\theta,\hat\theta)) = \int_{\mathcal{Y}} L(\theta,\hat\theta)\, p(\boldsymbol{y}\mid\theta)\,\textsf{d}\boldsymbol{y}\,, \] i.e., el valor medio de la función de perdida \(L(\theta,\hat\theta)\) a través de todos los valores de \(\boldsymbol{y} \in \mathcal{Y}\). De otra parte, el riesgo Bayesiano \(R_{\textsf{B}}(\theta,\hat\theta)\) se define como \[ 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\,, \] i.e., el valor medio del riesgo frecuentista \(R_{\textsf{F}}(\theta,\hat\theta)\) a priori a través de todos los valores de \(\theta\in\Theta\). Muestre que \[ R_{\textsf{B}}(\theta,\hat\theta) = \textsf{E}_{\boldsymbol{y}}(\textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta)))\,, \] donde \(\textsf{E}_{\boldsymbol{y}}(\cdot)\) denota el valor esperado respecto a la distribución marginal \(p(\boldsymbol{y})\).

Se tiene que \[ \begin{align*} R_{\textsf{B}}(\theta,\hat\theta) &= \int_{\Theta} R_{\textsf{F}}(\theta,\hat\theta)\, p(\theta)\,\textsf{d}\theta \\ &= \int_{\Theta} \left( \int_{\mathcal{Y}} L(\theta,\hat\theta)\, p(\boldsymbol{y}\mid\theta)\,\textsf{d}\boldsymbol{y} \right)\, p(\theta)\,\textsf{d}\theta \\ &= \int_{\Theta} \left( \int_{\mathcal{Y}} L(\theta,\hat\theta)\, p(\boldsymbol{y}\mid\theta)\, p(\theta)\,\textsf{d}\boldsymbol{y} \right)\,\textsf{d}\theta \\ &= \int_{\Theta} \left( \int_{\mathcal{Y}} L(\theta,\hat\theta)\, p(\theta\mid\boldsymbol{y})\, p(\boldsymbol{y})\,\textsf{d}\boldsymbol{y} \right)\,\textsf{d}\theta \\ &= \int_{\mathcal{Y}} \left( \int_{\Theta} L(\theta,\hat\theta)\, p(\theta\mid\boldsymbol{y})\, p(\boldsymbol{y})\, \textsf{d}\theta \right)\textsf{d}\boldsymbol{y} \\ &= \int_{\mathcal{Y}} \left( \int_{\Theta} L(\theta,\hat\theta)\, p(\theta\mid\boldsymbol{y})\, \textsf{d}\theta \right)\, p(\boldsymbol{y})\,\textsf{d}\boldsymbol{y} \\ &= \int_{\mathcal{Y}} \textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta))\, p(\boldsymbol{y})\,\textsf{d}\boldsymbol{y} \\ &= \textsf{E}_{\boldsymbol{y}}(\textsf{E}_{\theta\mid\boldsymbol{y}}(L(\theta,\hat\theta)))\,. \end{align*} \]

  1. Suponga que \(y_i\mid\boldsymbol{\theta}\stackrel{\text{iid}}{\sim} p(y_i\mid\boldsymbol{\theta})\), para \(i=1,\ldots,n\), con \(\boldsymbol{\theta}\sim p(\boldsymbol{\theta})\), donde \(\boldsymbol{\theta}=(\theta_1,\ldots,\theta_k)\) y \(\boldsymbol{y} = (y_1\ldots,y_n)\).

Sea \(\boldsymbol{\theta}_\text{MAP}\) el máximo a posteriori (MAP, por sus siglas en inglés), el cual se define como el valor de \(\boldsymbol{\theta}\) que maximiza la distribución posterior posterior, a saber, \[ \boldsymbol{\theta}_\text{MAP} = \textsf{arg max}_{\boldsymbol{\theta}}\,\log p(\boldsymbol{\theta}\mid\boldsymbol{y}) = \textsf{arg max}_{\boldsymbol{\theta}}\,\left(\log p(\boldsymbol{y}\mid\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})\right)\,, \] y además, sea \(\mathbf\Sigma_\text{MAP} = -\mathbf{H}^{-1}\), donde \(\mathbf{H}\) es la matriz Hessiana cuya \((i,j)\)-ésima entrada está dada por \[ H_{i,j} = \frac{\partial^2}{\partial\theta_i\,\partial\theta_j}\,\log p (\boldsymbol{\theta}\mid\boldsymbol{y}) \Bigg|_{\boldsymbol{\theta}=\boldsymbol{\theta}_\text{MAP}} = \frac{\partial^2}{\partial\theta_i\,\partial\theta_j}\,\left(\log p(\boldsymbol{y}\mid\boldsymbol{\theta}) + \log p(\boldsymbol{\theta})\right)\Bigg|_{\boldsymbol{\theta}=\boldsymbol{\theta}_\text{MAP}}\,. \] Se observa que la constante de normalización \(p(\boldsymbol{y})\) no depende de \(\boldsymbol{\theta}\), y por lo tanto, no interfiere en la maximización de \(\log p(\boldsymbol{\theta}\mid\boldsymbol{y})\). El Teorema del Límite Central Bayesiano (BCLT, por sus siglas en inglés) indica que \(\boldsymbol{\theta}\mid \boldsymbol{y} \approx \textsf{N}_k(\boldsymbol{\theta}_\text{MAP},\mathbf\Sigma_\text{MAP})\), cuando \(n\rightarrow\infty\).

Usando los datos del caso de victimas de violencia sexual, aproxime la distribución posterior de \(\theta\) por medio del BCLT. Dibuje la distribución exacta y la aproximada en un mismo gráfico.

El modelo se define como: \[ y \mid \theta \sim \textsf{Binomial}(n, \theta), \quad \theta \sim \textsf{Beta}(a, b). \]

En este caso se tiene que \(a = 1\), \(b = 1\), \(n = 80\), \(s = 69\), y por lo tanto la distribución posterior de \(\theta\) es: \[ p(\theta \mid y) \propto p(y \mid \theta)p(\theta) = \binom{n}{y} \theta^y (1-\theta)^{n-y} \cdot \frac{\theta^{a-1} (1-\theta)^{b-1}}{B(a, b)}. \]

Al simplificar, se observa que: \[ p(\theta \mid y) \propto \theta^{y+a-1}(1-\theta)^{n-y+b-1}. \] Esto corresponde a una distribución \(\textsf{Beta}(a + s, b + n - s)\), es decir, \(\theta \mid y \sim \textsf{Beta}(70, 12)\).

El MAP de \(\theta\) es el valor que maximiza \(\log p(\theta \mid y)\), que equivale a maximizar: \[ \log p(\theta \mid y) = (s + a - 1)\log(\theta) + (n - s + b - 1)\log(1-\theta). \] Derivamos con respecto a \(\theta\) y encontramos el máximo: \[ \frac{\partial}{\partial \theta} \log p(\theta \mid y) = \frac{s + a - 1}{\theta} - \frac{n - s + b - 1}{1-\theta}. \] Igualando a cero: \[ \frac{s + a - 1}{\theta} = \frac{n - s + b - 1}{1-\theta}. \] Resolviendo para \(\theta\), obtenemos: \[ \theta_\text{MAP} = \frac{s + a - 1}{n + a + b - 2}. \] Sustituyendo los valores, \(s = 69\), \(a = 1\), \(b = 1\), \(n = 80\): \[ \theta_\text{MAP} = \frac{69 + 1 - 1}{80 + 1 + 1 - 2} = \frac{69}{80} = 0.8625. \]

El Hessiano de \(\log p(\theta \mid y)\) es: \[ H = \frac{\partial^2}{\partial \theta^2} \log p(\theta \mid y) = -\frac{s + a - 1}{\theta^2} - \frac{n - s + b - 1}{(1-\theta)^2}. \] Evaluamos en \(\theta_\text{MAP}\): \[ H = -\frac{69}{(0.8625)^2} - \frac{11}{(1-0.8625)^2}. \] Calculando: \[ H = -92.66 - 92.36 = -185.02. \]

La varianza aproximada es \(-1/H\), y por tanto: \[ \text{Var}(\theta \mid y) \approx \frac{1}{185.02} \approx 0.0054. \]

Por lo tanto, la distribución aproximada es: \[ \theta \mid y \approx \textsf{N}(0.8625, 0.0054). \]

Para graficar ambas distribuciones, generamos puntos de la posterior exacta \(\textsf{Beta}(70, 12)\) y de la aproximación normal \(\textsf{N}(0.8625, 0.0054)\).

# parámetros
a <- 1; b <- 1; n <- 80; s <- 69
theta_map <- s / n
var_map <- 1 / ((s / theta_map^2) + ((n - s) / (1 - theta_map)^2))
# puntos para graficar
theta <- seq(0.6, 1, length.out = 1000)
posterior_exact <- dbeta(theta, a + s, b + n - s)
posterior_approx <- dnorm(theta, mean = theta_map, sd = sqrt(var_map))
# gráfico
plot(theta, posterior_exact, type = "l", col = "blue", lwd = 2,
     ylab = "Densidad", xlab = expression(theta),
     main = "Distribución exacta y aproximada")
lines(theta, posterior_approx, col = "red", lwd = 2, lty = 2)
legend("topleft", legend = c("Exacta (Beta)", "Aproximada (Normal)"),
       col = c("blue", "red"), lwd = 2, lty = c(1, 2), bty = "n")