tereom/fundamentos-2021
Notas del curso fundamentos de estadistica.

Package index
Search the tereom/fundamentos-2021 package
Vignettes
Functions
11
Source code
5
Man pages
0
Home
/
GitHub
/
tereom/fundamentos-2021
/

In tereom/fundamentos-2021: Notas del curso fundamentos de estadistica.

Introducción a inferencia bayesiana

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: 

```{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: 

## 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. 

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'}

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:

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
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:

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

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:

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))

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:

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:

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.

Observaciones:

$$\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:

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:

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:

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:

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:

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

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

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

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

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

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.

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:

Y más en general:

También aplicamos:

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:

```{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

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$:

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:

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.

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:

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:

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.

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:

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

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:

#### 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

muestra_sim <- simular_normal_invgamma(20, pars_posterior)
muestra_sim %>% round(0)

Podemos simular varias muestras y hacer una prueba de lineup:

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"}

### 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:

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:

sims_posterior <- sim_params(50000, pars_posterior) %>%
  mutate(y_pred = rnorm(n(), mu, sigma))
sims_posterior %>% head

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:

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'}

### 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.

  1. 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.

  1. 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.



tereom/fundamentos-2021 documentation built on Dec. 23, 2021, 8:46 a.m.

R Package Documentation

rdrr.io home R language documentation Run R code online

Browse R Packages

CRAN packages Bioconductor packages R-Forge packages GitHub packages

We want your feedback!

Note that we can't provide technical support on individual packages. You should contact the package authors for that.
GitHub issue tracker
ian@mutexlabs.com
Personal blog

 
Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.

Close