Introducción

Ejercicio 1: Se desea estimar la probabilidad, \(\theta\), de un evento, a partir del resultado de una sucesión de n ensayos Bernoulli, esto es, \(y_{1},...,y_{n}\) que son iguales a 1 si ocurre el evento (éxito) y cero si no ocurre. Sea y el número total de éxitos en el muestra de n ensayos. El modelo binomial simple es el siguiente:

\[ \begin{aligned} p\left(y|\theta\right)&=Bin\left(y|n , \theta\right)={n \choose y} \theta^{y}\left(1- \theta\right)^{n-y}\\ p\left(\theta\right)& \sim Uniforme\\ p\left(\theta\right | y)&= \propto \theta^{y}\left(1-\theta\right)^{n-y} \end{aligned} \]

La distribución no normalizada tiene un kernel equivalente a una distribución Beta, es decir:

\[\theta|y \propto Beta \left(y+1 , n-y+1\right )\]

Ejercicio 2: Sea \(y_{1},...,y_{n}\) kis regustris dek oesi de un objeto medido n veces. Sea \(\theta = N(\mu, \sigma^{2})\) el peso verdadero y la varianza de medición del peso respectivamente. Sea \(\tilde{y}\) el peso del objeto para una nueva pesada. Encuentre la distribucón predicitiva a posterior.

\[ \begin{aligned} p(\tilde{y} | y) &= \int p(\tilde{y},\theta|y) d\theta \\ &= \int p(\tilde{y},\theta|y) p(\theta|y)d\theta \\ &= \int p(\tilde{y},\theta) p(\theta|y) d\theta \\ \end{aligned} \]

Ejercicio 3: Sea X de dimensión 2 con \(X= (Y,Z)\) donde Y y Z son escalares. La función de densidad de X viene dada por. \[ \begin{aligned} f\left(x\right) = f\left(y,z\right) = \dfrac{1}{\sqrt{2\pi}\sigma^{2}}exp \lbrace -\dfrac{1}{2\sigma^{2}}[-(y-\theta_{y})^2 -(x-\theta_{z})^2] \rbrace, -\infty < y < \infty , -\infty < z < \infty \\ &\rightarrow X \sim N((\theta_{y},\theta_{x}),\sigma^{2}I) \\ Y& = log(U) \\ Z& =log(V) \end{aligned} \]

Por el teorema de la transformación, la densidad \(g\left(u,v\right)\) viene dada por:

\[ \begin{aligned} g\left(u,v\right) &= f\left(y\left(u,v\right),z\left(u,v\right)\right)|J\left(y,z\right)\rightarrow\left(u,v\right)| \\ &=\dfrac{1}{2\pi\sigma^{2}}exp \lbrace-\dfrac{1}{2\sigma^{2}}[(log u - \theta_{y})^2+(log v - \theta_{z})^2] \rbrace \times |J\left(y,z\right)\rightarrow\left(u,v\right)|\\ \end{aligned} \]

Luego el Jacobiano es:

