Байесовский анализ сложных моделей возможен с помощью алгоритмов Монте-Карло по схеме Марковской цепи (MCMC). По идее, исследователь должен не только сформулировать модель, априорное мнение о параметрах модели, но и самостоятельно сконструировать цепь сходящуюся к апостериорному распределению параметров. Однако оказалось, что конструирование цепи в большинстве ситуаций можно доверить компьютеру.
Программа STAN — это не пакет R, а отдельный язык описания байесовских моделей. Исследователь описывает модель и априорное распределение параметров, а STAN генерирует код C++, в котором запрограммирована требуемая марковская цепь. Среди альтернатив программе STAN можно назвать JAGS, R + пакет NIMBLE, python + пакет PyMC, Julia + змея Mamba и Julia + Klara.
Для работы требуется установить пакет rstan, который содержит и свежую версию программы STAN. В силу того, что здесь взаимодействует много программ (R, STAN, компилятор C++), установку следует выполнять аккуратно, строго по инструкции на сайте STAN.
library("rstan") # врубим на полную мощь Байеса!
library("shinystan") # и на полную громкость интерактивный анализ полученной цепиПара опций: определяем число ядер процессора, чтобы при симуляциях использовались все ядра и просим STAN сохранить скомпилированную модель, чтобы не перекомпилировать её для каждой цепи.
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())Модель 1: Байесовский МНК
model_bayesian_ls <- "
data {
int n;
int k;
vector[n] y;
matrix[n, k] X;
}
parameters {
vector[k] beta;
real<lower=0> sigma;
}
model {
y ~ normal(X * beta, sigma);
}
"Блок data описывает, какие данные нужно передать STAN. Блок model описывает, хм, неужели модель?
\[ y|\beta, \sigma, X \sim \mathcal{N}(X\beta, \sigma^2 \cdot I) \]
Блок parameters описывает априорное распределение. Если не указывать явно, то по умолчанию предполагаются равномерные распределения на допустимом множестве. Для вектора бет это означает: \[
\beta_j \sim U(-\infty, \infty)
\] Для константы \(\sigma\) мы подразумеваем: \[
\sigma \sim U(0, \infty)
\]
Конечно, настоящего равномерного распределения на неограниченном интервале не бывает, это несобственное распределение. Однако в данном случае апостериорное распределение оказывается собственным.
Можно явно указать любое другое априорное распределение.
Формируем данные:
n <- nrow(cars)
k <- 2
y <- cars$dist
X <- cbind(rep(1, n), cars$speed)
cars_data <- list(n = n, k = k, y = y, X = X)Запускаем 4 цепи:
fit <- stan(model_code = model_bayesian_ls, data = cars_data,
iter = 1000, chains = 4)Апостериорное распределение параметров:
fit## Inference for Stan model: dd48765ea01c4aabdbe455473e42eb1a.
## 4 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=2000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff
## beta[1] -17.53 0.26 6.84 -31.26 -21.94 -17.59 -13.18 -4.04 679
## beta[2] 3.93 0.02 0.42 3.07 3.67 3.94 4.20 4.77 681
## sigma 15.75 0.05 1.62 12.93 14.68 15.60 16.65 19.47 981
## lp__ -159.37 0.05 1.24 -162.68 -159.92 -159.00 -158.49 -158.01 746
## Rhat
## beta[1] 1
## beta[2] 1
## sigma 1
## lp__ 1
##
## Samples were drawn using NUTS(diag_e) at Wed Sep 28 08:59:18 2016.
## 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).Модель 2: Байесовский логит
Модель 3:
Описываем модель на языке STAN:
model_8_schools <- "
data {
int<lower=0> J; // number of schools
real y[J]; // estimated treatment effects
real<lower=0> sigma[J]; // s.e. of effect estimates
}
parameters {
real mu;
real<lower=0> tau;
real eta[J];
}
transformed parameters {
real theta[J];
for (j in 1:J)
theta[j] = mu + tau * eta[j];
}
model {
eta ~ normal(0, 1);
y ~ normal(theta, sigma);
}
"Задаем данные:
schools_data <- list(J = 8,
y = c(28, 8, -3, 7, -1, 1, 18, 12),
sigma = c(15, 10, 16, 11, 9, 11, 10, 18))Запускаем 4 цепи:
fit <- stan(model_code = model_8_schools, data = schools_data,
iter = 1000, chains = 4)Апостериорное распределение параметров:
fit## Inference for Stan model: 4cd96c47cbb67b368bab02377102f603.
## 4 chains, each with iter=1000; warmup=500; thin=1;
## post-warmup draws per chain=500, total post-warmup draws=2000.
##
## mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
## mu 8.08 0.15 4.74 -1.08 5.05 7.86 11.10 17.97 948 1.00
## tau 6.33 0.22 5.16 0.24 2.31 4.99 9.15 18.70 529 1.01
## eta[1] 0.35 0.03 0.93 -1.49 -0.26 0.38 1.02 2.08 1391 1.00
## eta[2] -0.03 0.02 0.84 -1.64 -0.59 -0.03 0.50 1.65 1529 1.00
## eta[3] -0.17 0.03 0.90 -1.93 -0.77 -0.19 0.41 1.68 1262 1.00
## eta[4] -0.04 0.02 0.90 -1.85 -0.62 -0.04 0.54 1.74 1535 1.00
## eta[5] -0.35 0.02 0.82 -1.92 -0.92 -0.37 0.18 1.31 1366 1.00
## eta[6] -0.24 0.02 0.88 -1.95 -0.83 -0.25 0.35 1.48 1406 1.00
## eta[7] 0.30 0.03 0.88 -1.45 -0.28 0.32 0.89 2.06 1173 1.00
## eta[8] 0.02 0.03 0.89 -1.75 -0.56 0.02 0.63 1.75 1024 1.00
## theta[1] 11.53 0.30 8.47 -1.10 6.03 9.82 15.60 33.39 800 1.00
## theta[2] 7.80 0.14 5.94 -3.87 4.17 7.88 11.43 19.81 1793 1.00
## theta[3] 6.33 0.20 7.35 -11.38 2.72 6.93 10.82 19.96 1311 1.00
## theta[4] 7.65 0.15 6.30 -4.69 3.75 7.49 11.62 20.56 1781 1.00
## theta[5] 5.31 0.14 6.04 -7.94 1.58 5.62 9.22 16.75 1808 1.00
## theta[6] 5.99 0.17 6.67 -9.28 2.42 6.43 10.19 18.19 1545 1.00
## theta[7] 10.62 0.20 7.04 -1.73 6.02 9.87 14.62 27.32 1228 1.00
## theta[8] 8.37 0.23 7.70 -6.49 3.67 8.04 12.50 25.17 1148 1.00
## lp__ -4.68 0.10 2.54 -9.79 -6.24 -4.44 -2.87 -0.30 594 1.00
##
## Samples were drawn using NUTS(diag_e) at Wed Sep 28 09:00:06 2016.
## 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).Графики:
….
Тесты на сходимость цепей:
….
Документация STAN содержит кучу встроенных примеров:
stan_demo()Визуализировать полученную цепь КРАЙНЕ удобно интерактивно с помощью команды
shiny_stan_object <- launch_shinystan(fit)Почиташки: