10-intro-bayesiana.Rmd

# Introducción a inferencia bayesiana

```{r setup, include=FALSE, message=FALSE}
library(tidyverse)
library(patchwork)
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE, fig.align = 'center', fig.width = 5, fig.height=3)
comma <- function(x) format(x, digits = 2, big.mark = ",")
theme_set(theme_minimal())
```

Consideremos la siguiente situación que se utiliza en análisis de
forense [@Lindley1977]. Cierto material se encuentra en la escena del
crimen y en la ropa de un sospechoso. La pregunta clave es si este
material tiene las mismas características en ambas muestras. Pues, si
así lo fuera habría evidencia que sugiere que el sospechoso es quién
incurrió en el crimen. Por ejemplo, si hay vidrios rotos en la escena
del crimen y se encuentran vidrios rotos en la ropa del sospechoso se
puede hacer un análisis forense para determinar si los vidrios
provienen del mismo lugar. Para esto se estudia el índice de
refracción de los vidrios.

```{block, type = "ejercicio"}

Dado lo que hemos visto hasta el momento, podemos pensar rápidamente en 
hacer inferencia estadística y podemos proponer un contraste de hipótesis. 
El más sencillo que podemos proponer es una prueba de diferencia en medias
para muestras de dos subpoblaciones. Sin embargo, ¿cuáles crees que serían 
las limitantes de este acercamiento? 

*Pista:* ¿qué pasa si los vidrios de la escena del crimen son vidrios 
comúnes y corrientes? ¿Qué pasa si los vidrios de la escena del crimen
son vidrios polarizados?

```

En esta sección veremos un marco para realizar inferencia estadística
que nos ayuda resolver situaciones como la planteada en el caso
anterior.

## Introducción {-}

Para esta sección seguiremos principalmente @Kruschke. Adicionalmente
puedes ver la sección correspondiente de @Chihara.

En las secciones anteriores estudiamos el método de máxima
verosimilitud y métodos de remuestreo. Esto lo hemos hecho para
estimar parámetros, y cuantificar la incertidumbre qué tenemos acerca
de valores poblacionales.  Independiente del método utilizado, la
incertidumbre que hemos aprendido a incorporar en nuestros estimadores
ha seguido el espirítu frecuentista. Es decir, nuestra incertidumbre
la hemos incorporado bajo distinats realizaciones hipotéticas de
nuestro conjunto de datos. Estás realizaciones las hemos simulado
considerando el proceso generador de datos, o aproximaciones razonables a éste.  

La inferencia bayesiana tiene por objetivo incorporar incertidumbre en
nuestras estimaciones. Sin embargo, la incertidumbre en el contexto
bayesiano se deriva de principios diferentes.

En particular tenemos las siguientes propiedades en el contexto bayesiano:  

- Igual que en máxima verosimilitud, la inferencia bayesiana comienza
con modelos probabilísticos y observaciones.
- En contraste con máxima verosimilitud, la inferencia bayesiana está
diseñada para incorporar información previa, o de expertos, que
tengamos acerca de los parámetros de interés.
- La inferencia bayesiana cubre como caso particular métodos basados
en máxima verosimilitud.

El concepto probabilístico básico que utilizamos para construir estos
modelos y la inferencia es el de probabilidad condicional: la
probabilidad de que ocurran ciertos eventos dada la información
disponible del fenómeno que nos interesa. Lo cual nos lleva a pensar
en la definición de probabilidad en sí.

```{block, type = "ejercicio"}

¿Cómo defines lo que es probabilidad? 

```

Posiblemente para contestar la pregunta anterior pensaste en lo que
hemos hecho durante el curso. Es decir, consideraste construir un
histograma donde la altura de las barras son proporcionales a la
frecuencia de cada valor. Y has pensado que la probabilidad de la
realización de un evento incierto es la frecuencia relativa de la
realización de dicho suceso. Esto hace referencia a la construcción
**frecuentista** del concepto de probabilidad.

Sin embargo, ¿cómo puedes definir la probabilidad de un evento que
nunca has observado? Es decir, para un evento el cual no hay registro
anterior o no hay evidencia histórica. La construcción de probabilidad
bajo el enfoque **bayesiano** considera que hay un agente
interactuando con un medio incierto. Bajo este enfoque el agente toma
decisiones de manera *coherente* entre distintas alternativas. Si el
evento al que se enfrenta el tomador de decisión es incierto entonces
el agente reflejará en la función de probabilidad sus *creencias*
sobre la realización de dicho evento. Es decir, bajo el enfoque
bayesiano la probabilidad se construye como el grado de creencia
asociado a la realización de un evento incierto. Lo anterior lo podemos 
escribir de manera formal, como sigue. 

```{block, type = "mathblock"}

**Definición.** (Elementos de un problema de decisión) Un problema 
de decisión bajo un ambiente de incertidumbre se define por la 
cuarteta $(\mathcal{D}, \mathcal{E}, \mathcal{C}, \leq)$, donde: 

- $\mathcal{D}$: Espacio de alternativas. 
- $\mathcal{E}$: Espacio de eventos inciertos. 
- $\mathcal{C}$: Espacio de consequencias. 
- $\leq$: Relación de preferencias$^*$. 

```

```{block2, type = "mathblock"}

**Definición.** (Cuantificación de los eventos inciertos y consecuencias)

- La información que el tomador de decisiones tiene sobre la posible 
ocurrencia de los eventos inciertos puede ser cuantificada a través 
de la *función de probabilidad* sobre el espacio $\mathcal{E}$. 

- Las preferencias sobre el espacio de consecuencias se 
pueden ordenar a través de una *función de utilidad* sobre el 
espacio $\mathcal{C}$ de tal forma que 
$$ c_{i,j} \leq c_{i',j'} \quad \text{si y sólo si} \quad u(c_{i,j}) \leq u(c_{i',j'})\,.$$ 

```

```{block2, type = "mathblock"}

**Definición.** (Probabilidad). Sea $E \subset \mathcal{E}$ un evento incierto *relevante*, se define la probabilidad de $E$ como 
$$P(E) = \textsf{Area}(R)\,,$$
donde $R$ cumple que $d_R = d_E$. 

De esta manera la función de probabilidad satisface:  

- $0 \leq P(E) \leq 1.$
- $P(\emptyset) = 0.$
- $P(\mathcal{E}) = 1.$
- Si $E \cap F = \emptyset$ entonces $P(E \cup F) = P(E) + P(F)$.

```

## Un primer ejemplo completo de inferencia bayesiana {-}

Consideremos el siguiente problema: Nos dan una moneda, y sólo sabemos
que la moneda puede ser una de dos posibilidades. Una opción es que la
moneda tiene probabilidad $3/5$ de tirar sol (está cargada a sol). O
alternativamente, puede ser una moneda cargada a águila, con
probabilidad $2/5$ de tirar sol. Es decir, tenemos información previa
al lanzamiento de la moneda que nos dice qué posibilidades tenemos
para la moneda que usaremos para nuestro experimento.

Vamos a lanzar la moneda dos veces y observamos su resultado (águila o
sol).  Queremos decir algo acerca de qué tan probable es que hayamos
tirado la moneda cargada a sol o la moneda cargada a águila. El proceso 
generador de nuestras observaciones (verosimilitud) lo podemos escribir como 

$$ P(X = x \,|\, \theta) = {2 \choose x} \theta^x (1- \theta)^{2-x}\,.$$ 

En este caso, tenemos dos variables: $X$, que cuenta el número de
soles obtenidos en el experimento aleatorio, y $\theta$, que da la
probabilidad de que un volado resulte en sol (por ejemplo, si la
moneda es justa entonces $\theta = 0.5$). 

¿Qué cantidades podríamos usar para evaluar qué moneda es la que
estamos usando? Si hacemos el experimento, y tiramos la moneda 2
veces, podríamos considerar la probabilidad
$$P(\theta = 0.4 \,|\, X = x)\,,$$
donde $x$ es el número de soles que obtuvimos en el experimento. Esta
es la probabilidad condicional de que estemos tirando la moneda con
probabilidad de sol 2/5 dado que observamos $x$ soles. Por ejemplo, si
tiramos 2 soles, deberíamos calcular
$$P(\theta=0.4 \,|\, X=2)\,.$$

¿Cómo calculamos esta probabilidad? ¿Qué sentido tiene?

Usando reglas de probabilidad condicional, podríamos calcular
$$P(\theta=0.4|X=2) = \frac{P(X=2 | \theta = 0.4) P(\theta =0.4)}{P(X=2)}\,.$$

Nota que en el numerador uno de los factores, $P(X=2 | \theta = 0.4),$ es la
verosimilitud. Así que primero necesitamos evaluarla:
$$P(X=2|\theta = 0.4) = (0.4)^2 = 0.16\,.$$

La novedad es que ahora tenemos que considerar la probabilidad
$P(\theta = 0.4)$. Esta cantidad no la habíamos encontrado
antes. Tenemos que pensar que bajo nuestro conocimiento previo, el
tipo de moneda que recibimos para nuestro juego es aleatoria y en este
caso sólo puede tomar dos valores.  

```{block, type = "comentario"}

En diversos textos esto se *simplifica* y se escribe que la moneda
con la que hacemos los lanzamientos es una variable aleatoria. ¡Pero esto es un
abuso de lenguaje! Una vez que recibimos la moneda, ésta la utilizamos para
generar nuestros datos. No hay ningún mecanismo aleatorio que nos cambie la
moneda entre lanzamientos. Entonces con $P(\theta)$ lo que estamos codificando
son nuestros grados de creencia de haber recibido una moneda cargada a sol o
cargada a águila.

```

Supongamos que en este caso, nos dicen que la moneda se escoge al
azar de una bolsa donde hay una proporción similar de los dos tipos de moneda
(0.4 ó 0.6). Es decir el espacio parametral es $\Theta = \{0.4, 0.6\},$  y las
probabilidades asociadas a cada posibilidad son las mismas. Es decir, tenemos
\begin{align}
P(\text{tener la moneda cargada a águila}) &= P(\theta = 0.4) \\
&= P(\theta = 0.6)\\
&= P(\text{tener la moneda cargada a sol}) \\
&= 0.5\,,
\end{align}
que representa la probabilidad de escoger de manera aleatoria la moneda con
una carga en particular.

Ahora, lo que nos falta es calcular $P(X=2)$, pero con el trabajo que
hicimos esto es fácil.  Pues requiere usar reglas de probabilidad
usuales para hacerlo. Podemos utilizar **probabilidad total** y después **probabilidad condicional**
\begin{align}
P(X) &= \sum_{\theta \in \Theta} P(X, \theta)\\
&= \sum_{\theta \in \Theta} P(X\, |\, \theta) P(\theta),
\end{align}
lo cual en nuestro ejemplo se traduce en escribir

$$ P(X=2) = P(X=2|\theta = 0.4)P(\theta = 0.4) + P(X=2|\theta=0.6)P(\theta =0.6)\,,$$
por lo que obtenemos
$$P(X=2) = 0.16(0.5) + 0.36(0.5) = 0.26\,.$$

```{block, type = "comentario"}

La función de probabilidad $P(X)$ tiene un lugar muy especial en inferencia 
bayesiana y dependiendo el uso recibe dos nombres. 

- Si la evaluamos en nuestros datos $P(X = 2)$ se conoce como la **evidencia**
de nuestros datos. 
- Si la evaluamos en un caso hipótetico $P(X = x^*)$ se conoce como la **predictiva previa**.

```

Finalmente la probabilidad de haber escogido la moneda con carga $2/5$ dado que
observamos dos soles en el lanzamiento es

$$P(\theta=0.4|X=2) = \frac{0.16(0.5)}{0.26} \approx  0.31\,.$$

Esto se conoce como la **probabilidad posterior** de que estemos
tirando la moneda con carga a águila.  Nota que nuestro conocimiento
inicial baja de 0.5 (nuestra información inicial) a 0.31. Es decir, 
los lanzamientos que realizamos aumentan nuestra creencia de haber recibido
una moneda cargada a soles.

Este es un ejemplo completo, aunque muy simple, de inferencia
bayesiana. La estrategia de **inferencia bayesiana implica tomar
decisiones basadas en probabilidades posteriores.** Y esto es por que
bajo el marco de teoría de decisión en situaciones con incertidumbre
tenemos el siguiente teorema.

```{block, type = "mathblock"}

**Teorema.** (Criterio de decisión Bayesiano)
Considérese el problema de decisión definido por $D = \{d_1, \ldots, d_k\}$
donde 
$$ d_i = \{ c_{i,j} | E_j, j = 1, \ldots, m_i\}\qquad i = 1, \ldots, k\,.$$
Sea $P(E_{i,j}| d_i)$ la probabilidad de que suceda $E_{i,j}$ si se elige
la opción $d_i$ y sea $u(c_{i,j})$ la utilidad de la consecuencia $c_{i,j}$. 
Entonces, la cuantificación de la opción $d_i$ es su utilidad esperada
$$\mathbb{E}[u(d_i)] = \sum_{j=1}^{m_i} u(c_{i,j}) P(E_{i,j}| d_i)\,.$$
La decisión óptima es aquella que 
$$ d^* = \arg \max_i \mathbb{E}[u(d_i)]\,.$$

```

Nota que la derivación de la probabilidad posterior arriba (que
argumentamos fue obtenida con reglas de probabilidad condicional) se
extiende a la regla de Bayes o teorema de Bayes en situaciones mas
generales.


```{block, type = "mathblock"}

**Teorema.** (Teorema de Bayes) Sea $\{E_j: j \in J\}$ una partición finita de 
$\mathcal{E}$. Es decir, que $E_i\cap E_j = \emptyset$ y $\cup_{j \in J} E_j = \mathcal{E}$. Sea $Z \neq \emptyset$, entonces
\begin{align}
P(E_i |Z) = \frac{P(Z|E_i)P(E_i)}{\sum_{j \in J}P(Z|E_j)P(E_j)}\,.
\end{align}


```


```{block, type='ejercicio'}
¿Cuál sería la estimación de máxima verosimilitud para este problema? ¿Cómo
cuantificaríamos la incertidumbre en la estimación de máxima verosimilitud?
```

Finalmente, podríamos hacer predicciones para lanzamientos nuevos
condicionando en nuestra información actualizada. Si ${X}_{nv}$ es una
nueva tirada de la moneda que estamos usando, nos
interesaría saber:
$$P({X}_{nv}=\mathsf{sol}\, | \, X=2)\,.$$

Notemos que un volado adicional es un resultado binario. Como tenemos las
probabilidades posteriores $P(\theta|X=2)$, podemos usar probabilidad total
condicionado en $X=2$:
\begin{align*}
P({X}_{nv}=\mathsf{sol}\, | \, X=2) & = \sum_{\theta \in \Theta} P({X}_{nv}=\mathsf{sol}, \theta \, | \, X=2) & \text{(probabilidad total)}\\
&= \sum_{\theta \in \Theta} P({X}_{nv}=\mathsf{sol}\, | \theta , X=2) P(\theta \, | \, X=2) & \text{(probabilidad condicional)}\\
&= \sum_{\theta \in \Theta} P({X}_{nv}=\mathsf{sol}\, | \theta ) P(\theta \, | \, X=2)\,, & \text{(independencia condicional)}
\end{align*}

lo que nos da el siguiente cálculo
$$P(X_{nv}=\mathsf{sol}\, |\, \theta=0.4) \,  P(\theta=0.4|X=2) \,  +\, P(X_{nv}=\mathsf{sol}|\theta = 0.6) \, P(\theta =0.6|X=2)\,.$$
Es decir, promediamos ponderando con las probabilidades posteriores.
Por lo tanto obtenemos
$$P(X_{nv} = \mathsf{sol}|X=2) =  0.4 ( 0.31) + 0.6 (0.69) = 0.538\,.$$

```{block, type = "comentario"}

La función de probabilidad $P(X_{nv} | X)$ recibe un nombre especial en 
inferencia bayesiana. Se conoce como la función de probabilidad *predictiva posterior*. 


```

```{block, type = "mathblock"}

**Definición.** La **función de probabilidad predictiva** es la función de probabilidad 
marginal que determina qué valores de X son los mas probables. 

Lo que se conoce sobre nuestras mediciones $X$ es que condicionado en $\theta$ resulta en el proceso
generador de datos (verosimilitud). Cuando no hemos visto realizaciones de este 
proceso. Es decir, que no tenemos información $X \sim P(X|\theta)$ entonces 
nuestro estado de conocimiento marginal, $P(X)$, se conoce como la 
**distribución predictiva previa** y se obtiene 
$$ P(X) = \sum_{\theta \in \Theta}^{} P(X | \theta) P(\theta)\,.$$

Cuando hemos recabado información, utilizamos probabilidad condicional para
determinar qué valores nuevos $X^*$ son los mas probables una vez que hemos observado mediciones. Esto es $P(X^*|X)$ y se obtiene al marginalizar utilizando 
nuestro estado actual de conocimiento. Es decir, 
$$ P(X^*|X) = \sum_{\theta \in \Theta}^{} P(X^* | \theta) P(\theta|X)\,.$$
Esta última es la **distribución predictiva posterior**.

```

```{block, type = "comentario"}

Nota que para definir las función de probabilidad predictiva previa (posterior)
marginalizamos con respecto a nuestro estado de conocimiento previo (posterior)
que denotamos por $P(\theta)$ (o respectivamente $P(\theta | \textsf{datos})$).
En este sentido decimos, que nuestra incertidumbre sobre el valor del parámetro
$\theta$ ha sido incorporada para expresar nuestras predicciones
(probabilisticas) de una realización hipotética.

```

#### Ejercicios {-}

Cambia distintos parámetros del número de soles observados, las probabilidades
de sol de las monedas, y las probabilidades iniciales de selección de las
monedas.

```{r}
n_volados <- 2
# posible valores del parámetro desconocido
theta = c(0.4, 0.6)
# probabilidades iniciales
probs_inicial <- tibble(moneda = c(1, 2),
                        theta = theta,
                        prob_inicial = c(0.5, 0.5))
probs_inicial
# verosimilitud
crear_verosim <- function(no_soles){
    verosim <- function(theta){
      # prob de observar no_soles en 2 volados con probabilidad de sol theta
      dbinom(no_soles, 2, theta)
    }
    verosim
}
# evaluar verosimilitud
verosim <- crear_verosim(2)
# ahora usamos regla de bayes para hacer tabla de probabilidades
tabla_inferencia <- probs_inicial %>%
  mutate(verosimilitud = map_dbl(theta, verosim)) %>%
  mutate(inicial_x_verosim = prob_inicial * verosimilitud) %>%
  # normalizar
  mutate(prob_posterior = inicial_x_verosim / sum(inicial_x_verosim))

tabla_inferencia %>%
  mutate(moneda_obs = moneda) %>%
  select(moneda_obs, theta, prob_inicial, verosimilitud, prob_posterior)
```

```{block, type='ejercicio'}

- ¿Qué pasa cuando el número de soles es 0? ¿Cómo cambian las probabilidades
posteriores de cada moneda?
- Incrementa el número de volados, por ejemplo a 10. ¿Qué pasa si observaste
8 soles, por ejemplo? ¿Y si observaste 0?
- ¿Qué pasa si cambias las probabilidades iniciales (por ejemplo incrementas
la probabilidad inicial de la moneda 1 a 0.9)?

```

Justifica las siguientes aseveraciones (para este ejemplo):

```{block, type='ejercicio'}

- Las probabilidades posteriores o finales son una especie de punto intermedio
entre verosimilitud y probablidades iniciales.
- Si tenemos pocas observaciones, las probabilidades posteriores son similares
a las iniciales.
- Cuando tenemos muchos datos, las probabilidades posteriores están más
concentradas, y no es tan importante la inicial.
- Si la inicial está muy concentrada en algún valor, la posterior requiere de
muchas observaciones para que se pueda concentrar en otros valores diferentes a
los de la inicial.

```

Ahora resumimos los elementos básicos de la inferencia bayesiana, que son
relativamente simples:

```{block, type='mathblock'}
**Inferencia bayesiana.** Con la notación de arriba:

- Como $X$ son los datos observados, llamamos a $P(X|\theta)$ la *verosimilitud*,
proceso generador de datos o modelo de los datos.
- El factor $P(\theta)$ le llamamos la distribución *inicial* o *previa*.
- La distribución que usamos para hacer inferencia $P(\theta|X)$ es la
distribución *final* o *posterior*

Hacemos inferencia usando la ecuación

$$P(\theta | X) = \frac{P(X | \theta) P(\theta)}{P(X)} = \frac{P(X | \theta) P(\theta)}{\sum_{\theta \in \Theta} P(X | \theta) P(\theta)}\,,$$

que también escribimos:

$$P(\theta | X) \propto P(X | \theta) P(\theta)\,,$$
donde $\propto$ significa "proporcional a". No ponemos $P(X)$ pues como vimos
arriba, es una constante de normalización (y en la mayoría de los casos es muy dificil poder derivar analiticamente).

```


En estadística Bayesiana, las probablidades posteriores $P(\theta|X)$ dan toda
la información que necesitamos para hacer inferencia. ¿Cuándo damos probablidad
alta a un parámetro particular $\theta$? Cuando su verosimilitud es alta y/o
cuando su probabilidad inicial es alta. De este modo, la posterior combina la
información inicial que tenemos acerca de los parámetros con la información en
la muestra acerca de los parámetros (verosimilitud). Podemos ilustrar como
sigue:

```{r, echo = FALSE, out.width= "60%"}
knitr::include_graphics('images/perros.png')
```


## Ejemplo: estimando una proporción {-}

Regresamos ahora a nuestro problema de estimar una proporción $\theta$ de una
población dada usando una muestra iid $X_1,X_2,\ldots, X_n$ de variables
Bernoulli.  Ya sabemos calcular la
verosimilitud (el modelo de los datos):

$$P(X_1=x_1,X_2 =x_2,\ldots, X_n=x_n|\theta) = {n \choose k} \theta^k(1-\theta)^{n-k}\,,$$

donde $k = x_1 + x_2 +\cdots + x_k$ es el número de éxitos que observamos.

Ahora necesitamos una distribución inicial o previa $P(\theta)$. Aunque esta
distribución puede tener cualquier forma, supongamos que nuestro conocimiento
actual podemos resumirlo con una distribución $\mathsf{Beta}(3, 3)$:

$$P(\theta) = \frac{1}{B(3,3)} \theta^{3-1}(1-\theta)^{3-1}\,,$$
donde $B(a, b)$ es la constante de normalización de una distribución 
$\mathsf{Beta}$. Podemos simular para examinar su forma:

```{r}
params.inicial <- list(a = 3, b = 3)
sim_inicial <- tibble(theta = rbeta(10000, params.inicial$a, params.inicial$b))
ggplot(sim_inicial) + geom_histogram(aes(x = theta, y = ..density..), bins = 15)
```
De modo que nuestra información inicial es que la proporción puede tomar
cualquier valor entre 0 y 1, pero es probable que tome un
valor no tan lejano de 0.5. Por ejemplo, con probabilidad 0.95 creemos
que $\theta$ está en el intervalo
```{r}
quantile(sim_inicial$theta, c(0.025, 0.975)) %>% round(2)
```
Es difícil justificar en abstracto por qué escogeriamos una inicial con esta
forma. *Aunque esto los detallaremos más adelante*, puedes pensar, por el
momento, que alguien observó algunos casos de esta población, y quizá vio tres éxitos y tres fracasos. Esto sugeriría que es poco probable que la probablidad
$\theta$ sea muy cercana a 0 o muy cercana a 1.

Ahora podemos construir nuestra posterior. Tenemos que

$$P(\theta| X_1=x_1, \ldots, X_n=x_n) \propto P(X_1 = x_1,\ldots
X_n=x_n | \theta)P(\theta) = \theta^{k+2}(1-\theta)^{n-k + 2}\,,$$

donde la constante de normalización no depende de $\theta$. Como
$\theta$ es un parámetro continuo, la expresión de la derecha nos debe
dar una función de densidad.

Supongamos entonces que hicimos la prueba con $n = 30$ (número de
prueba) y observamos 19 éxitos. Tendríamos que 

$$P(\theta | S_n = 19) \propto \theta^{19 + 2} (1-\theta)^{30 - 19 +2} = \theta^{21}(1-\theta)^{13}\,.$$

La cantidad de la derecha, una vez que normalizemos por el número $P(X=19)$, nos
dará una densidad posterior (tal cual, esta expresion no integra a 1). Podemos
obtenerla usando cálculo, pero recordamos que una distribución
$\mathsf{\mathsf{Beta}}(a,b)$ tiene como fórmula
$$\frac{1}{B(a,b)} \theta^{a-1}(1-\theta)^{b-1}\,.$$
Concluimos entonces que la posterior tiene una distribución $\mathsf{Beta}(22,
14)$. Podemos simular de la posterior usando código estándar para ver cómo luce:

```{r}
dbinom_verosimilitud <- function(data){
  function(theta, scale = 1){
    scale * dbinom(x = data$x, size = data$n, prob = theta)
  }
}

verosimilitud <- dbinom_verosimilitud(tibble(x = 19, n = 30))
```

```{r}
sim_inicial <- sim_inicial %>% mutate(dist = "inicial")
sim_posterior <- tibble(theta = rbeta(10000, 
                                      params.inicial$a + 19, 
                                      params.inicial$b + 30 - 19)) %>% 
  mutate(dist = "posterior")

sims <- bind_rows(sim_inicial, sim_posterior)
ggplot(sims, aes(x = theta, fill = dist)) +
  geom_histogram(aes(x = theta), bins = 30, 
                 alpha = 0.5, position = "identity") + 
  geom_function(aes(group = 'verosimilitud'), fun = verosimilitud, 
                args = list(scale = 1e4), lty = 2, show.legend = TRUE)
```
La posterior nos dice cuáles son las *posibilidades* de dónde puede estar
el parámetro $\theta$. Nótese que ahora excluye prácticamente valores más chicos
que 0.25 o mayores que 0.9. Esta distribución posterior es el objeto con el
que hacemos inferencia: nos dice dónde es creíble que esté el parámetro $\theta$.

Podemos resumir de varias maneras. Por ejemplo, si queremos un estimador puntual
usamos la media posterior:

```{r}
sims %>% group_by(dist) %>%
  summarise(theta_hat = mean(theta) %>% round(3))
```
Nota que el estimador de máxima verosimilitud es $\hat{p} = 19/30 = 0.63$, que
es ligeramente diferente de la media posterior. ¿Por qué?

Y podemos construir intervalos de percentiles, que en esta situación
suelen llamarse *intervalos de credibilidad*, por ejemplo:

```{r}
f <- c(0.025, 0.975)
sims %>% group_by(dist) %>%
  summarise(cuantiles = quantile(theta, f) %>% round(2), f = f) %>%
  pivot_wider(names_from = f, values_from = cuantiles)
```
El segundo renglón nos da un intervalo posterior para $\theta$ de *credibilidad*
95\%. En inferencia bayesiana esto sustituye a los intervalos de confianza.

- El **intervalo de la inicial** expresa nuestras creencias a priori
acerca de $\theta$. Este intervalo es muy amplio (va de 0.15 a 0.85)
- El **intervalo de la posterior** actualiza nuestras creencias acerca
de $\theta$ una vez que observamos los datos, y es considerablemente
más angosto y por lo tanto informativo.


**Observaciones**:

- Nótese que escogimos una forma analítica fácil para la inicial, pues resultó
así que la posterior es una distribución beta. No siempre es así, y veremos qué
hacer cuando nuestra inicial no es de un tipo "conveniente".
- Como tenemos la forma analítica de la posterior, es posible hacer los cálculos
de la media posterior, por ejemplo, **integrando** la densidad posterior a mano. Esto
generalmente no es factible, y en este ejemplo preferimos hacer una aproximación numérica. En este caso
particular es posible usando cálculo, y sabemos que la media de una $\mathsf{\mathsf{Beta}}(a,b)$ es
$a/(a+b)$, de modo que nuestra media posterior es

$$\hat{\mu} = (19 + 3)/(30 + 6) = 22/36 = 0.611\,, $$
que podemos interpretar como sigue: para calcular la media posterior, a nuestras
$n$ pruebas iniciales agregamos 6 pruebas adicionales fijas, con 3 éxitos y 3
fracasos, y calculamos la proporción usual de éxitos.

```{block, type='ejercicio'}
Repite el análisis considerando en general $n$ pruebas, con $k$ éxitos. Utiliza la
misma distribución inicial.
```
- Lo mismo aplica para el intervalo de 95% (¿cómo se calcularía integrando?). También
puedes usar la aproximación de R, por ejemplo:

```{r}
qbeta(0.025, shape1 = 22, shape2 = 14) %>% round(2)
qbeta(0.975, shape1 = 22, shape2 = 14) %>% round(2)
```

**Predictivas**:

Ahora escribiremos la función predictiva previa para este caso. La
derivación de la predictiva posterior se deja como ejercicio
adicional. Primero notemos lo siguiente. La previa corresponde a una
variable aleatoria Beta (la posterior también). Esto aliviará nuestros
cálculos.

Para la predictiva previa tenemos que calcular la integral
$$P(X = k) = \int_0^1 P(X = k|\theta) P(\theta) d\theta\,,$$
lo cual cual nos lleva a
\begin{align}
P(X = k) &= \int_{0}^{1} {n \choose k} \theta^{k}(1- \theta)^{n-k} \frac{1}{B(a,b)} \theta^{a-1} (1-\theta)^{b-1} d\theta\\
&= {n \choose k} \frac{1}{B(a,b)} \int_{0}^{1} \theta^{k + a -1}(1- \theta)^{n-k + b -1} d\theta\\
&= {n \choose k} \frac{B(k + a, n-k+b)}{B(a,b)} \,,
\end{align}
donde hemos utilizado lo que sabemos del *núcleo* de una
$\mathsf{Beta}$.  La distribución resultante es parte de un catálogo
*extenso* de distribuciones. De hecho, decimos que $X \in \mathbb{N}$
condicionado en $n$ tiene una distribución $\mathsf{Beta-Binomial}$ si
su función de densidad es igual a la expresión de arriba. Como vimos
esta distribución tiene su origen al *mezclar* una $\mathsf{Beta}$ y
una $\mathsf{Binomial}$.

En el contexto de nuestro ejemplo, al usar una $\mathsf{Beta}(3,3)$ podemos
explorar el número de éxitos que creemos probables bajo esta distribución
previa. Para esto podemos simular como sigue:

```{r}

rbetabinom <- function(parametros){
  tibble(theta  = rbeta(n = 1, parametros$a, parametros$b), 
         exitos = rbinom(n = 1, size = parametros$n, prob = theta))
}

tibble(id = seq(1, 5000), 
       n = 30,             # Número de lanzamientos
       a = 3, b = 3) |>    # Parámetros de la beta
  nest(-id) |> 
  mutate(samples = map(data, rbetabinom)) |> 
  unnest(samples) |> 
  ggplot(aes(x = exitos)) + 
    geom_histogram(binwidth = 2)

```

Lo cual nos da otro mecanismo de validar (antes de analizar los datos) si 
nuestras creencias previas tienen sentido **en el contexto** del problema. 
Es decir, podríamos verificar si tiene sentido que el promedio de éxitos se 
concentre alrededor de 15 cuando estamos recolectando 30 experimentos totales. 
También podemos expresar los intervalos de credibilidad sobre el número de
éxitos en 30 experimentos:

```{r}

tibble(id = seq(1, 5000), 
       n = 30,             # Número de lanzamientos
       a = 3, b = 3) |>    # Parámetros de la beta
  nest(-id) |> 
  mutate(samples = map(data, rbetabinom)) |> 
  unnest(samples) |> 
  pull(exitos) |> 
  quantile(probs = c(.025, .975))

```

```{block, type = "ejercicio"}

Repite lo anterior con la distribución predictiva posterior considerando que hemos observado 
19 éxitos en 30 experimentos, y que queremos realizar 30 experimentos más. Utiliza como 
previa una $\mathsf{Beta}(3,3)$. ¿Cómo se ve el histograma y el intervalo de credibilidad 
posterior?

```


## Ejemplo: observaciones uniformes {-}

Ahora regresamos al problema de estimación del máximo de una distribución uniforme.
En este caso, consideraremos un problema más concreto. Supongamos que hay una lotería
(tipo tradicional)
en México y no sabemos cuántos números hay. Obtendremos una muestra iid de $n$ números,
ya haremos una aproximación continua, suponiendo que

$$X_i \sim U[0,\theta]\,.$$

La verosimilitud es entonces (recuerda la ecuación (6.5))
$$P(X_1,\ldots, X_n|\theta) = \theta^{-n} 1_{[X_\max,\infty)}(\theta)\,,$$
cuando $\theta$ es mayor que todas las $X_i$ (o bien, cuando $\theta$ es mayor
que el máximo de las observaciones $X_\max$), y cero en otro caso. Necesitaremos
una inicial $P(\theta)$.

Por la forma que tiene la verosimilitud, podemos intentar una distribución Pareto,
que tiene la expresióñ
$$P(\theta) = \frac{\alpha \theta_0^\alpha}{\theta^{\alpha + 1}} 1_{[\theta_0,\infty]}(\theta) \,,$$
con soporte en $[\theta_0,\infty]$ (esto se refleja al incorporar la función
indicadora). Tenemos que escoger entonces el mínimo $\theta_0$ y el parámetro
$\alpha$. En primer lugar, como sabemos que es una lotería nacional, creemos que
no puede haber menos de unos 300 mil números, así que $\theta_0 = 300$. La
función acumulada de la pareto es $1- (300/\theta)^\alpha$, así que el cuantil
99% es

```{r}
alpha <- 1.1
(300/(0.01)^(1/alpha))
```
es decir, alrededor de 20 millones de números.  Creemos que es un poco probable
que el número de boletos sea mayor a esta cota.
Nótese ahora que la posterior cumple (multiplicando verosimilitud por inicial):

$$P(\theta|X_1,\ldots, X_n |\theta) \propto \theta^{-(n + 2.1)} \left(1_{[X_\max,\infty]}(\theta)  \times 1_{[\theta_0,\infty]}(\theta)\right) \,,$$

para $\theta$ mayor que el máximo de las $X_n$'s y $\theta_0 = 300$, y cero en
otro caso. **Recuerda** que $1_A(x) \times 1_B(x) = 1_{A\cap B}(x)$. Esta
distribución es Pareto con $\theta_0' = \max\{300, X_1,\ldots, X_n\}$ y $\alpha
= n + 1.1$

Una vez planteado nuestro modelo, veamos los datos. Obtuvimos la siguiente
muestra de números:
```{r}
loteria_tbl <- read_csv("data/nums_loteria_avion.csv", col_names = c("id", "numero")) %>%
  mutate(numero = as.integer(numero))
set.seed(334)
muestra_loteria <- sample_n(loteria_tbl, 25) %>%
  mutate(numero = numero/1000)
muestra_loteria %>% as.data.frame %>% head
```

Podemos simular de una Pareto como sigue:

```{r}
rpareto <- function(n, theta_0, alpha){
  # usar el método de inverso de distribución acumulada
  u <- runif(n, 0, 1)
  theta_0 / (1 - u)^(1/alpha)
}
```

Simulamos de la inicial:

```{r}
sims_pareto_inicial <- tibble(
  theta = rpareto(20000, 300, 1.1 ),
  dist = "inicial")
```

Y con los datos de la muestra, simulamos de la posterior:

```{r}
sims_pareto_posterior <- tibble(
  theta = rpareto(20000,
                  max(c(300, muestra_loteria$numero)),
                  nrow(muestra_loteria) + 1.1),
  dist = "posterior")
sims_theta <- bind_rows(sims_pareto_inicial, sims_pareto_posterior)
ggplot(sims_theta) +
  geom_histogram(aes(x = theta, fill = dist),
                 bins = 70, alpha = 0.5, position = "identity",
                 boundary = max(muestra_loteria$numero))  +
  xlim(0, 15000) + scale_y_sqrt() +
  geom_rug(data = muestra_loteria, aes(x = numero))
```
Nótese que cortamos algunos valores de la inicial en la cola derecha: un defecto
de esta distribución inicial, con una cola tan larga a la derecha, es que
pone cierto peso en valores que son poco creíbles y la vuelve poco apropiada para
este problema. Regresamos más adelante a este problema.

