XS3310 Teoría Estadística

.title[
# XS3310 Teoría Estadística
]
.subtitle[
## I Semestre 2023
]
.author[
### Escuela de Estadística
]
.date[
### 16-06-2023
]

---

# ¿Qué hemos visto hasta ahora?

Introducción a la Estadística Bayesiana: filosofía, historia y un poco de cálculo. Distribuciones previas (a priori).

# ¿Qué vamos a discutir hoy?

Más distribuciones previas (informativa, no informativa).

---

## Estadística Bayesiana

El **teorema de Bayes** establece que para dos eventos `$A$` y `$B$` se cumple

`$$P(A|B) = \frac{P(B|A)P(A)}{P(B)}$$`

**Estadística Bayesiana** considera `$\theta$` como un evento aleatorio y el evento `$X$` como la información empírica brindada por los datos:

`$$\pi(\theta|x) = \frac{\mathcal{L}(\theta|x)\pi(\theta)}{f(x)}$$`
**Teorema**. Bajo las condiciones anteriores:

$$
\pi(\theta|X) =
\dfrac{f(X_1|\theta)\cdots f(X_n|\theta)\pi(\theta)}{g_n(X)} 
$$

para `$\theta \in \Omega$`, donde `$g_n$` es una constante de normalización.

---

## Selección de distribuciones previas
	
### Previas informativas
	
Una forma de seleccionar una previa es utilizando una previa informativa. Claramente, para poder hacer uso de esta debemos tener bastante conocimiento sobre el fenómeno de interés que estamos estudiando (como en el ejemplo que hicimos al inicio sobre la proporción de desempleo en el país). Para ello es imperativo tener información de varios expertos para poder llegar a construir el modelo estadístico que mejor represente esa información. Es importante destacar algo del uso de previas muy informativas; entre más informativa sea la previa que estamos utilizando entonces más datos se necesitan para tratar de observar algo nuevo. Imaginen que están casi `$100\%$` seguros sobre cuánto debería ser el valor de cierto parámetro desconocido. Para que ustedes piensen cambien de opinión entonces van a requerir de una gran cantidad de evidencia que apunte a lo contrario. Eso es lo que pasa si se usa una previa muy informativa; la posterior se va a parecer mucho a esta y solo podrá cambiar si existen muchos datos que apuntan a lo contrario.

---
## Familias conjugadas

**Definición**. Sea `$X_1,\dots, X_n$` i.i.d. condicional dado `$\theta$` con
densidad `$f(X|\theta)$`. Sea `$\psi$` la familia de posibles densidades previas
sobre `$\Omega$`. Si, sin importar los datos, la posterior pertenece a `$\psi$`,
entonces decimos que `$\psi$` es una familia conjugada de previas.

**Ejemplos**:

* La familia Beta es familia conjugada para muestras según una Bernoulli.
*  La familia Gama es familia conjugada para muestras exponenciales.
* Para el caso Poisson, si `$X_1,\dots,X_n\sim Poi(\lambda)$`,entonces la familia
  Gamma es familia conjugada.

La función de densidad de una Poisson es `$P(X_i = k) = e^{-\lambda}\dfrac{\lambda^k}{k!}$`. La verosimilitud corresponde a 
`$$f_n(X|\lambda) = \prod_{i=1}^{n}e^{-\lambda}\dfrac{\lambda^{X_i}}{X_i!} = \dfrac{e^{-n\lambda}\lambda^Y}{\prod_{i=1}^n X_i!}~~~ \text{ donde } Y=\sum_{i=1}^n X_i.$$`
La previa de `$\lambda$` está definida por `$\pi(\lambda)\propto\lambda^{\alpha-1}e^{-\beta\lambda}$`. Por lo tanto, la posterior es
$$ \pi(\lambda|X) \propto \lambda^{Y+\alpha-1}e^{-(\beta+n)\lambda} \implies
 \lambda|X \sim \Gamma(Y+\alpha,\beta+n)$$

---
* En el caso normal, si `$X_1,\dots,X_n\sim N(\theta,\sigma^2)$`,entonces la familia normal es conjugada si `$\sigma^2$` es conocido.

Si `$\theta \sim N(\mu_0,V_0^2) \implies \theta|X \sim N(\mu_1, V_1^2)$` donde,
`$$\mu_1 = \dfrac{\sigma^2\mu_0 + nV_0^2 \bar X_n}{\sigma^2 + nV_0^2}  = \dfrac{\sigma^2}{\sigma^2 + nV_0^2}\mu_0 + \dfrac{nV_0^2}{\sigma^2 + nV_0^2}\bar X_n$$`

Combina de manera ponderada la previa y la de los datos.

**Ejemplo**

Considere una verosimilitud Poisson `$(\lambda)$` y una previa

`$$\pi(\lambda) = \begin{cases}2e^{-2\lambda} & \lambda> 0 \\ 0 & \lambda \leq 0\end{cases} \quad \lambda \sim \Gamma(1,2)$$`

Supongamos que es una muestra aleatoria de tamaño `$n$`. ¿Cuál es el número de observciones para reducir la varianza, a lo sumo, a 0.01?

Por teorema de Bayes, la posterior `$\lambda|x \sim \Gamma(y+1,n+2)$`. Luego, la varianza de la Gamma es

`$$\dfrac{\alpha}{\beta^2} = \dfrac{\sum x_i + 1}{(n+2)^2} \leq 0.01 \implies \dfrac{1}{(n+2)^2} \leq \dfrac{\sum x_i + 1}{(n+2)^2} \leq 0.01 \implies 100 \leq (n+2)^2 \implies n\geq 8$$`

---
**Teorema**. Si `$X_1,\dots,X_n \sim N(\theta, \sigma^2)$` con `$\sigma^2$` conocido y la previa es `$\theta \sim N(\mu_0,V_0^2)$`, entonces `$\theta|X\sim N(\mu_1,V_1^2)$` donde
$$ \mu_1 =  \dfrac{\sigma^2\mu_0 + nV_0^2 \bar X_n}{\sigma^2 + nV_0^2}, \quad V_1^2 = \dfrac{\sigma^2V_0^2}{\sigma^2 + nV_0^2}$$
*Prueba*:

* **Verosimilitud**:

`$$f_n(X|\theta) \propto \exp\left[- \dfrac{1}{2\sigma^2} \sum_{i=1}^{n}(X_i-\theta)^2\right]$$`
Luego, 
`\begin{align*}
\sum_{i=1}^n (X_i-\theta)^2 & = \sum_{i=1}^n (X_i-\bar X + \bar X - \theta)^2 \\
& = n(\bar X - \theta)^2 + \sum_{i=1}^n (X_i-\bar X)^2 + \underbrace{2 \sum_{i=1}^n (X_i-\bar X)(\bar X - \theta)}_{= 0 \text{ pues } \sum Xi = n\bar X)}
\end{align*}`

Entonces
$$ f_n(X|\theta) \propto \exp\left[-\dfrac{n}{2\sigma ^2}(\bar X - \theta )^2\right].$$

---
* **Previa**:

$$ \pi(\theta) \propto \exp\left[-\dfrac{1}{2V_0^2}(\theta - \mu_0)^2\right].$$

* **Posterior**:

$$ \pi(\theta|X) \propto \exp\left[-\dfrac{n}{2\sigma ^2}(\bar X - \theta )^2-\dfrac{1}{2V_0^2}(\theta - \mu_0)^2\right].$$

Con `$\mu_1$` y `$V_1^2$` definidos anteriormente, se puede comprobar la siguiente identidad:

`$$-\dfrac{n}{\sigma ^2}(\bar X - \theta )^2-\dfrac{1}{V_0^2}(\theta - \mu_0)^2= \dfrac{1}{V_1^2}(\theta-\mu_1)^2 + \underbrace{\dfrac{n}{\sigma^2 + nV_0^2}(\bar X_n- \mu_0)^2}_{\text{Constante con respecto a }\theta}$$`
Por lo tanto, `$$\pi(\theta|X) \propto \exp\left[-\dfrac{n}{2V_1^2}(\theta -\mu_1)^2\right]$$`

---
*Media posterior*:

`$$\mu_1 = \underbrace{\dfrac{\sigma^2}{\sigma^2 + nV_0^2}}_{W_1}\mu_0 + \underbrace{\dfrac{nV_0^2}{\sigma^2 + nV_0^2}}_{W_2}
\bar X_n$$`

**Afirmaciones**:

1) Si `$V_0^2$` y `$\sigma^2$` son fijos, entonces `$W_1 \xrightarrow[n\to \infty]{}0$` (la importancia de la media empírica crece conforme aumenta `$n$`).

2) Si `$V_0^2$` y `$n$` son fijos, entonces `$W_2 \xrightarrow[\sigma^2\to \infty]{}0$` (la importancia de la media empírica decrece conforme la muestra es menos precisa).

3) Si `$\sigma^2$` y `$n$` son fijos, entonces `$W_2 \xrightarrow[V_0^2\to \infty]{}1$` (la importancia de la media empírica crece conforma la previa es menos precisa).

---

**Ejemplo (determinación de n)**

Sean `$X_1,\dots, X_n \sim N(\theta,1)$` y `$\theta\sim N(\mu_0,4)$`. Sabemos que $$V_1^2 = \dfrac{\sigma^2V_0^2}{\sigma^2 + nV_0^2}. $$
Buscamos que `$V_1\leq 0.01$`, entonces
$$ \dfrac{4}{4n+1}\leq 0.01 \implies n\geq 99.75 \text{ (al menos 100 observaciones)}$$