\[ \begin{aligned} \left|J\right|&=\left| \begin{array}{cc} \dfrac{\partial y}{\partial u} & \dfrac{\partial y}{\partial v} \\ \dfrac{\partial z}{\partial u} & \dfrac{\partial z}{\partial v} \end{array} \right| \\ &=\left| \begin{array}{cc} \dfrac{1}{u} & 0\\ 0 & \dfrac{1}{u} \end{array} \right| \\ &= \dfrac{1}{uv} \\ &\rightarrow g\left(u,v\right)+\dfrac{1}{2\pi\sigma^{2}}exp\lbrace [-\dfrac{1}{2\sigma^{2}}[(log u - \theta_{y})^2+(log v - \theta_{z})^2]\rbrace , \text{ } 0 < u < \infty , \text{ }0 < v < \infty \end{aligned} \]

Ejercicio 4: Considere el modelo de probabilidad de Poisson \(p\left(y| \theta\right)\), esto es una variable discreta, y que toma valores \(\left\lbrace 0, 1, \cdots \right\rbrace\) con probabibilidad de que \(y=k\), igual a \(exp \left(- \theta\right)\theta^{k}/k!\)

Si se obtienen n observaciones, \(y_{1} , y_{2}, \cdots, y_{n}\) y se asume independencia, la verosimilitud estará dada por

\[ \begin{aligned} &= exp\left(- \theta\right)\dfrac{\theta^{y_{1}}}{y_{1}!} \times exp\left(- \theta\right)\dfrac{\theta^{y_{2}}}{y_{2}!} \times \cdots \times exp\left(- \theta\right)\dfrac{\theta^{y_{n}}}{y_{n}!} \times \\ &= exp\left(- n\theta\right)\dfrac{\theta^{\sum y_{i}}}{\prod y_{i}!} \\ &= exp\left(- n\theta\right)\dfrac{exp({\sum y_{i}ln \theta})}{\prod y_{i}!} \\ & \Rightarrow p(y|\theta) \propto exp\left(- n\theta\right)exp({\sum y_{i}ln \theta}) \end{aligned} \] Examinando la forma de esta función es fácil ver que una distribución con función de densidad proporcional a \[ exp\left(-\theta \tau_{0}\right)exp \left(\tau_{1} \log \theta\right) \] es una distribución conjugada, y en particular, haciendo \[ \begin{aligned} \tau_{0}&=\beta \\ \tau_{0}&=\alpha - 1 \\ \end{aligned} \]

Se puede reconocer como una distribución previa a una distribución Gamma, con parámetros \(\alpha , \beta\). Consecuentemente, la distribución posterior es tambien una Gamma con parámetros

\[ \begin{aligned} \alpha^{\ast}&=\alpha + \sum y_{i}\\ \beta^{\ast}&= \beta + n \end{aligned} \] Ejercicio 5: Si \(Y\sim Bin(n,\theta)\), su log verosimil es:

\[ \begin{aligned} Bin(n,\theta) &= {n\choose y}\theta^y(1-\theta)^{n-y} \Rightarrow\\ \lg p(y|\theta)&=\lg{n\choose y}+\lg \theta^y+\lg (1-\theta)^{n-y}= cte+y\lg \theta+(n-y)\lg(1-\theta)\\ \dfrac{\partial \lg p(y|\theta)}{\partial \theta}&=\dfrac{y}{\theta}-\dfrac{(n-y)}{1-\theta}=\dfrac{y-y\theta-n\theta+y\theta}{\theta(1-\theta)}=\dfrac{y-n\theta}{\theta(1-\theta)}\Rightarrow\\ J(\theta)&= E \left[\left(\dfrac{\partial \lg p(y|\theta)}{\partial \theta}\right)^2|\theta\right]= E\left[\dfrac{(y-n\theta)^2}{\theta^2(1-\theta^2)}\right]\\ J(\theta)&=\dfrac{1}{\theta^2(1-\theta^2)}E[(y-E(y))^2]=\dfrac{1}{\theta^2(1-\theta^2}Var(y)\\ J(\theta)&=\dfrac{n}{\theta(1-\theta)} \end{aligned} \]

Ejercicio 6: El conjunto de datos normales, normaldata, está relacionado con el famoso experimento de Michelson-Morlay que abrió el comino a la teoría de la relatividad de Einstein en 1887. El experimento pretendía detector el “flujo del éter” y, por tanto, la existencia del éter, esto medio teórico que los fisicos postulaban en osa época como necesario para la transmisión de la luz. El método de medición de Michelson consistía en medir la diferencia de velocidad de dos haces de luz que recorrían la misma distancia en dos direcciones ortogonales.

Como suele ocurrir con física, la medición se hacía por interferometría y las diferencias en el tiempo de viaje se iuferían a partir del desplazamiento de las franjas del espectro luminoso, Sin embargo, el experimento. produjo mediciones muy pequeñas que no fueron concluyentes para la detección del éter. Los experimentos. posteriores intentaron alcanzar tuna mayor precisión, como el realizado por Illingworth en 1927, utilizado aquí como datos normales, sólo para obtener lhmites superiores cada vez más pequeños de la velocidad del éter, Aunque el conjunto de datos original está disponible en R como morley, las entradas se aproximan al múltiplo de diez más cercano y, por tanto, son dificiles de analizar como observaciones normales. Los 64 puntos de. datos en normaldata están asociados a números de sesión (primera columna), que corresponden a diferentes. momentos del día, y los valores de la segunda columna representan el desplazamiento medio de la franja debida a la orientación tomado sobre diez mediciones realizadas por Alinguorth, que asumió un modelo de. error normal. La figura siguiento produce un histograma de los datos medianto los sencillos comandos de R.

Este histograma parece compatible con ula distribución simétrica unimodal como la distribución normal. Como se muestra en la Figura un gráfico qq.

library(bayess)
## Loading required package: MASS
## Loading required package: mnormt
## Loading required package: gplots
## 
## Attaching package: 'gplots'
## The following object is masked from 'package:stats':
## 
##     lowess
## Loading required package: combinat
## 
## Attaching package: 'combinat'
## The following object is masked from 'package:utils':
## 
##     combn
data("normaldata")
shift = normaldata[,2]
hist(shift , nclass = 10, col = "skyblue",prob = TRUE, main="")

library(bayess)
qqnorm((shift - mean(shift))/sd(shift),pch = 19, col="skyblue")
abline(a=0,b=1,lty = 2, col="black",lwd=2)

El ajuste de \(N(\mu,\sigma^2)\) puede no ser perfecto, sin embargo, debido a(a) una posible bimodalidad del histograma y (b) posibles valores atípicos.

n = length(shift)
mmu = sum(shift)/(n+1)
mmu
## [1] -0.01461538
vmu = 0.75^2/(n+1);vmu
## [1] 0.008653846

Ejercicio 7: Como ilustración, consideremos el caso más sencillo de la distribución normal con varianza conocida, \(N(\mu,\sigma^2)\). Si la distribución a priori de \(\mu\), \(p(\mu)\), es la normal \(N(0,\sigma^2)\), la distribución posterior se obtiene fácilmente mediante la teoría de Bayes.

\[ \begin{aligned} p(y|\theta) \propto & \ p(u)p(y|\theta) \\ \propto \ & \dfrac{1}{\sqrt{2\pi}\sigma^{2}}exp \lbrace -\dfrac{1}{2\sigma^{2}} \mu^2\rbrace \times \dfrac{1}{\sqrt2\pi\sigma^{2}} exp \lbrace \dfrac{1}{2\sigma^2}(y_1 - \mu)^2 \rbrace \times \cdots \times \dfrac{1}{\sqrt{2\pi}\sigma^{2}} exp \lbrace -\dfrac{1}{2\sigma^{2}}(y_n - \mu)^2 \rbrace \\ \propto \ & \dfrac{1}{\sqrt{2\pi}\sigma^{2}}exp \lbrace -\dfrac{1}{2\sigma^{2}} \mu^2\rbrace \left(\dfrac{1}{\sqrt{2\pi}\sigma^{2}}\right)^{n} exp \lbrace -\dfrac{1}{2\sigma^{2}} \sum_i ^n (y_i - \mu)^2\rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} \mu^2 - \dfrac{1}{2\sigma^{2}} \sum_{i=1}^{n} [(y_i-\bar{y})-(\mu-\bar{y})]^2\rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} \mu^2 -\dfrac{1}{2\sigma^{2}} [\sum_{i=1}^{n}(y_i-\bar{y})^2 - 2(\mu-\bar{y})\sum_{i=1}^n(y_i - \bar{y})+n(\mu-\bar{y})^2] \rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} \mu^2 -\dfrac{1}{2\sigma^{2}} \left[S_{y}^2+n(\mu-\bar{y})^2]\right \rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} S_y^{2} \rbrace \times exp \lbrace -\dfrac{1}{2\sigma^{2}} [\mu^2 + n\mu^2 -2n\mu^2+n\bar{y}^2] \rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} [(n+1)\mu^2-2n\mu\bar{y}+n\bar{y}^2]\rbrace \\ \propto \ & exp \lbrace -\dfrac{1}{2\sigma^{2}} (n+1)[\mu^2-\dfrac{2n\mu\bar{y}}{(n+1)}+\dfrac{n\bar{y}^2}{(n+1)}] \rbrace \\ \propto \ & exp \left\lbrace -\dfrac{1}{2\sigma^{2}}(n+1) \left[\left(\mu - \dfrac{n \bar{x}}{(n+1)}\right)^2 + \dfrac{n \bar{x}}{(n+1)} - \dfrac{(n \bar{x})^{2}}{(n+1)^{2}} \right]\right \rbrace \\ \end{aligned} \]

