Distribución Normal

Retomemos el caso de observaciones normales, supongamos que tengo una muestra \({x}_{1},...,{x}_{N}\) de observaciones independientes e identicamente distribuidas, con \({x}_{i}∼N({\mu},{\sigma}^{2})\), veremos el caso de media desconocida, varianza desconocida y de ambas desconocidas.

Normal con media desconocida. Supongamos que \({\sigma}^{2}\) es conocida, por lo que nuestro parámetro de interés es únicamente \({\mu}\) entonces si describo mi conocimiento inicial de \({\mu}\) a través de una distribución normal:

\[{\mu}{\sim}{N}(m,{\tau}^{2})\]

resulta en una distribución posterior:

\[ {\mu}|{x}∼{N}\left(\frac{σ2}{{\sigma}^{2}+N{\tau}^{2}}m+\frac{N{\sigma}^{2}}{{\sigma}^{2}+N{\tau}^{2}}\bar{x},\frac{{\sigma}^{2}{\tau}^{2}}{{\sigma}^{2}+N{\tau}^{2}}\right) \]

Normal con varianza desconocida. Supongamos que μ es conocida, por lo que nuestro parámetro de interés es únicamente \({\sigma}^{2}\). En este caso una distribución conveniente para describir nuestro conocimiento inicial es la distribución Gamma Inversa.

La distribución Gamma Inversa es una distribución continua con dos parámetros y que toma valores en los positivos. Como su nombre lo indica, esta distribución corresponde al recírpoco de una variable cuya distribución es Gamma, recordemos que si \({x}{\sim}{Gamma}(\alpha,\beta)\) entonces:

\[ p(x)=\frac{1}{{\beta}^{\alpha}{\Gamma}({\alpha})}{x}^{\alpha−1}{e}^{−\frac{x}{\beta}} \]

donde \(x>0\). Ahora si \(y\) es la variable aleatoria recírpoco de \(x\) entonces:

\[ p(y)=\frac{{\beta}^{\alpha}}{{\Gamma}({\alpha})}{y}^{-\alpha−1}{e}^{−\frac{\beta}{y}} \]

con media

\[ \frac{\beta}{\alpha-1} \]

y varianza

\[ \frac{{\beta}^{2}}{({\alpha-1})^{2}({\alpha-2})} \]

Debido a la relación entre las distribuciones Gamma y Gamma Inversa, podemos utilizar la función rgamma de R para generar valores con distribución gamma inversa.

x_gamma <- rgamma(2000, shape = 5, rate = 1)

x_igamma <- 1 / x_gamma


library(MCMCpack)

## Loading required package: coda

## Loading required package: MASS

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:dplyr':
## 
##     select

## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)

## ## Copyright (C) 2003-2023 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park

## ##
## ## Support provided by the U.S. National Science Foundation

## ## (Grants SES-0350646 and SES-0350613)
## ##

x_igamma <- data.frame(x_igamma)

ggplot(x_igamma, aes(x = x_igamma)) +
  geom_histogram(aes(y = ..density..), binwidth = 0.05, fill = "gray") + 
  stat_function(fun = dinvgamma, args = list(shape = 5, scale = 1), 
    color = "darkgreen")  

## Warning: The dot-dot notation (`..density..`) was deprecated in ggplot2 3.4.0.
## ℹ Please use `after_stat(density)` instead.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

Volviendo al problema de inferir acerca del parámetros \({\sigma}^{2}\), si resumimos nuestro conocimiento inicial a través de una distribución Gamma Inversa tenemos

\[ p({\sigma}^{2})=\frac{{\beta}^{\alpha}}{{\Gamma}({\alpha})}\frac{1}{{y}^{\alpha+1}}{e}^{−\frac{\beta}{{\sigma}^{2}}} \]

la verosimiltud:

\[ p({x}|{\mu},{\sigma}^{2})=\frac{1}{({2{\pi}{\sigma}^{2}})^{\frac{N}{2}}}exp\left(−\frac{1}{2{\sigma}^{2}}{\sum}_{j=1}^{N}({x}_{j}−{\mu})2\right) \]

y calculamos la posterior:

\[ p(σ2)∝p(x|{\mu},{\sigma}^{2})p({\sigma}^{2}) \]

obtenemos que \({{\sigma}^{2}|{x}}{\sim}{GI}\left[{\frac{N}{2}+{\alpha},{\beta}+\frac{1}{2}{\sum}({{x}_{i}−{\mu}})^{2}}\right]\)

.

Por tanto tenemos que la inicial Gamma con verosimilitud Normal es una familia conjugada.

Normal con media y varianza desconocidas

Sea \({\theta}=({\mu},{\sigma}^{2})\) especificamos la siguiente inicial para \(\theta\):

\[p({\theta})=N({\mu}|{m},{\tau}^{2}){IG}({\sigma}^{2}|{\alpha},{\beta})\]

suponemos hiperparámetros \(m\),\({\tau}^{2}\),\(\alpha\),\(\beta\) conocidos. Entonces, la distribución posterior es:

\[p({\theta}|{x})∝p(x|{\theta})p({\theta})=\frac{1}{{({\sigma}^{2}})^{\frac{N}{2}}}exp\left[−\frac{1}{2{\sigma}^{2}}{\sum}_{i=1}^{N}({{x}_{i}−{\mu}})^{2}\right]exp\left[−\frac{1}{2{\tau}^{2}}({{\mu}−{m}})^{2}\right]\frac{1}{({\sigma}^{2})^{{\alpha}+1}}exp\left(−\frac{{\beta}}{{\sigma}^{2}}\right)\]

en esta última distribución no reconocemos el núcleo de niguna distribución conocida (existe una distribución normal-gamma inversa) pero si nos concenteramos únicamente en los términos que involucran a \({\mu}\) tenemos:

\[exp\left\{−\frac{1}{2}\left[{\mu}^{2}\left(\frac{N}{{\sigma}^{2}}+\frac{1}{{\tau}^{2}}\right)−2{\mu}\left(\frac{{{\sum}_{i=1}^{n}}{{x}_{i}}}{σ2}+\frac{m}{{\tau}^{2}}\right)\right]\right\}\]

esta expresión depende de \({\mu}\) y \({\sigma}^{2}\), sin embargo condicional a \({\sigma}^{2}\) observamos el núcleo de una distribución normal,

\[{{\mu}|{\sigma}^{2},{x}}{\sim}N\left(\frac{{n}{\tau}^{2}}{{n}{\tau}^{2}+{\sigma}^{2}}\bar{x}+\frac{{\sigma}^{2}}{{N}{\tau}^{2}+{\sigma}^{2}}m,\frac{{\tau}^{2}{\sigma}^{2}}{{n}{\tau}^{2}+{\sigma}^{2}}\right)\]

Si nos fijamos únicamente en los tárminos que involucran a \({\sigma}^{2}\) tenemos:

\[\frac{1}{{({\sigma}^{2}})^{\frac{N}{2}+{\alpha}+{1}}}exp\left[−\frac{1}{{\sigma}^{2}}{{\sum}_{i=1}^{N}}{\frac{({{x}_{i}−{\mu}})^{2}}{2}}+{\beta}\right]\]

y tenemos

\[{{\sigma}^{2}|{\mu},{x}}{\sim}{GI}\left[{\frac{N}{2}+{\alpha},{\beta}+\frac{1}{2}{\sum}({{x}_{i}−{\mu}})^{2}}\right]\]

Obtenemos así las densidades condicionales completas \(p({{\mu}|{\sigma}^{2},{x}})\) y \(p({{\sigma}^{2}|{\mu},{x}})\) cuyas distribuciones conocemos y de las cuales podemos simular.

Implementaremos un muestreador de Gibbs.

Comenzamos definiendo las distrbuciones iniciales:

• \({\mu}{\sim}N(1.5,16)\), esto es \(m=1.5\) y \({\tau}^{2}=16\)

• \({\sigma}^{2}{\sim}GI(3,3)\) esto es \(\alpha=\beta=3\).

Ahora supongamos que observamos 20 realizaciones provenientes de la distribución de interés:

N <- 50
set.seed(122)
x <- rnorm(N, 2, 2) 
x

##  [1]  4.62140176  0.24829384  2.39904749  2.93190885 -1.60411353  4.89748691
##  [7]  2.59770769  2.72362329 -0.01388084  1.48600171  1.73574244  0.31673052
## [13]  2.54850449 -2.92518071 -2.30679198  4.31835150  3.37948021  3.76050265
## [19]  0.11325957  3.43814572  0.92434912  0.95470268 -0.10584378  2.20303449
## [25]  5.72700211  1.96078181 -0.15661507  2.34520855  3.06610819  5.90452895
## [31]  4.82270939  3.20273052  0.17200480  5.16085186  1.06016874  5.20367778
## [37]  2.74547989  2.06775571  2.20820677 -2.03674850  0.68398061 -1.21700489
## [43] -0.68818260  1.60665220  4.75132535  2.34663000  1.12144484 -1.19064581
## [49]  0.51144488  5.63041621

Escribimos el muestreador de Gibbs.

m <- 1.5; tau2 <- 16; alpha <- 3; beta <- 3

pasos <- 20000
camino <- matrix(0, nrow = pasos + 1, ncol = 2)
camino[1, 1] <- 0

# Generamos la caminata aleatoria
for (j in 2:(pasos + 1)){
  # sigma^2
  mu <- camino[j - 1, 1]
  a <- N / 2 + alpha
  b <- sum((x  - mu) ^ 2) / 2 + beta
  camino[j - 1, 2] <- 1/rgamma(1, shape = a, rate = b)
  
  # mu
  sigma2 <- camino[j - 1, 2]
  media <- (N * tau2 * mean(x) + sigma2 * m) / (N * tau2 + sigma2)
  var <- sigma2 * tau2 / (N * tau2 + sigma2)
  camino[j, 1] <- rnorm(1, media, sd = sqrt(var))
}

caminata <- data.frame(pasos = 1:pasos, mu = camino[1:pasos, 1], 
  sigma2 = camino[1:pasos, 2])

caminata_g <- caminata %>%
  gather(parametro, val, mu, sigma2) %>%
  arrange(pasos)

ggplot(filter(caminata_g, pasos > 15000), aes(x = pasos, y = val)) +
  geom_path(alpha = 0.3) +
  facet_wrap(~parametro, ncol = 1, scales = "free") +
  scale_y_continuous("")

ggplot(filter(caminata_g, pasos > 5000), aes(x = val)) +
  geom_histogram(fill = "gray") +
  facet_wrap(~parametro, ncol = 1, scales = "free") 

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Algunos resúmenes de la posterior:

caminata_g %>%
  filter(pasos > 1000) %>% # eliminamos la etapa de calentamiento
  group_by(parametro) %>%
  summarise(
    mean(val), 
    sd(val), 
    median(val)
    )

## # A tibble: 2 × 4
##   parametro `mean(val)` `sd(val)` `median(val)`
##   <chr>           <dbl>     <dbl>         <dbl>
## 1 mu               1.91     0.305          1.91
## 2 sigma2           4.74     0.933          4.62

Predicción. Para predecir el valor de una realización futura \(y\) recordemos que:

\[p(y)={\int}{p(y|\theta)p(\theta|x)}{d\theta}\]

Por tanto podemos aproximar la distribución predictiva posterior como:

caminata_f <- filter(caminata, pasos > 5000)

caminata_f$y_sims <- rnorm(1:nrow(caminata_f), caminata_f$mu, caminata_f$sigma2)

ggplot(caminata_f, aes(x = y_sims)) + 
  geom_histogram(fill = "gray") +
  geom_vline(aes(xintercept = mean(y_sims)), color = "darkgreen")

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

En estadística bayesiana es común parametrizar la distribución Normal en términos de precisión (el inverso de la varianza). Si parametrizamos de esta manera \({\nu}=\frac{1}{{\sigma}^{2}}\) podemos repetir el proceso anterior con la diferencia de utilizar la distribución Gamma en lugar de Gamma inversa.