---

## Selección de distribuciones previas

### Previas conjugadas

Un tipo de previas informativas que son sumamente convenientes de utilizar son las **previas conjugadas**. Se dice que una previa es conjugada si la distribución de la posteriori pertenece a la misma familia que la distribución de la previa. El ejemplo anterior donde utilizamos una previa Gamma y terminamos con una posteriori Gamma es un ejemplo de una previa conjugada. Este tipo de previas son muy convenientes pues nos aseguran que vamos a tener una distribución a posteriori conocida, lo único que cambiaría son los parámetros de la distribución. Sin embargo, esto no quiere decir que la Gamma siempre vaya a ser una previa conjugada; esto solo va a pasar si los datos son Poisson. Por lo tanto, una parte importante que permite que la previa sea o no sea conjugada es la distribución de la población.

---

## Distribuciones conjugadas conocidas

Distribución             | previa conjugada                  |
-------------------------|-----------------------------------|
`$$\text{Bernoulli}(p)$$`    | `$$p \sim Beta(\alpha,\beta)$$`  |
`$$\text{Binomial}(n,p)$$`   | `$$p \sim Beta(\alpha,\beta)$$`  |
`$$\text{Binomial Negativa}(n,p)$$`| `$$p \sim Beta(\alpha,\beta)$$`|
`$$\text{Poisson}(\lambda)$$`      | `$$\lambda \sim Gamma(\alpha,\beta)$$`|
`$$\text{Exponencial}(\theta)$$`| `$$\theta \sim Gamma(\alpha,\beta)$$`|
`$$\text{Normal}(\mu,\sigma^2)$$`| `$$\mu \sim N(\mu_0,\sigma^{2}_{0}) \quad \text{ y } \quad \sigma^{2} \sim GammaInversa(\alpha,\beta)$$`|

---

## Distribuciones conjugadas conocidas

https://en.wikipedia.org/wiki/Conjugate_prior

https://www.johndcook.com/blog/conjugate_prior_diagram/

---

## Previas no informativas
	
Este tipo de previa es posiblemente el más utilizado en la práctica. Con modelos relativamente simples utilizar una previa no informativa por lo general brinda resultados muy similares a los resultados frecuentistas, mientras que con modelos jerárquicos más complejos los resultados sí pueden ser más distintos. No obstante, la mayoría del tiempo en que se quiere hacer inferencia sobre parámetros desconocidos no se tiene mucha información al respecto, aparte de un posible rango en donde se puedan encontrar; por esto es más atractivo utilizar una previa no informativa. Diremos que una previa es no informativa si le da la libertad a los datos de encontrar los mejores valores de los parámetros. 
	
Por mucho tiempo se tuvo la idea de que las previas no informativas eran exclusivamente las uniformes, pues asignaban probabilidades iguales a cualquier parámetro. Sin embargo, hay quienes criticaron esto, como Fisher, diciendo que no es posible que la uniforme sea siempre no informativa.

---

## Previas no informativas

El argumento de Fisher era algo así:
	
>	Supongamos que tenemos un parámetro desconocido `$\theta$` que representa la probabilidad de éxito de un experimento Bernoulli. Supongamos que no sabemos nada de `$\theta$` entonces decimos que `$\theta \sim Unif(0,1)$`. Ahora, si no sabemos nada de `$\theta$` entonces tampoco sabemos nada de `$\lambda = -\ln(\theta)$`, por lo que también podríamos preferir una previa uniforme para `$\lambda$`. Sería lógico pensar que, mediante las técnicas de transformaciones, la transformación aplicada a la previa de `$\theta$` llegue a la previa de `$\lambda$`. Sin embargo hay un problema lógico pues si aplicamos las técnicas de transformación a la previa de `$\theta$` no vamos a llegar a la misma previa de `$\theta$`. Esto sugiere que una distribución uniforme no es un buen ejemplo de una previa no informativa.

---

## Previas no informativas
	
Un tipo de previa no informativa es la **previa de Jeffreys**. Estas hacen uso de la Información de Fisher, previamente utilizada en el Tema 1. Por lo tanto, si tenemos una muestra aleatoria `$X_1 , X_2 , ... , X_n$` de una población con función de densidad `$f_{X}(x|\theta)$` entonces la información de Fisher se define como:
	
`$$I(\theta) = -E\left[ \frac{\partial^{2} \ln(f_{X}(x|\theta)) }{\partial \theta^{2}} \right]$$`
	
Por lo tanto, la previa de Jeffreys se define como:
	
`$$\pi_{J}(\theta) = c\sqrt{I(\theta)}$$`
	
Donde `$c$` es una constante positiva mayor a cero que asegura que esta función de densidad integra a uno. Como `$c$` es constante entonces podemos decir que `$\pi_{J}(\theta) \propto \sqrt{I(\theta)}$`, por lo que muchas veces no nos importa el valor de `$c$` para encontrar la distribución a posteriori a partir de la previa de Jeffreys.

---

## Previas no informativas
	
Ejemplo: Sea `$X_1 , X_2 , ... , X_n$` una muestra aleatoria tal que `$X_j \sim N(\mu,1)$`. Encuentre la previa de Jeffreys para `$\mu$`. 
	
Solución: Recordemos la función de densidad de una Normal:
	
`$$f_{X}(x|\mu) = \sqrt{2\pi} e^{-\frac{(x-\mu)^{2}}{2}}$$`
	
Por lo que el logaritmo natural de esta sería:
	
`$$\ln(f_{X}(x|\mu))  = \frac{1}{2}\ln(2\pi) - \frac{(x-\mu)^{2}}{2}$$`
	
Derivamos dos veces con respecto a `$\mu$`:
	
`$$\frac{\partial \ln(f_{X}(x|\mu)) }{\partial \mu} = x - \mu$$`
	
`$$\Rightarrow \frac{\partial^{2} \ln(f_{X}(x|\mu)) }{\partial \mu^{2}} = -1$$`
	
	
	
---

## Previas no informativas

Finalmente obtenemos
	
`$$I(\mu) = 1$$`
	
Por lo tanto, podemos concluir que `$\pi_{J}(\theta) \propto 1$`, es decir, es proporcional a una constante. Por lo tanto, la previa de Jeffreys para estimar a `$\mu$` sería una distribución Uniforme, escogida en un rango bastante amplio.

Hay un par de puntos importantes de destacar cuando se usa la previa de Jeffreys. El primero de ellos es que no siempre se va a llegar a una función de densidad propia (es decir, una función de densidad que integre a 1 en su dominio). Por lo general se ignora este problema si la posteriori si es propia. Por lo tanto, siempre que se vaya a utilizar una previa impropia hay que revisar que la posteriori sea propia, si no los resultados no tendrían sentido. El segundo punto es un poco más filosófico y trata con el hecho de que las previas se deben elegir antes de ver los datos. La previa de Jeffreys usa la distribución de los datos para encontrar una previa lo que contradice lo que muchos dicen sobre el planteamiento de la previa.

---

## Densidades previas impropias

**Definición**. Sea `$\pi$` una función positiva cuyo dominio está en `$\Omega$`. Suponga que `$\int\pi(\theta)\;d\theta = \infty$`. Entonces decimos que `$\pi$` es una **densidad impropia**.

**Ejemplo**: `$\theta \sim \text{Unif}(\mathbb{R})$`, `$\lambda \sim \text{Unif}(0,\infty)$`.

Una técnica para seleccionar distribuciones impropia es sustituir los hiperparámetros previos por 0.

---

## Densidades previas impropias

**Ejemplo**:

Se presenta el número de soldados prusianos muertos por una patada de caballo (280 conteos, unidades de combate en 20 años).

| Unidades | Ocurrencias |  
|----------|-------------|
| 144      | 0           |   
|    91    | 1           |  
|    32    | 2           |   
|       11 | 3           |  
|         2| 4           |

* Muestra de Poisson: `$X_1 = 0, X_2 = 1, X_3 = 1,\dots, X_{280} = 0 \sim \text{Poi}(\lambda)$`.

* Previa: `$\lambda \sim \Gamma(\alpha, \beta)$`.

* Posterior: `$\lambda|X \sim \Gamma(y+\alpha, n+\beta) = \Gamma(196 + \alpha, 280 + \beta)$`.

---

Sustituyendo, `$\alpha=\beta = 0$`
`\begin{align*}
\pi(\lambda) &= \dfrac{1}{\Gamma(\alpha)}\beta^\alpha\lambda^{\alpha-1}e^{\beta\lambda}  \\
& \propto \lambda^{\alpha-1}e^{-\lambda\beta} \\
&=\dfrac{1}{\lambda}
\end{align*}`
donde `$\displaystyle\int_{0}^{\infty}\dfrac{1}{\lambda} d\lambda = \infty$`.

Por teorema de Bayes, `$$\theta|X \sim \Gamma(196,280)$$`
	
---

```r
alfa1 <- 196 
beta1 <- 280
theta <- seq(0,2,by=0.01)
posterior <- dgamma(theta,shape = alfa1, rate = beta1)
plot(theta,posterior, type="l")
abline(v=196/280) # media a posteriori = media muestral
```

---

## ¿Qué discutimos hoy?

Distribuciones previas (a priori).

## ¿Qué nos falta para terminar el curso?

Estadística Bayesiana: inferencia (estimación puntual, intervalos de credibilidad y factor de Bayes).