Lo que significa que esta distribución posterior en \(\mu\) es una distribución normarl con media \(\dfrac{n\bar{y}}{(n+1)}\) y varianza \(\dfrac{\sigma^2}{(n+1)}\). La media (y la moda) de la posterior es, por tanto, diferente al estimador clásico \(\bar{X}\), lo que puede parecer una característica de este análisis bayesiano. La razón de la diferencia es que la información a prior de que \(\mu\) está lo suficientemente cerca de ceca es tenida en cuenta por la distribución posterior, que de este modo escoge la estimación original hacia cero.

Ejercicio 8 y 9: Siendo el caso de una distribución normal con varianza conocida bastante poco realista, considermos ahora el caso general de una muestra iid \(x= (x_1 ,\cdots , x_n)\) de la distribución normal \(N(\mu,\sigma^2)\) y \(\theta = (\mu,\sigma^2)\). Manteniendo la misma distribución a priori \(N(\mu,\sigma^2)\) sobre \(\mu\), que aparece entocnes como una distribución condicional de \(\mu\) dada \(\sigma^2\), es decir, se basa en la descomposición genérica. \[ \begin{aligned} p(\mu,\sigma^2) = p(\mu|\sigma^2)p(\sigma^2) \end{aligned} \]

Tenemos que introducir otra distribución a priori sobre \(\sigma^2\). Para simplificar los cáculos en esta fase inicial elegimos una distribución exponencial E(1) sobre \(\sigma^{-2}\).Esto significa que la variable aleatoria \(w=\sigma^{-2}\) se distribuye a partir de una distribución exponencial E(1), derivándose de la distribución sobre \(\sigma^2\) mediante la técnica habitual de cambio de variable,

\[ p(\sigma^2)=exp(-\sigma^{-2})\left|\dfrac{d\sigma^{-2}}{d\sigma^{2}}\right| = exp(-\sigma^{-2})(\sigma^2)^{-2} \] Esta distribución es un caso especial de una distribución gamma inversa, concretamente IG(1.1). La correspondiente desndiad posterior de \(\theta\) viene dada entonces por:

\[ \begin{aligned} p((\mu,\sigma^2)|y) & \propto p(\mu,\sigma^2) \ell ((\mu,\sigma^2)|y) \\ & \propto exp(-\sigma^{-2})(\sigma^2)^{-2} \dfrac{1}{\sigma}\left(-\dfrac{1}{2\sigma^2}\right)(\sigma^{-2})^{n/2}exp \left(-\dfrac{1}{2\sigma^2}[S_{y}^2+n(\mu-\bar{y})^2]\right) \\ & \propto (\sigma^{-2})^{2-1/2} exp \left[-\dfrac{1}{\sigma^2}(\mu^2+2)\right](\sigma^{-2})^{n/2} exp \left(-\dfrac{1}{2\sigma^{2}}\left[S_{y}^2+n(\mu-\bar{y})^2\right]\right) \\ & \propto (\sigma^2)^{-(n+5)/2} exp \left[-\dfrac{1}{2\sigma^{2}} \left\lbrace (n+1)\left(\mu- \dfrac{n\bar{y}}{n+1}\right)^2 + \dfrac{n\bar{y}^2}{n+1}- \dfrac{(n\bar{y})^2}{n+1} + \ 2 \ + S_{y}^2\right\rbrace\right] \\ & \propto (\sigma^2)^{-1/2} exp \left[-\dfrac{1}{2\sigma^{2}} \left\lbrace (n+1)\left(\mu-\dfrac{n\bar{y}}{n+1}\right)^2\right\rbrace\right](\sigma^2)^{-(n+2)/2-1} exp \left[-\dfrac{1}{2\sigma^{2}}(2+S_{y}^2)\right]\\ \end{aligned} \]