Si obtenemos percentiles,
obtenemos el intervalo

```{r}
f <- c(0.01, 0.5, 0.99)
sims_theta %>% group_by(dist) %>%
  summarise(cuantiles = quantile(theta, f) %>% round(2), f = f) %>%
  pivot_wider(names_from = f, values_from = cuantiles)
```
Estimamos entre 5.8 millones y 6.7 millones de boletos. El máximo en la muestra
es de
```{r}
max(muestra_loteria$numero)
```

Escoger la distribución pareto como inicial es conveniente y nos permitió
resolver el problema sin dificultad, pero por su forma vemos que no
necesariamente es apropiada para el problema por lo que señalamos arriba.
Nos gustaría, por ejemplo, poner una inicial como
la siguiente

```{r}
qplot(rgamma(2000, 5, 0.001), geom="histogram", bins = 20) +
  scale_x_continuous(breaks = seq(1000, 15000, by = 2000))
```

Sin embargo, los cálculos no son tan simples en este caso, pues la posterior
no tiene un forma reconocible. Tendremos que usar otras estrategias de simulación
para ejemplos como este (integración Monte Carlo por medio de Cadenas de Markov,
que veremos en el próximo curso).

## Probabilidad a priori {-}

La inferencia bayesiana es conceptualmente simple: siempre hay que
calcular la posterior a partir de verosimilitud (modelo de datos) y
distribución inicial o a priori. Sin embargo, una crítica usual que se
hace de la inferencia bayesiana es precisamente que hay que tener esa
información inicial, y que distintos analistas llegan a distintos
resultados si tienen información inicial distinta.

Eso realmente no es un defecto, es una ventaja de la inferencia
bayesiana. Los datos y los problemas que queremos resolver no viven en
un vacío donde podemos creer que la estatura de las personas, por
ejemplo, puede variar de 0 a mil kilómetros; o el número de boletos de
una lotería puede variar de 2 a 3 boletos o quizá 500 millones de
boletos; incluso la proporción de personas infectadas de una
enfermedad puede ser de unos cuantos hasta miles de millones.

- En todos estos casos tenemos cierta información inicial que podemos
usar para informar nuestras estimaciones. Esta información debe
usarse.
- Antes de tener datos, las probabilidades iniciales deben ser
examinadas en términos del conocimiento de expertos.
- Las probabilidades iniciales son supuestos que hacemos acerca del
problema de interés, y también están sujetas a críticas y
confrontación con datos.

Las fuentes usuales de información que podemos utilizar varía de
acuerdo al estudio. Por ejemplo, puede provenir de datos
históricos. Incluso puede provenir de datos en la literatura o de
experimentos relacionados.

## Análisis conjugado {-}

Los dos ejemplos que hemos visto arriba son ejemplos de análisis conjugado:

- ($\textsf{Beta-bernoulli}$) Si las observaciones $X_i$ son $\mathsf{Bernoulli}(p)$ ($n$ fija)
 queremos estimar $p$, y tomamos como distribución inicial para $p$ una $\mathsf{Beta}(a,b)$,
entonces la posterior para $p$ cuando $S_n=k$ es $\mathsf{Beta}(k + a, n - k + b)$,
donde $S_n = X_1 + X_2 +\cdots +X_n$.

Y más en general:

- ($\textsf{Beta-binomial}$) Si las observaciones $X_i, i=1,2,\ldots, m$
son $\mathsf{Binomial}(n_i, p)$ ($n_i$'s fijas) independientes,
 queremos estimar $p$, y tomamos como distribución inicial para $p$ una $\mathsf{Beta}(a,b)$,
entonces la posterior para $p$ cuando $S_m=k$ es $\mathsf{Beta}(k + a, n - k + b)$,
donde $S_m = X_1 + X_2 +\cdots +X_m$ y $n= n_1+n_2+\cdots+n_m$

También aplicamos:

- ($\textsf{Uniforme-Pareto}$) Si el modelo de datos $X_i$ es uniforme $\mathsf{U}[0,\theta]$ ($n$ fija),
queremos estimar $\theta$, y tomamos como distribución inicial para $\theta$ una
Pareto $(\theta_0, \alpha)$, entonces la posterior para $p$ si el máximo de las $X_i$'s
es igual a $M$ es Pareto con parámetros $(\max\{\theta_0, M\}, \alpha + n)$.

Nótese que en estos casos, dada una forma de la verosimilitud, tenemos una
familia conocida de iniciales tales que las posteriores están en la misma
familia. Estos modelos son convenientes porque podemos hacer simulaciones de la
posterior de manera fácil, o usar sus propiedades teóricas.

```{block, type = 'mathblock'}

**Definición.** Un *modelo bayesiano conjugado* es aquel que mantiene la
distribución posterior dentro de la misma familia que la distribución previa. 

Es decir, sea una muestra $X_1, \ldots, X_n \overset{iid}{\sim} f(X\,|\,
\theta),$ y sea $p (\theta ) = \pi(\theta \, | \, \varphi_0)$ 
la distribución previa para nuestro parámetro desconocido $\theta$ donde
$\varphi_0$ denota los hiper-parámetros que definen dicha distribución. 

Bajo un modelo conjugado, la distribución posterior $P(\theta | X_1, \ldots, X_n)$ 
satisface 
$$P(\theta \, | \, X_1, \ldots, X_n) = \pi(\theta \, | \, \varphi_1)\,,$$
donde $\varphi_1 = h(\varphi_0, X_1, \ldots, X_n)$ son los hiper-parámetros 
actualizados del modelo que se usó para la previa.

```


Otro ejemplo típico es el modelo normal-normal:

- ($\textsf{Normal-normal}$) Si $X_i\sim \mathsf{N}(\mu,\sigma)$, con $\sigma$ conocida, y tomamos
como distribución inicial para $\mu \sim \mathsf{N}(\mu_0,\sigma_0)$, y definimos
la *precisión* $\tau$ como el inverso de la varianza $\sigma^2$, entonces la posterior
de $\mu$ es Normal con media $(1-\lambda) \mu_0 + \lambda\overline{x}$,
y precisión $\tau_0 + n\tau$, donde $\lambda = \frac{n\tau}{\tau_0 + n\tau}$

```{block, type='ejercicio'}
Completa cuadrados para mostrar las fórmulas del modelo normal-normal con
varianza conocida.
```

Más útil es el siguiente modelo:

- ($\textsf{Normal-Gamma}$ inverso) Sean $X_i\sim \mathsf{N}(\mu, \sigma)$. Queremos estimar $\mu$ y $\sigma$. Tomamos
como distribuciones iniciales (dadas por 4 parámetros: $\mu_0, n_0, \alpha,\beta$):

  - $\tau = \frac{1}{\sigma^2} \sim \mathsf{Gamma}(\alpha,\beta)$
  - $\mu|\sigma$ es normal con media $\mu_0$ y varianza $\sigma^2 / {n_0}$ , y
  - $P(\mu, \sigma) = P(\mu|\sigma)P(\sigma)$

- Entonces la posterior es:
  - $\tau|x$ es $\mathsf{Gamma}(\alpha', \beta')$, con $\alpha' = \alpha + n/2$,
  $\beta' = \beta + \frac{1}{2}\sum_{i=1}^{n}(x_{i} - \bar{x})^2 + \frac{nn_0}{n+n_0}\frac{({\bar{x}}-\mu_{0})^2}{2}$.
  - $\mu|\sigma,x$ es normal con media $\mu' = \frac{n_0\mu_{0}+n{\bar{x}}}{n_0 +n}$ y varianza $\sigma^2/({n_0 +n})$.

  - $P(\mu,\sigma|x) = P(\mu|x,\sigma)P(\sigma|x)$

```{block , type = 'mathblock'}

**Definición.** Decimos que la distribución conjunta de dos variables aleatorias 
$\theta = (\mu, \tau)$ sigue una distribución Normal--Gamma-Inversa, si podemos
establecer la siguiente estructura de dependencia condicional
$$\pi(\mu, \tau ) = \pi(\mu \, | \, \tau) \, \pi(\tau),$$
donde los factores siguen las siguientes distribuciones
\begin{align*}
\mu \, | \, \tau , \mu_0, n_0   \sim \mathsf{N}\left(\mu_0, \frac{1}{\tau n_0}\right), \qquad \text{y} \qquad \tau \, | \, \alpha_0, \beta_0  \sim \textsf{Gamma}(\alpha_0, \beta_0)\,.
\end{align*}

```



**Observaciones**

1. Nótese que este último ejemplo tiene más de un parámetro. En estos casos, el
objeto de interés es la **posterior conjunta** de los parámetros
$P(\theta_1,\theta_2,\cdots, \theta_p|x)$. Este último ejemplo es relativamente
simple pues, por la selección de iniciales, para simular de la conjunta de $\mu$
y $\tau$ podemos simular primero $\tau$ (o $\sigma$), y después usar este valor
para simular de $\mu$: el par de valores resultantes son una simulación de la
conjunta.

2. Los parámetros $a,b$ para la inicial de $\tau$ pueden interpretarse como
sigue: $\sqrt{b/a}$ es un valor "típico" a priori para la varianza poblacional,
y $a$ indica qué tan seguros estamos de este valor típico.

3. Nótese que para que funcionen las fórmulas de la manera más simple, escogimos
una dependencia a priori entre la media y la precisión: $\tau = \sigma^{-2}$
indica la escala de variabilidad que hay en la población, la incial de la media
tiene varianza $\sigma^2/n_0$. Si la escala de variabilidad de la población es
más grande, tenemos más incertidumbre acerca de la localización
de la media.

4. Aunque esto tiene sentido en algunas aplicaciones, y por convenviencia usamos
esta familia conjugada, muchas veces es preferible otro tipo de especificaciones
para las iniciales: por ejemplo, la media normal y la desviación estándar
uniforme, o media normal, con iniciales independientes. Sin embargo, estos casos
no son tratables con análisis conjugado (veremos más adelante cómo tratarlos con
MCMC).


### Ejemplo {-}
Supongamos que queremos estimar la estatura de los cantantes de tesitura tenor
con una muestra iid de tenores de Estados Unidos. Usaremos el modelo normal de
forma que $X_i\sim \mathsf{N}(\mu, \sigma^2)$.

Una vez decidido el modelo, tenemos que poner distribución inicial para los
parámetros $(\mu, \sigma^2)$.

Comenzamos con $\sigma^2$. Como está el modelo, esta inicial debe estar dada
para la precisión $\tau$, pero podemos simular para ver cómo se ve nuestra
inicial para la desviación estándar. En la población general la desviación
estándar es alrededor de 7 centímetros.

```{r}
# Comenzamos seleccionando un valor que creemos típico para la desviación estándar
sigma_0 <- 7
# seleccionamos un valor para a, por ejemplo: si es más chico sigma tendrá más
# disperisón
a <- 3
# ponemos 8 = sqrt(b/a) -> b = a * 64
b <- a * sigma_0^2
c(a = a, b = b)
```

Ahora simulamos para calcular cuantiles

```{r}
tau <- rgamma(1000, a, b)
quantile(tau, c(0.05, 0.95))
sigma <- 1 / sqrt(tau)
c(media = 1/sqrt(mean(tau)), media_sigma = mean(sigma), quantile(sigma, c(0.05, 0.95)))
```
Que es dispersión considerable: con poca probabilidad la desviación estándar
es menor a 4 centímetros, y también creemos que es poco creíble la desviación
estándar sea de más de 13 centímetros.


Comenzamos con $\mu$. Sabemos, por ejemplo, que con alta probabilidad la media
debe ser algún número entre 1.60 y 1.80. Podemos investigar: la media nacional
en estados unidos está alrededor de 1.75, y el percentil 90% es 1.82. Esto es
*variabilidad en la población*: debe ser muy poco probable, por ejemplo, que la
media de tenores sea 1.82 Quizá los cantantes tienden a ser un poco más altos o
bajos que la población general, así que podríamos agregar algo de dispersión.

Podemos establecer parámetros y simular de la marginal a partir de las fórmulas
de arriba para entender cómo se ve la inicial de $\mu$:

```{r}
mu_0 <- 175 # valor medio de inicial
n_0 <- 5    # cuánta concentración en la inicial
tau <- rgamma(1000, a,b)
sigma <- 1/sqrt(tau)
mu <- map_dbl(sigma, ~ rnorm(1, mu_0, .x / sqrt(n_0)))
quantile(mu, c(0.05, 0.5, 0.95))
```
Que consideramos un rango en el que con alta probabilidad debe estar
la media poblacional de los cantantes.

Podemos checar nuestros supuestos simulando posibles muestras usando
sólo nuestra información previa:

```{r, fig.width = 7, fig.height = 6}
simular_normal_invgamma <- function(n, pars){
  mu_0 <- pars[1]
  n_0 <- pars[2]
  a <- pars[3]
  b <- pars[4]
  # simular media
  tau <- rgamma(1, a, b)
  sigma <- 1 / sqrt(tau)
  mu <- rnorm(1, mu_0, sigma/sqrt(n_0))
  # simular sigma
  rnorm(n, mu, sigma)
}
set.seed(3461)
sims_tbl <- tibble(rep = 1:20) %>%
  mutate(estatura = map(rep, ~ simular_normal_invgamma(500, c(mu_0, n_0, a, b)))) %>%
  unnest(cols = c(estatura))
ggplot(sims_tbl, aes(x = estatura)) + 
  geom_rect(inherit.aes = FALSE,
            aes(xmin=150, xmax=180, ymin=0, ymax=+Inf), 
            fill='gray', alpha=0.1) + 
  geom_histogram() +
  facet_wrap(~ rep) +
  ggtitle("Realizaciones de datos bajo la previa")
```

Pusimos banda de referencia entre alturas de 150 y 180. Vemos que nuestras
iniciales no producen simulaciones totalmente fuera del contexto, y parecen
cubrir apropiadamente el espacio de posiblidades para estaturas de los tenores.
Quizá hay algunas realizaciones poco creíbles, pero no extremadamente. En este
punto, podemos regresar y ajustar la inicial para $\sigma$, que parece tomar
valores demasiado grandes (produciendo por ejemplo una simulación con estatura
de 220 y 140, que deberían ser menos probables).

Ahora podemos usar los datos para calcular nuestras posteriores.

```{r}
set.seed(3413)
cantantes <- lattice::singer %>%
  mutate(estatura_cm = round(2.54 * height)) %>%
  filter(str_detect(voice.part, "Tenor")) %>%
  sample_n(20)
cantantes
```

Los cálculos son un poco tediosos, pero podemos construir una función apropiada:

```{r}
calcular_pars_posterior <- function(x, pars_inicial){
  # iniciales
  mu_0 <- pars_inicial[1]
  n_0 <- pars_inicial[2]
  a_0 <- pars_inicial[3]
  b_0 <- pars_inicial[4]
  # muestra
  n <- length(x)
  media <- mean(x)
  S2 <- sum((x - media)^2)
  # sigma post
  a_1 <- a_0 + 0.5 * n
  b_1 <- b_0 + 0.5 * S2 + 0.5 * (n * n_0) / (n + n_0) * (media - mu_0)^2
  # posterior mu
  mu_1 <- (n_0 * mu_0 + n * media) / (n + n_0)
  n_1 <- n + n_0
  c(mu = mu_1, n = n_1, alpha = a_1, beta = b_1)
}
pars_posterior <- calcular_pars_posterior(cantantes$estatura_cm, c(mu_0, n_0, a, b))
pars_posterior
```

¿Cómo se ve nuestra posterior comparada con la inicial? Podemos hacer simulaciones:

```{r}
sim_params <- function(m, pars){
  mu_0 <- pars[1]
  n_0 <- pars[2]
  a <- pars[3]
  b <- pars[4]
  # simular sigmas
  sims <- tibble(tau = rgamma(m, a, b)) %>%
    mutate(sigma = 1 / sqrt(tau))
  # simular mu
  sims <- sims %>% mutate(mu = rnorm(m, mu_0, sigma / sqrt(n_0)))
  sims
}
sims_inicial <- sim_params(5000, c(mu_0, n_0, a, b)) %>%
  mutate(dist = "inicial")
sims_posterior <- sim_params(5000, pars_posterior) %>%
  mutate(dist = "posterior")
sims <- bind_rows(sims_inicial, sims_posterior)
ggplot(sims, aes(x = mu, y = sigma, colour = dist)) +
  geom_point()
```

Y vemos que nuestra posterior es consistente con la información inicial
que usamos, hemos aprendido considerablemente de la muestra. La posterior se
ve como sigue. Hemos marcado también las medias posteriores de cada parámetro:
media y desviación estándar.

```{r}
medias_post <- sims %>% filter(dist == "posterior") %>%
  select(-dist) %>%
  summarise(across(everything(), mean))
ggplot(sims %>% filter(dist == "posterior"),
    aes(x = mu, y = sigma)) +
  geom_point(colour = "#00BFC4") +
  geom_point(data = medias_post, size = 5, colour = "black") +
  coord_equal()
```
Podemos construir intervalos de credibilidad del 90% para estos dos parámetros,
por ejemplo haciendo intervalos de percentiles:

```{r}
f <- c(0.05, 0.5, 0.95)
sims %>%
  pivot_longer(cols = mu:sigma, names_to = "parametro") %>%
  group_by(dist, parametro) %>%
  summarise(cuantil = quantile(value, f) %>% round(1), f= f) %>%
  pivot_wider(names_from = f, values_from = cuantil)
```

Como comparación, los estimadores de máxima verosimlitud son

```{r}
media_mv <- mean(cantantes$estatura_cm)
sigma_mv <- mean((cantantes$estatura_cm - media_mv)^2) %>% sqrt
c(mu_mle = media_mv, sigma_mle = sigma_mv)
```

Ahora solo resta checar que el modelo es razonable. Veremos más adelante cómo
hacer esto, usando la distribución predictiva posterior.


## Pasos de un análisis de datos bayesiano {-}

```{block, type='comentario'}
Como vimos en los ejemplos, en general un análisis de datos bayesiano sigue los
siguientes pasos:

- Identificar los datos releventes a nuestra pregunta de investigación, el tipo
de datos que vamos a describir, que variables queremos estimar.

- Definir el modelo descriptivo para los datos. La forma matemática y los
parámetros deben ser apropiados para los objetivos del análisis.

- Especificar la distribución inicial de los parámetros.

- Utilizar inferencia bayesiana para reubicar la credibilidad a lo largo de
los posibles valores de los parámetros.

- Verificar que la distribución posterior replique los datos de manera
razonable, de no ser el caso considerar otros modelos descriptivos para los datos.
```



#### Elicitando probablidades subjetivas (opcional) {-}

No siempre es fácil elicitar probabilidades subjetivas de manera que capturemos
el verdadero conocimiento de dominio que tenemos. Una manera clásica de hacerlo
es con apuestas

Considera una pregunta sencilla que puede afectar a un viajero: ¿Qué tanto
crees que habrá una tormenta que ocasionará el cierre de la autopista
México-Acapulco en el puente del $20$ de noviembre? Como respuesta debes dar
un número entre $0$ y $1$ que refleje tus creencias. Una manera de seleccionar
dicho número es calibrar las creencias en relación a otros eventos cuyas
probabilidades son claras.

Como evento de comparación considera una experimento donde hay canicas en una
urna: $5$ rojas y $5$ blancas. Seleccionamos una canica al azar. Usaremos esta urna
como comparación para considerar la tormenta en la autopista. Ahora, considera
el siguiente par de apuestas de las cuales puedes elegir una:

* A. Obtienes $\$1000$ si hay una tormenta que ocasiona el cierre de la autopista
el próximo $20$ de noviembre.

* B. Obtienes $\$1000$ si seleccionas una canica roja de la urna que contiene
$5$ canicas rojas y $5$ blancas.

Si prefieres la apuesta B, quiere decir que consideras que la probabilidad de
tormenta es menor a $0.5$, por lo que al menos sabes que tu creencia subjetiva de
una la probabilidad de tormenta es menor a $0.5$. Podemos continuar con el proceso
para tener una mejor estimación de la creencia subjetiva.

* A. Obtienes $\$1000$ si hay una tormenta que ocasiona el cierre de la autopista
el próximo $20$ de noviembre.

* C. Obtienes $\$1000$ si seleccionas una canica roja de la urna que contiene
$1$ canica roja y $9$ blancas.

Si ahora seleccionas la apuesta $A$, esto querría decir que consideras que la
probabilidad de que ocurra una tormenta es mayor a $0.10$. Si consideramos ambas
comparaciones tenemos que tu probabilidad subjetiva se ubica entre $0.1$ y $0.5$.

## Verificación predictiva posterior {-}

Una vez que ajustamos un modelo bayesiano, podemos simular nuevas observaciones
a partir del modelo. Esto tiene dos utilidades:

- Hacer predicciones acerca de datos no observados.
- Confirmar que nuevas producidas simuladas con el modelo son similares a las
que de hecho observamos. Esto nos permite confirmar la calidad del ajuste del
modelo, y se llama **verificación predictiva posterior**.

Supongamos que tenemos la posterior $P(\theta | x)$. Podemos generar una nueva
*replicación* de los datos como sigue:

```{block, type='mathblock'}
La distribución **predictiva posterior** genera nuevas observaciones a partir de
la información observada. La denotamos como $P(\tilde{x}|x)$.

Para simular de ella:

- Muestreamos un valor $\tilde{\theta}$ de la posterior $P(\theta|x)$.
- Simulamos del modelo de las observaciones $\tilde{x} \sim P(\tilde{x}|\tilde{\theta})$.
- Repetimos el proceso hasta obtener una muestra grande.
- Usamos este método para producir, por ejemplo, **intervalos de predicción** para nuevos datos.

Si queremos una replicación de las observaciones de la predictiva posterior,

- Muestreamos un valor $\tilde{\theta}$ de la posterior $P(\theta|x)$.
- Simulamos del modelo de las observaciones $\tilde{x}_1, \tilde{x}_2,\ldots, \tilde{x}_n \sim P(\tilde{x}|\tilde{\theta})$,
done $n$ es el tamaño de muestra de la muestra original $x$.
- Usamos este método para producir conjuntos de datos simulados que comparamos con los observados para verificar nuestro modelo.

```


### Ejemplo: estaturas de tenores {-}

En este ejemplo, usaremos la posterior predictiva para checar nuestro modelo.
Vamos a crear varias muestras, del mismo tamaño que la original, según nuestra
predictiva posterior, y compararemos estas muestras con la observada.

Y ahora simulamos otra muestra

```{r}
muestra_sim <- simular_normal_invgamma(20, pars_posterior)
muestra_sim %>% round(0)
```

Podemos simular varias muestras y hacer una prueba de lineup:

```{r}
library(nullabor)
sims_obs <- tibble(.n = 1:19) %>%
  mutate(estatura_cm = map(.n, ~ simular_normal_invgamma(20, pars_posterior))) %>%
  unnest(estatura_cm)
set.seed(9921)
pos <- sample(1:20, 1)
lineup_tbl <- lineup(true = cantantes %>% select(estatura_cm),
                     samples = sims_obs, pos = pos)
ggplot(lineup_tbl, aes(x = estatura_cm)) + geom_histogram(binwidth = 2.5) +
  facet_wrap(~.sample)

attr(lineup_tbl, "pos")
```
Con este tipo de gráficas podemos checar desajustes potenciales de nuestro modelo.

```{block, type="ejercicio"}

- ¿Puedes encontrar los datos verdaderos? ¿Cuántos seleccionaron los
datos correctos?
- Prueba hacer pruebas con una gráfica de cuantiles. ¿Qué problema
ves y cómo lo resolverías?

```


### Ejemplo: modelo Poisson {-}

Supongamos que pensamos el modelo para las observaciones es Poisson con
parámetro $\lambda$. Pondremos como inicial para $\lambda$ una exponencial con
media 10.

Nótese que la posterior está dada por

$$P(\lambda|x_1,\ldots, x_n) \propto e^{-n\lambda}\lambda^{\sum_i x_i} e^{-0.1\lambda} = \lambda^{n\overline{x}}e^{-\lambda(n + 0.1)}\,,$$
que es una distribución gamma con parámetros $(n\overline{x} + 1, n+0.1)$

Ahora supongamos que observamos la siguiente muestra, ajustamos nuestro modelo
y hacemos replicaciones posteriores de los datos observados:

```{r}
x <- rnbinom(250, mu = 20, size = 3)
crear_sim_rep <- function(x){
  n <- length(x)
  suma <- sum(x)
  sim_rep <- function(rep){
    lambda <- rgamma(1, sum(x) + 1, n + 0.1)
    x_rep <- rpois(n, lambda)
    tibble(rep = rep, x_rep = x_rep)
  }
}
sim_rep <- crear_sim_rep(x)
lineup_tbl <- map(1:5, ~ sim_rep(.x)) %>%
  bind_rows() %>%
  bind_rows(tibble(rep = 6, x_rep = x))
ggplot(lineup_tbl, aes(x = x_rep)) +
  geom_histogram(bins = 15) +
  facet_wrap(~rep)
```
Y vemos claramente que nuestro modelo no explica apropiadamente la variación
de los datos observados. Contrasta con:

```{r}
set.seed(223)
x <- rpois(250, 15)
crear_sim_rep <- function(x){
  n <- length(x)
  suma <- sum(x)
  sim_rep <- function(rep){
    lambda <- rgamma(1, sum(x) + 1, n + 0.1)
    x_rep <- rpois(n, lambda)
    tibble(rep = rep, x_rep = x_rep)
  }
}
sim_rep <- crear_sim_rep(x)
lineup_tbl <- map(1:5, ~ sim_rep(.x)) %>%
  bind_rows() %>%
  bind_rows(tibble(rep = 6, x_rep = x))
ggplot(lineup_tbl, aes(x = x_rep)) +
  geom_histogram(bins = 15) +
  facet_wrap(~rep)
```
Y verificamos que en este caso el ajuste del modelo es apropiado.


## Predicción {-}

Cuando queremos hacer predicciones particulares acerca de datos
que observemos en el futuro, también podemos usar la
posterior predictiva. En este caso, tenemos que considerar

1. La variabilidad que produce la incertidumbre en la estimación de los parámetros
2. La variabilidad de las observaciones dados los parámetros.

Es decir, tenemos que simular sobre todos las combinaciones factibles de los
parámetros.

### Ejemplo: cantantes {-}

Si un nuevo tenor llega a un coro, ¿cómo hacemos una predicción de su estatura? Como
siempre, quisiéramos obtener un intervalo que exprese nuestra incertidumbre acerca
del valor que vamos a observar. Entonces haríamos:

```{r}
sims_posterior <- sim_params(50000, pars_posterior) %>%
  mutate(y_pred = rnorm(n(), mu, sigma))
sims_posterior %>% head
```
```{r}
f <- c(0.025, 0.5, 0.975)
sims_posterior %>% summarise(f = f, y_pred = quantile(y_pred, f))
```

Y con esto obtenemos el intervalo (163, 189), al 95%, para una nueva
observación. Nótese que este intervalo no puede construirse con una simulación
particular de la posterior de los parámetros, pues sería demasiado corto.

Es posible demostrar que en este caso, la posterior predictiva tiene una forma
conocida:

- La posterior predictiva para el modelo normal-gamma inverso es una
distribución $t$ con $2\alpha'$ grados de libertad, centrada en $\mu'$, y con
escala $s^2 = \frac{\beta'}{\alpha'}\frac{n + n_0 + 1}{n +n_0}$

```{r}
mu_post <- pars_posterior[1]
n_post <- pars_posterior[2]
alpha_post <- pars_posterior[3]
beta_post <- pars_posterior[4]
s <- sqrt(beta_post/alpha_post) * sqrt((n_post + 1)/n_post)
qt(c(0.025, 0.5, 0.975), 2 * alpha_post) * s + mu_post
```




```{block, type='ejercicio'}

- Calcula la posterior predictiva del modelo Beta-Bernoulli y Beta-Binomial.
- (Más difícil) Calcula la posterior predictiva del modelo Poisson-Gamma.

```

### Ejemplo: posterior predictiva de Pareto-Uniforme. {-}

 La posterior predictiva del modelo Pareto-Uniforme no tiene un nombre estándar, pero
 podemos aproximarla usando simulación. Usando los mismos datos del ejercicio de la lotería, haríamos:

```{r}
rpareto <- function(n, theta_0, alpha){
  # usar el método de inverso de distribución acumulada
  u <- runif(n, 0, 1)
  theta_0 / (1 - u)^(1/alpha)
}
# Simulamos de la posterior de los parámetros
lim_inf_post <- max(c(300, muestra_loteria$numero))
k_posterior <- nrow(muestra_loteria) + 1.1
sims_pareto_posterior <- tibble(
  theta = rpareto(100000, lim_inf_post, k_posterior))
# Simulamos una observación para cada una de las anteriores:
sims_post_pred <- sims_pareto_posterior %>%
  mutate(x_pred = map_dbl(theta, ~ runif(1, 0, .x)))
# Graficamos
ggplot(sims_post_pred, aes(x = x_pred)) +
  geom_histogram(binwidth = 50) +
  geom_vline(xintercept = lim_inf_post, colour = "red")
```

Que es una mezcla de una uniforme con una Pareto.

## Inferencia estadística como un problema de decisión {-}

Hasta ahora hemos argumentado que el tipo de resultado de
un análisis de datos bayesiano considera la distribución posterior completa. Se
puede resumir esta distribución y reportar ciertas cantidades
puntuales (para realizar estimación puntual) o se pueden reportar
intervalos de probabilidad (para realizar estimación por intervalos). 
Ambas ideas también las hemos explorado bajo un enfoque frecuentista. 

Para **estimación puntual** consideremos que $\mathcal{D} = \mathcal{E} =
\Theta$.  Para la función de utilidad consideraremos $v(\hat \theta,
\theta)$ la pérdida de estimar $\hat \theta$ el verdadero parámetro de
interés $\theta$. Ahora consideraremos casos puntuales: 

1. **Pérdida cuadrática**:
$$v(\hat \theta,\theta) = (\hat \theta- \theta)^2\,.$$

En este caso la decisión óptima que minimiza la pérdida esperada es 
$$\hat \theta = \mathbb{E}(\theta)\,.$$
Es decir, la media de la distribución que contempla nuestro estado de 
conocimiento. 

2. **Pérdida en valor absoluto**:
$$v(\hat \theta,\theta) = |\hat \theta- \theta|\,.$$

En este caso la decisión óptima que minimiza la pérdida esperada es 
$$\hat \theta = \mathsf{mediana}(\theta)\,.$$
Es decir, la mediana de la distribución que contempla nuestro estado de 
conocimiento. 

3. **Pérdida vecindad**:
$$v(\hat \theta,\theta) = 1 - \mathbb{I}_{B_\epsilon(\hat \theta)}(\theta)\,,$$
denota la vecindad con radio $\epsilon$ alrededor de $\hat \theta$. 
En este caso, la decisión óptima que minimiza la pérdida esperada 
cuando $\epsilon \rightarrow 0$ es 

$$\hat \theta = \mathsf{moda}(\theta)\,.$$
Es decir, la moda media de la distribución que contempla nuestro estado de 
conocimiento. 

Para **estimación por intervalos** consideramos que $\mathcal{D} = \{
I : I \subset \Theta\}$, donde $I$ es un intervalo de probabilidad
$(1-\alpha)$ si 
$$ \int_{I}^{} P(\theta)d\theta = 1 - \alpha\,.$$

**Nota** que para una $\alpha$ fija no existe un intervalo único de
probabilidad. 

El espacio de eventos es el espacio parametral y la función de pérdida es
$$v(I, \theta) = \|D\| - \mathbb{I}_I(\theta)\,.$$ 
De lo anterior, observamos que para un nivel dado $\alpha$ es preferible
reportar el intervalo de probabilidad que sea de tamaño mínimo. 

Para **contraste de hipótesis** consideremos por simplicidad, dos hipótesis que
denotaremos por $H_0$ y $H_1$. Estas pueden representar dos modelos
distintos. En este caso (en el marco de decisión) el espacio de
alternativas está definido por
$$ \mathcal{D} = \{H_0, H_1\}\,.$$
El espacio parametral (espacio de eventos) está definido en términos
de estas posibilidades. Es decir,   
$$ \mathcal{E}= \{H_0, H_1\}\,.$$
La función de pérdida la definimos en términos de la decisión que tomamos
y lo que en realidad sucede. Esta función de pérdida la representamos como

\begin{center}
\begin{tabular}{l|l|l|}
$  v(d, \theta)$ & $H_0$     & $H_1$     \\
    \hline
$H_0$ & $v_{0,0}$ & $v_{0,1}$ \\
$H_1$ & $v_{1,0}$ & $v_{1,1}$\\
\hline
\end{tabular}
\end{center}

donde $v_{0,0}$, $v_{1,1}$ son la pérdida de tomar una decisión correcta 
(usualmente igual a 0). El costo $v_{1,0}$ es la pérdida de rechazar $H_0$
cuando $H_0$ es cierta y $v_{0,1}$ es la pérdida de no rechazar $H_0$ cuando 
$H_0$ es falsa. 

Se denotamos por $p_0 = P(H_0)$ la probabildiad asociada a la hipótesis $H_0$ 
al momento de tomar una decisión. Entonces la pérdida esperada para cada decisión
es 
\begin{gather}
\mathbb{E}[v(H_0)] = v_{0,0} p_0 + v_{0,1} (1 - p_0) = v_{0,1} - (v_{0,1}  - v_{0,0} ) p_0\,, \\
\mathbb{E}[v(H_1)] = v_{1,0} p_0 + v_{1,1} (1 - p_0) = v_{1,1} - (v_{1,1}  - v_{1,0} ) p_0\,. \\
\end{gather}

Que podemos pensar son dos rectas que intersectan en el punto 
\begin{align}
p^* = \frac{v_{0,1} - v_{1,1}}{v_{1,0}- v_{1,1} + v_{0,1} - v_{0,0}} \,.
\end{align}

De donde se sigue que si 
\begin{align}
\mathbb{E}[v(H_0)] < \mathbb{E}[v(H_1)] \Leftrightarrow \frac{p_0}{1 - p_0} < \frac{ v_{0,1} - v_{0,0}}{v_{1,0}- v_{1,1}} \Leftrightarrow p_0 > p^*\,,
\end{align}
que implica escoger $H_0$. Lo cual indica que si $p_0$ es
*suficientemente* grande comparada con $1 - p_0$ entonces decidiríamos
por $H_0$.

El caso contrario es análogo y podemos ver si $p_0$ es *suficientemente*
pequeña comparada con $1-p_0$ entonces decidimos por $H_1$. 

Desde el punto de vista bayesiano tenemos una tercera opción que no teníamos bajo 
el enfoque frecuentista de contraste de hipótesis. Si $p_0 = p^*$ entonces somos 
indiferentes y no podemos decir por una u otra alternativa. 

### Incorporando observaciones para contraste hipótesis {-}

Siguiendo con el marco de contraste de hipótesis bayesiana podemos considerar que
en las ecuaciones de arriba sólo hemos utilizado $P(H_0)$. Es decir, 
la probabilidad a priori de que ocurra $H_0$. Sin embargo, siguiendo un poco la 
intuición de nuestra discusión con las función predictiva, es natural considerar 
la distribución posterior de $H_0$ una vez observados los datos $P(H_0|X)$. 

Suponiendo que los costos de no equivocarnos son 0, entonces escogeríamos $H_0$ si 
y sólo si 
\begin{align}
\frac{P(H_0|X)}{P(H_1|X)} < \frac{v_{0,1}}{v_{1,0}}\,.\\
\end{align}
Lo cual es equivalente a 
\begin{align}
\frac{P(X|H_0)}{P(X|H_1)} < \frac{P(H_1)}{P(H_0)} \times \frac{v_{0,1}}{v_{1,0}}\,,\\
\end{align}

que nos indica que la decisión (o indecisión) se toma de acuerdo al
contraste del ajuste de un modelo sobre otro (cociente de
verosimilitudes) contra un factor que incorpora un balance entre los
costos asociados a la decisión incorrecta y los pesos relativos de
nuestras creencias iniciales.