10.12 Reparametrización

El muestreador de Stan puede ser muy lento cuando la geometría de la distribución posterior es complicada, una manera de acelerar la convergencia es reparametrizando.

Ejemplo: El embudo de Neal

Veamos como ejemplo el embudo de Neal, tomado de la guía de usuarios de Stan, Neal (2003) define una distribución en la que ejemplifica las dificultades de muestrear de modelos jerárquicos, es un ejemplo extremo pero es fácil ver como la reparametrización hace que el muestreo de la posterior se simplifique.

En este ejemplo la densidad es,

\[p(y,x)=normal(y|0,3)*\prod_{n=1}^9normal(x_n|0,exp(y/2))\]

Y las curvas de nivel de probabilidad tienen la forma de embudos de 10 dimensiones. El cuello de los embudos es muy estrecho por la transformación de la variable \(y\).

Y podemos implementar en Stan:

Cuando el modelo está expresado de esta manera, Stan tiene problemas para muestrear del cuello del embudo, pues cuando \(y\) es chica \(x\) está restringida a valores cercanos a cero, esto se debe a que la escala de la densidad cambia con \(y\) de tal manera que el tamaño de paso que funciona bien en el cuerpo de la densidad será muy grande en eñ cuello y viceversa, un tamaño de paso que funciona bien en el cuello será ineficiente en el cuerpo.

En este ejemplo particular podemos reparametrizar y escribir el modelo de la siguiente forma que hace el muestreador más eficiente:

En este segundo modelo, x_raw y y_raw se muestrean de manera independiente de normales estándar, lo cuál es fácil para Stan. En un segundo paso se transforman a muestras del embudo.

no_funnel_cpp <- stan_model("src/stan_files/no_funnel.stan")
no_funnel_sims <- sampling(object = no_funnel_cpp, chains = 3 , iter = 1000)
no_funnel_sims
#> Inference for Stan model: no_funnel.
#> 3 chains, each with iter=1000; warmup=500; thin=1; 
#> post-warmup draws per chain=500, total post-warmup draws=1500.
#> 
#>           mean se_mean    sd   2.5%   25%   50%   75% 97.5% n_eff Rhat
#> y_raw     0.01    0.02  0.98  -1.96 -0.63  0.01  0.68  1.87  2626 1.00
#> x_raw[1] -0.01    0.02  1.04  -2.10 -0.77 -0.01  0.69  1.99  2440 1.00
#> x_raw[2]  0.00    0.02  0.98  -1.91 -0.67  0.01  0.66  1.85  2237 1.00
#> x_raw[3] -0.02    0.02  0.97  -1.86 -0.68 -0.05  0.65  1.96  2489 1.00
#> x_raw[4] -0.01    0.02  1.01  -1.91 -0.72 -0.01  0.70  1.99  3252 1.00
#> x_raw[5]  0.02    0.02  1.02  -1.93 -0.65  0.00  0.70  1.99  2682 1.00
#> x_raw[6]  0.00    0.02  1.01  -1.97 -0.68  0.00  0.70  1.91  3389 1.00
#> x_raw[7] -0.01    0.02  1.01  -1.97 -0.70  0.00  0.63  2.03  3093 1.00
#> x_raw[8]  0.01    0.02  0.99  -1.89 -0.66  0.02  0.69  1.95  2651 1.00
#> x_raw[9] -0.02    0.02  1.00  -2.01 -0.71  0.00  0.67  1.88  2605 1.00
#> y         0.02    0.06  2.94  -5.89 -1.90  0.03  2.03  5.60  2626 1.00
#> x[1]      0.09    0.33 11.08  -9.60 -0.67  0.00  0.66  8.63  1110 1.00
#> x[2]      0.16    0.18  6.58  -8.21 -0.55  0.00  0.62  9.59  1300 1.00
#> x[3]     -0.51    0.22  7.36 -10.37 -0.58 -0.01  0.51  7.33  1166 1.00
#> x[4]      0.01    0.17  6.03  -8.65 -0.65  0.00  0.57  9.26  1302 1.00
#> x[5]     -0.39    0.34 10.42 -10.30 -0.58  0.00  0.59  9.04   916 1.00
#> x[6]      0.00    0.19  6.42  -9.02 -0.58  0.00  0.64  8.82  1171 1.00
#> x[7]     -0.24    0.23  8.14 -11.23 -0.63  0.00  0.54  8.85  1283 1.00
#> x[8]     -0.09    0.22  8.17  -8.58 -0.57  0.00  0.58  9.51  1359 1.00
#> x[9]     -0.10    0.29 10.12  -9.66 -0.58  0.00  0.57 10.03  1190 1.00
#> lp__     -5.00    0.09  2.25 -10.44 -6.25 -4.66 -3.31 -1.75   617 1.01
#> 
#> Samples were drawn using NUTS(diag_e) at Tue Dec 10 23:29:28 2019.
#> For each parameter, n_eff is a crude measure of effective sample size,
#> and Rhat is the potential scale reduction factor on split chains (at 
#> convergence, Rhat=1).

Modelos jerárquicos y parametrización no centrada

En modelos jerárquicos es común encontrarse con geometrías complejas, una parametrización que suele ayudar para acelerar convergencia es la parametrización no centrada.

Parametrización centrada

Un modelo jerárquico usual seleccionaría muestras de un vector de coeficientes \(\beta\) como sigue:

parameters {
  real mu_beta;
  real<lower=0> sigma_beta;
  vector[K] beta;
  ...
model {
  beta ~ normal(mu_beta, sigma_beta);
  ...

Y tendría ineficiencias como el embudo de Neal, debido a que el valor de \(\beta\) (vector), \(\mu_{\beta}\) y \(\sigma_{\beta}\) tienen alta correlación en la posterior. El nivel de correlación dependerá de la cantidad de datos disponibles siendo mayor cuanto menos información se tenga. En estos casos, de datos limitados, es más eficiente usar una parametrización no centrada.

Parametrización no centrada

parameters {
  vector[K] beta_raw;
  ...
transformed parameters {
  vector[K] beta;
  // implies: beta ~ normal(mu_beta, sigma_beta)
  beta = mu_beta + sigma_beta * beta_raw;
model {
  beta_raw ~ std_normal();
  ...

Las iniciales de mu_beta y sigma_beta permanecen sin cambios respecto al modelo original.

Ahora veremos un ejemplo práctivo de parametrización no centrada en un modelo de predicción de resultados del conteo rápido, que además ejemplifica un flujo de trabajo bayesiano, con pasos de ajuste, inferencia, calibración y evaluación de robustez. El ejemplo está en esta liga.