Por lo tanto, la posterior sobre \(\theta\) puede descomponerse como el producto de una distribución gamma inversa sobre \(\sigma^2\), \(\mathcal{IG}\left(\dfrac{n+2}{2},\dfrac{2+S_{y}^2}{2}\right)\) que es la distribución de la inversa de una variable aleatoria \(\gamma\) \(\mathcal{G}\left(\dfrac{n+2}{2},\dfrac{2+S_{y}^2}{2}\right)\)y, condicionada a \(\sigma^2\), una distribución normal sobre \(\mu\), \(N(\dfrac{n\bar{y}}{n+1},\dfrac{\sigma^2}{n+1})\). La interpretación de esta posterior es bastante similar al caso en que se conoce \(\sigma\), con la diferencia de que la variabilidad en \(\sigma\) induce más variabilidad en \(\mu\), siendo entonces la posterior marginal en \(\mu\) una distribución en t de Student.

\[ \mu | x \sim \mathcal{T}((n+2),n\bar{y}/(n+1),(2+S_{y}^2)/(n+1)(n+2)) \]

Con n+2 grados de libertad, un parámetro de localización proporcional; a \(\bar{y}\) y un parámetro de escala (casi) proporcional a S.

Para los datos normales, con una previa \(\mathcal{E}\S(1)\) sobre \(\sigma^{-2}\) es compatible con el valor observado en el experimento de Michelson-Morley, por lo que los parámetros de la distribución t sobre \(\mu\) son \(n=64\).

mtmu = sum(shift)/(n+1)
mtmu
## [1] -0.01461538
stmu = (2+(n-1)*var(shift))/((n+2)*(n+1)) 
stmu
## [1] 0.001084149

Comparamos la posterior resultante con la basada en la hipótesis \(\sigma = 0.75\) en la Figura.

library(mnormt)
curve(dmt(x, mean = mmu, S=stmu,df=n+2),col="skyblue",lwd = 2,xlab="x",ylab = "",xlim = c(-.5,.5))
curve(dnorm(x,mean=mmu,sd= sqrt(vmu)),col="black",lwd=2,add=TRUE,lty=2)

Aunque esto puede parecer contradictorio, en este mismo caso, la estimación de la varianza produce una reducción de la variabilidad de la distribución posterior \(\mu\). Esto se debe a que el valor postulado de \(\sigma^2\) es en realidad inapropiado para el experimento de Illingworth, siendo demasiado grande. Dado que la distribución posterior de \(\sigma^2\) e una distribución \(\mathcal{IG}(33,1.82)\) para los datos normales, la probabilidad de que \(\sigma\) sean tan grande como 0.75 puede evaluarse como:

diagmma = function(x,shape,scale){dgmma(1/x,shape,scale)/x^2} > curve(digmma(x,shape=33,scale = (1+(n+1)*var(shift))))
pgamma(1/(.75)^2,shape = 33,scale = (1+(n+1)*var(shift))/2)
## [1] 8.994529e-39

que muestra que 0.75 es bastante irreal, siendo diez veces mayor que la moda de la densidad posterior en \(\sigma^2\).

Estimaciones Bayesianas

Un concepto que está en la base del análisis byaesiano es que se debe proporcionar una evaluación inderencial condicionada al valor realizado en los datos. El análisis bayesiano de un significado probabilístico adecuado a este condicionamiento asignado a \(\theta\) una distribución de probabilidad. Una vez seleccionada la distribución a priori, la inferencia bayesiana está formalmente “terminada” ; es decir, está completamente determinada, yya que los procedimientos de estimación, prueba y evaluación, prueba y evaluación son proporcionados automáticamente por la a priori y la forma en que comparan (o penalizan) los procedimientos. Por ejemplo, si las estimaciones \(\bar{\theta}\) de \(\theta\) se comparan a través de errores al cuadradp,

\[ \mathcal{L} = ||\theta-\bar{\theta}||^2 \]

el correspondiente de Bayes es el valor esperado de \(\theta\) bajo la distribución posterior.

\[ \bar{\theta} = \dfrac{\int\theta\mathcal{l}(\theta|y)p(\theta)d\theta}{\int\mathcal{l}(\theta|y)p(\theta)d\theta} \] Cuando no se dispone de un criterio de penalización específico, el estimador suele utilizarse como el estimador por defecto, aunque existan alternativas. Por ejemppl, el estimador máximo a posterior MAP se define como, \[ \bar{\theta} = arg \max_{\theta} p(\theta | y) = arg \max_{\theta}p(\theta) \mathcal{l}(\theta|y) \]

donde la función a maximizar se suele proporcionar en forma cerrada. Sin embargo, los problema numéricos a menudo hacen que optimización para encontrar el MAP no sea trivial. Nótese también aquí la similitud de la ecuación anterior con el estimador de máxima verosimilitud (MLE): La influencia de la distribución a priori \(p(\theta)\) en la estimación desaparece progresivamente a medida que aumenta el número de observaciones n, y el estimador MAP suele recuperar las propiedades asintóticas del MLE.

\[ E[\sigma^2|y] = 1.82/(33-1) = 0.057 \] * Del mismo modo, la estimación MAP viene dada aquí por

\[ arg \max_{\theta} p(\theta | y) = 1.82/(33+1) = 0.054 \]

Distribuciones previas conjugadas

La selección de la distribución a priori es una cuestión importante en la estadística bayesiana. Cuando se dispone de información previa sobre los datos o el modelo, puede (y debe) utilizarse para construir la distribución a priori, y veremos algunas aplicaciones de esta recomendación en los siguientes capítulos. En muchas situaciones, sin embargo, la selección de la distribución a priori es bastante delicada, debido a la ausencia de información fable a priori, y en su lugar deben elegirse soluciones por defecto. Dado que la elección de la distribución a priori tiene una influencia considerable en la in! ultante, este paso inferencial debe realizarse con el máximo cuidado.

Desde un punto de vista computacional, la elección más conveniente de las distribucio estructura de la verosimilitud dentro de la prioridad. En los casos más ventajosos, los priores y los posteriors permanecen dentro de la misma familia paramétrica. Tales prioritarios se denominan conjugados. Anmque los fundamentos de este principio son demasiado avanzados para procesarlos aquí (véase, por ejemplo, Robert, 2007, Cap. 3), dichos priores existen para la mayoría de las familias habituales, incluida la distribución normal. De hecho, como se ha visto, cuando la prioridad de una media normal es normal, la posterior correspondiente también es normal.

Dado que las previas conjugadas son tales que las densidades a priori y a posteriori pertenecen a la misma familia paramétrica, el uso de las observaciones se reduce a una actualización de los parámetros de la priori. Para evitar confusiones, los parámetros que intervienen en la distribución a priori sobre el parámetro del modelo suelen llamarse hiperparámetros. (Pueden asociarse a su vez a distribuciones a priori, denominándose entonces hiperparámetros)

Para la mayoría de los propósitos prácticos, es suficiente considerar las previas conjugadas descritas en la Tabla siguiente. La derivación de cada fila es sencilla y procede de la misma aplicación de la fórmula de Bayes que para el caso normal anterior. Para las distribuciones que no están dentro de esta tabla, se puede disponer o no una previa conjugada. Una catacterística importante de las previas conjugadas es que uno tiene que seleccionar a priori dos hiperámetros, por ejemplo, una media y varianza en el caso normal. Por un lado, esto es una ventaja cuando se utiliza una prioridad conjugada, ya que sólo hay que seleccionar unos poco parámetros para determinar la distribución a priori. Por otro lado, es un inconveniente, ya que la información conocida a priori. Por otro lado, es un incoveniente,ya que la información conocida a priori sobre \(\mu\) puede ser insuficiente para determinar ambos parámetros o incompatible con la estructura impuesta por la conjugación.

Previas no informativas

No hay ninguna razón de peso para elegir las previas conjugadas como nuestras previas, excepto por su simplicidad, pero el aspecto restrictivo de las previas conjugadas puede atenuarse utilizando hiper-previas ex los propios liperparámetros. El mensaje principal es, por tanto, que es agradable trabajar con las previas conjugadas, pero que requieren una determinación de los hiperparámetros que puede resultar incómoda en álgunos entornos y que, además, puede tener un impacto duradero en la inferencia resultante,

En lugar de utilizar previas conjugados, se puede optar por una perspectiva completamente diferente y confiar en las amadas previas no informativos que tieuen como objetivo atenuar el impacto de la previa vn la inforencia resultamto. Estas previas se definen fimdamentalmente como extensiones coherentes de la distribución uniforme. Su objetivo es proporcionar wma medida de referencia que infhuya lo menos posible en la inferencia (cn relación con la información aportada por la verosimilitud). Por ejemplo, los modelos de localización

\[ y \sim p(y-\theta) \]

suelen estar asociados a previas planas \(p(\theta)=1\) (nótese que estos modelos incluyen la normal \(N(\theta,1)\) como un caso especial), mientras que los modelos de escala

\[ y \sim \dfrac{1}{\theta}f(\dfrac{y}{\theta}) \]

suelen asociarse a la transformada logarítimica de una previa plana, es decir,

\[ p(\theta)= \dfrac{1}{\theta} \]

En un entorno mgeneral, la previa (no informativa) favorecita por la maoyoría de los bayesianos es la llamada previa de Jeffreys, que está relacionada con la matriz de informacion de Fisher.

\[ \begin{aligned} J(\theta) & = E\left[\left(\dfrac{d \lg p(y|\theta) }{d\theta}\right)^2|\theta\right] = -E\left[\left(\dfrac{d^2 \lg p(y|\theta) }{d\theta^2}\right)^2|\theta\right] \\ p(\theta) & \propto [J(\theta)]^{-1/2} \end{aligned} \] \[ \begin{aligned} I^F (\theta) & = Var_\theta\left(\dfrac{\partial \log f(Y|\theta)}{\partial \theta}\right) \\ p^J(\theta) & = |I^F (\theta)|^{1/2} \end{aligned} \] Dado que la media \(\mu\) de un modelo normal es un parámetro de localización, cuando la varianza \(\sigma^2\) es conocida, la elección estándar de parámetro no informativo es una constante arbitraria \(p(\theta)\) (que se toma como t por defecto). Dado que esta previa plana corresponde formalmente al caso límite \(\tau = \infty\) en la previa normal conjugada, es fácil verificiar que esta previa no informativa está asociada a la distribución posterio \(N(x,1)\), que resulta ser la función de verosimilitud en ese caso. Una consecuencia interesante de sta observación es que el estimador MAP es tambien el estimador de máxima verosimilitud en ese caso (especial). Para el caso cuando

\[ \begin{aligned} p((\mu,\sigma^2)|y) & \propto (\sigma^{-2})^{(3+n/2)} exp \left \lbrace -\left(n(\mu-\bar{y})^2+s^2\right) \right\rbrace \\ & \propto \sigma^{-1} exp \left \lbrace -n(\mu-\bar{y})^2/2\sigma^2 \right\rbrace \times\lbrace(\sigma^2)^{-(n+2)/2}\rbrace exp \left \lbrace \dfrac{-s^2}{2\sigma^2} \right\rbrace\\ \theta & \sim N(\bar{y},\sigma^2/n)\times \mathcal{IG}(n/2,s^2/2) \end{aligned} \]

un producto de una normal condicional sobreb \(\mu\) por una gamma inversa sobre \(\sigma^2\). Por tanto, la distribución marginal posterior sobre \(\mu\) es una distribución t.

\[ \mu | y \sim \mathcal{T}(n,\bar{y},s^2/n^2) \]

\[ \int p(\theta)\mathcal{l}(\theta|y)d\theta < \infty \]

Ejercicio 10: En un estudio llevado a cabo en cierto país se encuentra que de 980 nacimientos con la condición de placenta previa, 437 eran niñas. ¿Cuánta evidencia proporcionan estos datos sobre la hipótesis de que la proporción de nacimientos femeninos en la población de placenta previa es menor que la proporción 0.458 de mujeres en la población?

\[ \begin{aligned} p\left(y|\theta \right) &= \mathcal{B}eta \left(\alpha + y,\beta + n-y\right) \\ E \left(\theta|y\right) &= \dfrac{\alpha + y}{\alpha + \beta + n}\\ Var(\theta|y) &= \dfrac{(\alpha + y)(\alpha + \beta - n)}{(\alpha\beta+n)^2(\alpha+\beta+n+1)} \end{aligned} \]

Si se usa una distribución inicial \(\mathcal{U}\left(0,1\right)\) equivalente a una \(\mathcal{B}eta \left(1,1\right)\) la distribución posterior de \(\theta\), la proporción de nacimientos niñas con la condición de placenta previa, es \[ \begin{aligned} \mathcal{B}eta \left(\theta | y+1, n-y\right) &= \mathcal{B}eta \left(\theta|437+1,980-437\right) \\ &=\mathcal{B}eta\left(\theta|438,544\right) \end{aligned} \]

Se puede calcular la media posteriori, desviación etandar, media y cuantiles del \(2.5\%\) y \(98.5\%\) de probabilidad

ls()
## [1] "diagmma"    "mmu"        "mtmu"       "n"          "normaldata"
## [6] "shift"      "stmu"       "vmu"
rm()
(media <- 438/(438+543)) 
## [1] 0.4464832
(sd <- (438*543)/((438*544)^{2})*(438+543+1))
## [1] 0.004113764
qbeta(0.5,438,543)
## [1] 0.4464468
qbeta(0.025,438,543)
## [1] 0.4155005
qbeta(0.975,438,543)
## [1] 0.4776725

Con esto, nos queda claro que el intervalo central del \(95\%\) de probabilidad a posteriori es: \(\left[0.415,0.477\right]\) Nótese que este intervalo no contiene a al proporción de niñas en la poblaci{on, con lo que se puede inferir que la proporción de placenta previa en nacimientos de niñas es menor que la proporción poblacional.

simula <- rbeta(1000,1,1)
hist(rbeta(10000,438,543))

Ejercicio 11: Ejemplo hipotético: Tasa de incidencia de malaria en el municipio de Guayaquil.

Se desea estudiar la tasa de incidencia de malaria en el municipio de Guayaquil. Se observa por ejemplo que en el año _2001: se registraron 500 casos en este municipio de 20.000 personas (250 casos por cada 10.000 persnas en promedio por año.)

Si se usa una distribución de muestreo Poisson para el y el número de casos en el municipio de 20.000 habitantes en un año, la distribución de muestreo puede extender como una Poisson \(\mathcal{P}\left(2\theta\right)\), donde \(\theta\) es la verdadera tasa de incidencia de la enfermedad a largo plazo en el municipio por cada 10.000 habitantes.

Considerando opinión de expertos, se tiene un conocimiento previo (a priori) que se espera un valor de \(\theta\) alrededor de 10, se usa una distribción previa para \(\theta=\mathcal{G}amma \left(5,0.5\right)\) la cual tiene media a 10. la distribución posterior resultante es una \(\mathcal{G}amma \left(505,2.5\right)\) con una media de 202.

Esta distribución se actualiza luego de que se dispone de observaciones adicionales. Suponiendo que la población no cambia en 5 años (20.000 habitantes), se observan 1.500 casos en 5 años. En este caso la distribucion de muestreo es una Poisson \(\mathcal{P} \left(10 \theta \right)\) y la distribución posterior es \(\mathcal{G}amma \left(1505,10.5\right)\) con media 143.3.

theta.prior <- rgamma(1000,5,scale = 2)
summary(theta.prior)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   1.604   6.849   9.483  10.000  12.644  29.424
theta.posterior <- rgamma(1000,510,scale = 2/5)
summary(theta.posterior)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   171.4   197.9   203.6   203.8   209.9   235.1
length(theta.posterior[theta.posterior > 150])/1000
## [1] 1
theta.posterior.2 <- rgamma(1000,1505,scale = 2/21)
summary(theta.posterior.2)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   134.0   141.2   143.5   143.6   145.9   154.8
length(theta.posterior.2[theta.posterior.2 > 150])/1000
## [1] 0.042
par(mfrow = c(3,1), bty = 'n')
hist(theta.prior,xlab = expression(theta),ylab = "",main = "Previa",xlim =  c(0,250))
hist(theta.posterior,xlab = expression(theta),ylab = "",main = "Distribución Posterior",xlim =  c(0,250))
hist(theta.posterior.2,xlab = expression(theta),ylab = "",main = "Distribución Posterior II",xlim =  c(0,250))

Ejercicio 12: Datos discretos En el cuadro se marca el número de accidentes mortales y fallecidos en vuelos regulares de aerolíneas por año durante un periodo de diez años. Utilizaremos los siguientes datos para ajustar un modelo de datos discreto.

library(readxl)
datos <- read_excel("tasa_aereos.xlsx")
library(knitr)
df <- datos 
kable(df, caption = "Accidentes fatales aéreos", align = c('c','c','c','c'), col.names = c("Año", "Accidentes fatales","Pasajeros fallecidos","Tasa de fallecidos"), row.names = FALSE, digits = 1, format.args = list(decimal.mark= "."))
Accidentes fatales aéreos
Año Accidentes fatales Pasajeros fallecidos Tasa de fallecidos
1976 24 734 0.2
1977 25 516 0.1
1978 31 754 0.1
1979 31 877 0.2
1980 22 814 0.1
1981 21 362 0.1
1982 26 764 0.1
1983 20 809 0.1
1984 16 223 0.0
1985 22 1066 0.1

\[ y_{i}| \theta \sim Poisson\left(\theta\right) \]

\[ Gamma \left(\alpha +10\bar{y}, \beta+10\right ) \]

\[ 1976: \dfrac{746}{0.19}\times 100=3.863 \times 10^{1} \]

\[ \theta|y \sim Gamma \left(10\bar{y},10\bar{x}\right) = Gamma\left(238,5.716 \times 10^{12}\right) \]

\[ Poisson(x\theta) = Poisson(8\times10^1) \]

\[ \begin{aligned} E\left(\theta|y\right)&= \dfrac{238}{10} = 23.8 \\ sd\left(\theta|y\right) &= \dfrac{\sqrt{238}}{10} \ = 1.54 \end{aligned} \]

\[ \begin{aligned} E\left(\tilde{y}\right|\theta) &= \theta \\ sd\left(\tilde{y}\theta \right) &= \sqrt{\theta} \end{aligned} \] - La media y la varianza de la distribución predictiva para \(\tilde{y}\) son:

\[ \begin{aligned} E\left(\tilde{y}\right|y) &= E\left(E\left(\tilde{y}|\theta, y\right)\right) = E\left(\theta|y \right)= 23.8 \\ Var\left(\tilde{y}|y\right) &= E\left(Var\left(\tilde{y}|\theta,y \right)\right) + Var\left(E\left(\tilde{y}|\theta,y \right)\right) = E\left(\theta|y \right) + Var\left(\theta|y \right) = 23.8 + (1.54)^2 = 26.2 \\ \end{aligned} \]

# Gamma y Poisson
set.seed(1)
theta <- rgamma(1000, 238)/10
y1986 <- rpois(1000, theta)
y <- sort(y1986)

print(sort(y1986)[c(25,976)])
## [1] 14 35
# aproximación normal

alfa <- 238
beta <- 10

media.theta <- alfa/beta
var.theta <- alfa/beta^{2}

media.ytilde <- media.theta
var.ytilde <- media.theta + var.theta

signif <- 0.05

cuantil <- qnorm(1-signif/2)
(LI <- media.ytilde-(cuantil * sqrt(var.ytilde)))
## [1] 13.77157
(LS <- media.ytilde+(cuantil * sqrt(var.ytilde)))
## [1] 33.82843
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
#número de pasajeros por millas en cada año 

millas <- 100000000
datos1 <- mutate(df, millas.pasa = as.numeric(format((pasa.falle/tasa.falle)*millas, digits = 4 , scientific = TRUE)))

#estimacion del numero de millas por pasajeros

suma <- apply(datos1, 2, sum)

sIMULACIÓN

set.seed(1)
theta <- rgamma(1000,238)/suma[5]
y1986 <- rpois(1000, theta*8e11)
print(sort(y1986)[c(25,976)])
## [1] 22 46
set.seed(1)
alfa <- suma[2]
beta <- suma[5]
y.tilde <- 8e11
media.theta <- alfa*beta
var.theta <- alfa / beta^{2}

media.poisson <- y.tilde * media.theta
var.poisson <- y.tilde * media.theta

media.ytilde <- media.poisson
var.ytilde <- (y.tilde*media.theta)+(y.tilde^{2}*var.theta)

signif <- 0.05
cuantil <- qnorm(1-signif/2)

(LI <- media.ytilde-(cuantil * sqrt(var.ytilde)))
##  acc.fatales 
## 1.088288e+27
(LS <- media.ytilde+(cuantil * sqrt(var.ytilde)))
##  acc.fatales 
## 1.088288e+27
set.seed(1)

theta <- rgamma(1000, suma[3])/10
y1986 <- rpois(1000, theta)
y <- sort(y1986)

print(sort(y1986)[c(25,976)])
## [1] 640 749
set.seed(1)
suma
##        anio acc.fatales  pasa.falle  tasa.falle millas.pasa 
##  1.9805e+04  2.3800e+02  6.9190e+03  1.2600e+00  5.7158e+12
theta <- rgamma(1000, suma[3])/suma[5]
y1986 <- rpois(1000, theta*8e11)
y <- sort(y1986)

print(sort(y1986)[c(25,976)])
## [1]  904 1034

Ejercicio 13 Distribución normal con media desconocida: Se extrae una muestra aleatoria de n estudiantes de una gran población y se mide sus pesos. El peso promedio de los n estudiantes de la muestra es y=150 libras. Supongamos de que los pesos la población se distribuyen normalmente con media desconocida \(\theta\) y una desviación estandar conocida de 20 libras. Suponga que su distribución previa \(\theta\) es normal con media 180 libras y una desviación estándar de 40 libras.

  1. Dé la distribución posterior para \(\theta\). Su respuesta debe estar en función de n.

\[ \begin{aligned} p\left(\theta | y_{1}, \cdots, y_{n}\right) &= p\left(\theta | \bar{y} \right) = N\left(\theta | \mu_{n}, \tau_{n}^{2}\right) \\ \mu_{n} &= \dfrac{\dfrac{1}{\tau_{0}^{2}} \mu_{0} + \dfrac{n}{\sigma^{2}}\bar{y}}{\dfrac{1}{\tau_{0}^{2}}+\dfrac{n}{\sigma^{2}}} \\ \dfrac{1}{\tau_{n}^{2}} &= \dfrac{1}{\tau_{0}^{2}}+\dfrac{n}{\sigma^{2}} \\ \theta|y &\sim N\left(\dfrac{\dfrac{1}{40^2}180+\dfrac{n}{20^2}150}{\dfrac{1}{40^2}+\dfrac{n}{20^2}},\dfrac{1}{\dfrac{1}{40^2}+\dfrac{n}{20^2}}\right) \end{aligned} \]

  1. Se toma una muestra al azar de un nuevo estudiante de la misma población y tiene un peso de \(\tilde{y}\) libras. Construya una distribución predicitva posterior para \(\tilde{y}\). Su respuesta debe estar en función de n.

\[ \begin{aligned} E\left(\tilde{y}\right|y) &= E\left(\theta| y\right) = \mu_{1} \\ Var\left(\tilde{y}|y\right) &= E\left(\sigma^2| y\right) + Var\left(\theta|y\right) + \sigma^2 + \tau_{1}^2 \\ \tilde{y}|y &\sim N\left(\dfrac{\dfrac{1}{40^2}180+\dfrac{n}{20^2}150}{\dfrac{1}{40^2}+\dfrac{n}{20^2}},\dfrac{1}{\dfrac{1}{40^2}+\dfrac{n}{20^2}}+20^2\right) \end{aligned} \]

  1. Para n = 10, construya un intervalo a psoterior al \(95\%\) para \(\theta\) y un intervalo predictivo posterio al \(\95%\) para \(\tilde{y}\).
a <- ((1/40^2)*180)+((10/20^{2})*150)
b <- (1/(40^2))+(10/(20^{2}))
a/b
## [1] 150.7317
sqrt(1/b)
## [1] 6.24695
(a/b) - (1.96*sqrt(1/b))
## [1] 138.4877
(a/b) + (1.96*sqrt(1/b))
## [1] 162.9757

\[ 150.7317 \pm 1.96 \times 6.24 = \left[138.4877 , 162.9757\right] \]

a <- ((1/40^2)*180)+((10/20^{2})*150)
b <- (1/(40^2))+(10/(20^{2}))
a/b
## [1] 150.7317
sqrt(1/b)
## [1] 6.24695
(a/b) - (1.96*sqrt((1/b)+20^{2}))
## [1] 109.664
(a/b) + (1.96*sqrt((1/b)+20^{2}))
## [1] 191.7994

\[ 150.7317 \pm 1.96 \times 20.95 = \left[109.664 , 191.7994\right] \]

  1. Haga lo mismo para n=100.
a <- ((1/40^2)*180)+((100/20^{2})*150)
b <- (1/(40^2))+(100/(20^{2}))
a/b
## [1] 150.0748
sqrt(1/b)
## [1] 1.997505
(a/b) - (1.96*sqrt((1/b)+20^{2}))
## [1] 110.6798
(a/b) + (1.96*sqrt((1/b)+20^{2}))
## [1] 189.4698

\[ 150.7317 \pm 1.96 \times 20.95 = \left[110.6798 , 189.4698\right] \] Ejercicio Modelo Multinomial:

En 1988 se hizo una encuesta pre-electoral en la elecciíon presidencial de USA. De 1447 personas encuestadas, \(y_{1}=727\) apoyaron a Bush; \(y_{2}=583\) apoyaron a Michales Dukakis y \(y_{3}=173\) apoyaron a otros candidatos.

  1. Observación: La paramertización de la distribución Gamma en R y la usada en el ejercicio difieren en el parámetro de escala.↩︎