Обзор библиотеки Stan в R
Приветствую!
Stan — это библиотека на C++, предназначенная для байесовского моделирования и вывода. Она использует сэмплер NUTS, чтобы создавать апостериорные симуляции модели, основываясь на заданных пользователем моделях и данных. Так же Stan может использовать алгоритм оптимизации LBFGS для максимизации целевой функции, к примеру как логарифмическое правдоподобие.
Для облегчения работы с Stan из языка программирования R доступен пакет rstan, который предоставляет интерфейс R для Stan.
Сегодня мы и рассмотрим этот пакет.
Установка
Открываем RStudio и ющаем команду в консоли для установки пакета:
нinstall.packages("rstan", dependencies = TRUE)Stan использует C++ для высокопроизводительных вычислений, поэтому нужно убедиться, что система настроена для компиляции C++ кода. Для этого может потребоваться установить Rtools для винды или Xcode для мака. После установки, выполняем команду:
pkgbuild::has_build_tools(debug = TRUE)Ну и на всякий случай чекнем, все ли ОК, это делается таким образом:
example(stan_model, package = "rstan", run.dontrun = TRUE)Если нет ошибок, то вы большой молодец.
Пакет rstan также сильно зависит от нескольких других пакетов R:
StanHeaders (заголовки Stan C++)
BH (заголовки Boost C++)
RcppEigen (собственные заголовки C++)
Rcpp (облегчает использование C++ из R)
встроенный (компилирует C++ для использования с R)
Эти зависимости должны быть установлены автоматически, если вы устанавливаете пакет rstan с помощью одного из обычных механизмов.
Немного про сам Stan
Stan представляет собой предметно-ориентированный язык, оптимизированный под задачи статистического моделирования
Stan использует декларативный синтаксис, т.е, юзер описывает модель в терминах статистических распределений без необходимости указывать алгоритмы для вычисления выводов. В этом есть только одни плюсы, и как минимум так — понятней и удобней.
Автоматическое дифференцирование позволяет Stan автоматически вычислять градиенты функции потерь, что мастхев для алгоритмов оптимизации и Монте-Карло методов марковских цепей, используемых для вывода.
Stan дает доступ к ряду алгоритмов для байесовского вывода, включая:
NUTS, который является расширением алгоритма HMC.
Вариационные методы, которые предлагают более быструю, но приближенную альтернативу MCMC для вывода.
Stan поддерживает модульность, позволяя юзерам определять собственные функции, что облегчает повторное использование кода и упрощает создание сложных моделей.
Stan имеет множество интеграций с различными ЯПами: R, Python, Julia и MATLAB.
Так же есть очень хорошее и активное сообщество на их форуме.
Структура программы Stan
Программа на Stan состоит из нескольких блоков, каждый из которых выполняет определенную роль в процессе моделирования:
Блок данных (
data): тут объявляются все входные данные, которые используются в модели. Это могут быть скаляры, векторы, матрицы и так далее. Данные в этом блоке считаются фиксированными и неизменными на протяжении всего процесса моделирования.Блок параметров (
parameters): здесь объявляются раметры модели, которые будут оцениваться на основе данных. Это могут быть, например, коэффициенты регрессии, дисперсии, масштабные коэффициенты и так далее. Stan автоматически применяет методы Монте-Карло с цепями Маркова для оценки этих параметров.Блок модели (
model): база программы на Stan. Здесь определяется вероятностная модель данных и параметров. Можно юзать встроенные распределения и функции Stan для описания априорных распределений параметров и вероятностей наблюдаемых данных при заданных параметрах.Блоки преобразований (
transformed data,transformed parameters): блоки используются для преобразования данных и параметров соответственно. Преобразованные данные могут быть юзабельны для упрощения модели или для предварительной обработки данных, а преобразованные параметры часто используются для выведения вторичных величин, которые вычисляются на основе оцененных параметров.Блок генерации (
generated quantities): тут можно вычислять величины, которые зависят от параметров после их оценки. Часто используются для генерации прогнозов и для вычисления статистик постериорного распределения.
Перейдем сразу к созданию, к примеру, модели линейной регрессии:
// определяем модель линейной регрессии
data {
int N; // количество наблюдений
vector[N] y; // зависимая переменная
int K; // колво предикторов
matrix[N, K] X; // матрица предикторов
}
parameters {
vector[K] beta; // коэф. регрессии
real alpha; // перехват (интерсепт)
real sigma; // отклонение ошибок
}
model {
// априорные распределения
beta ~ normal(0, 10); // нормальное априорное распределение для коэффициентов
alpha ~ normal(0, 10); // нормальное априорное распределение для перехвата
sigma ~ cauchy(0, 5); // распределение Коши для стандартного отклонения ошибок
// вероятностная модель данных
y ~ normal(X * beta + alpha, sigma); // нормальное распределение зависимой переменной
}
generated quantities {
vector[N] y_pred; // вектор предсказанных значени
for (i in 1:N) {
y_pred[i] = normal_rng(X[i] * beta + alpha, sigma); // генерация предсказаний
}
} Итак:
В блоке данных (data): объявляем необходимые данные для модели: количество наблюдений N, зависимую переменную y, количество предикторов K, и матрицу предикторов X.
Блок параметров (parameters): объявляются параметры модели, которые будут оцениваться: вектор коэффициентов beta, перехват alpha, и стандартное отклонение ошибок sigma.
Блок модели (model): определяются априорные распределения для параметров и вероятностная модель данных. Используем нормальное распределение для коэффициентов и перехвата, а распределение Коши для стандартного отклонения ошибок. Затем указываем, что наблюдаемые данные y распределены нормально с математическим ожиданием, равным линейной комбинации предикторов, и стандартным отклонением sigma.
Блок генерации (generated quantities): генерируем предсказанные значения y_pred для каждого наблюдения. Это делается с использованием функции normal_rng, которая генерирует случайные значения из нормального распределения с заданными параметрами.
Можно расширирять функциональность Stan с помощью пользовательских функций, написанных на R.
Допустим, у нас есть набор данных, и мы хотим использовать некоторые фичи из пакета R для предварительной обработки этих данных перед их использованием в модели Stan. Мы можем сделать это, выполнив предварительную обработку в R, а затем передав обработанные данные в модель Stan:
library(rstan)
library(dplyr)
# наш дф
df <- data.frame(
x = rnorm(100),
y = rnorm(100, 2, 4)
)
# dplyr для предварительной обработки данных
processed_data <- df %>%
mutate(x_scaled = scale(x), y_scaled = scale(y))
# заготовка данных для Stan
stan_data <- list(
N = nrow(processed_data),
x = processed_data$x_scaled,
y = processed_data$y_scaled
)
# stan
stan_code <- "
data {
int N;
vector[N] x;
vector[N] y;
}
parameters {
real alpha;
real beta;
real sigma;
}
model {
y ~ normal(alpha + beta * x, sigma);
}
"
# запуск модели
fit <- stan(model_code = stan_code, data = stan_data) Использование функций C++
Можно написать функции на C++, скомпилировать их в динамически подключаемые библиотеки (DLL или .so файлы в зависимости от вашей ос), а затем вызывать эти функции непосредственно из модели Stan.
Допустим есть следующая C++ функция, которую хотим использовать в модели Stan. Важно, чтобы функция была «чистой», т.е. не имела побочных эффектов, и ее вывод был полностью определен ее входными данными:
#include
extern "C" {
double my_custom_function(double x) {
return std::exp(x); // простой пример, возвращающий экспоненту входного значения
}
} Далее юзаем компилятор C++, к примеру g++, для компиляции функции в DLL или .so файл:
g++ -shared -fPIC -o my_custom_function.so my_custom_function.cppЧтобы использовать внешнюю функцию в модели Stan, нужно объявить ее в блоке functions:
functions {
// объявление внешней функции
real my_custom_function(real x);
}
model {
// использование функции в модели
real y = my_custom_function(0.5);
}При компиляции и запуска модели Stan из R, надо убедиться, что динамически подключаемая библиотека доступна для сеанса R. Юзаем функцию dyn.load в R для загрузки .so или DLL файла перед компиляцией модели Stan:
# загрузка динамически подключаемой библиотеки в R
dyn.load("path/to/my_custom_function.so")Модели временных рядов
Авторегрессионная модель первого порядка (AR (1)) предполагает, что текущее значение временного ряда зависит от его предыдущего значения с добавлением случайного шума:
data {
int N; // колво наблюдений
vector[N] y; // временной ряд
}
parameters {
real alpha; // параметр авторегрессии
real beta; // коэффициент уровня
real sigma; // стандартное отклонение шума
}
model {
for (n in 2:N)
y[n] ~ normal(alpha + beta * y[n-1], sigma);
} Запускаем:
library(rstan)
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# грузим данные
N <- 100
y <- rnorm(N, 0, 1) # данные временного ряда
# данные для Stan
stan_data <- list(N = N, y = y)
# комп. и запуск модели
fit_ar1 <- stan(file = 'ar1_model.stan', data = stan_data, iter = 2000, chains = 4)
print(fit_ar1)Модель с локальным уровнем и сезонностью подходит для анализа временных рядов с трендом и сезонными колебаниями:
data {
int N; // колво наблюдений
int K; // период сезонности
vector[N] y; // временной ряд
}
parameters {
vector[N] level; // лок. уровень
vector[K] season; // сезонные эффекты
real sigma_level; // шум уровня
real sigma_season; // шум сезонности
}
model {
for (n in 2:N)
level[n] ~ normal(level[n-1], sigma_level); // мдель уровня
for (n in (K+1):N)
y[n] ~ normal(level[n] + season[n mod K + 1], sigma_season); // сезонная модель
// приоры для сезонности
season ~ normal(0, sigma_season);
} Запускаем:
# допустим, что данные и файл модели уже подготовлены
K <- 12 # период сезонности
stan_data <- list(N = N, K = K, y = y)
fit_seasonal <- stan(file = 'seasonal_model.stan', data = stan_data, iter = 2000, chains = 4)
print(fit_seasonal)Модель с локальным трендом позволяет анализировать временные ряды с изменяющимся во времени трендом.
data {
int N;
vector[N] y;
}
parameters {
vector[N] mu; // локальный тренд
real sigma_mu; // шум тренда
}
model {
for (n in 2:N)
mu[n] ~ normal(mu[n-1], sigma_mu);
y ~ normal(mu, sigma_mu);
} Запускаем:
# также допускаем, что данные и файл модели уже подготовлены
stan_data <- list(N = N, y = y)
fit_trend <- stan(file = 'local_trend_model.stan', data = stan_data, iter = 2000, chains = 4)
print(fit_trend)Обработка отсутствующих данных и частично известных параметров
Моделирование отсутствующих данных как параметров. В этом подходе отсутствующие данные рассматриваются как неизвестные параметры, которые нужно оценить в процессе моделирования:
Рассмотрим случай, когда есть набор данных с отсутствующими значениями в переменной y:
data {
int N_obs; // колво наблюдаемых данных
int N_miss; // колво отсутствующих данных
real y_obs[N_obs]; // наблюд. данные
}
parameters {
real mu; // ср. значение распределения
real sigma; // дефолтное отклонение распределения
real y_miss[N_miss]; // отсутствующие данные, моделируемые как параметры
}
model {
y_obs ~ normal(mu, sigma); // модель для наблюдаемых данных
y_miss ~ normal(mu, sigma); // модель для отсутствующих данных
} y_miss моделируется так же, как и наблюдаемые данные, с теми же параметрами mu и sigma.
Использование предиктивных величин для обработки отсутствующих данных. Этот подход заключается в использовании предиктивных величин для заполнения отсутствующих данных:
data {
int N; // общее количество данных (наблюдаемые + отсутствующие)
int missing[N]; // индикатор отсутствующих данных
real y[N]; // данные, где отсутствующие значения обозначены как NaN
}
parameters {
real mu;
real sigma;
real y_miss[N];
}
model {
for (n in 1:N) {
if (missing[n] == 1) {
y_miss[n] ~ normal(mu, sigma); // модель для отсутствующих данных
} else {
y[n] ~ normal(mu, sigma); // модель для наблюдаемых данных
}
}
} missing является индикатором массива, который определяет, является ли значение в y отсутствующим. Если значение отсутствует, то для соответствующего элемента в y_miss применяется модель.
Иногда параметры модели могут быть частично известны из предыдущих исследований или экспертных оценок. Stan позволяет интегрировать эту информацию в модель:
data {
int N;
real y[N];
real sigma_prior; // частично известное стандартное отклонение
}
parameters {
real mu;
real sigma;
}
model {
sigma ~ normal(0, sigma_prior); // использование частично известного параметра
y ~ normal(mu, sigma);
} sigma_prior используется для задания априорного распределения параметра sigma
Модели выживаемости
Экспоненциальная модель предполагает, что интенсивность события (или риск) является постоянной во времени. Это самая простая форма модели выживаемости:
data {
int N; // колво наблюдений
vector[N] T; // тайм до события или цензурирования
int delta[N]; // индикатор события (1 если событие произошло, 0 если данные цензурированы)
}
parameters {
real lambda; // параметр интенсивности события
}
model {
// априорное распределение для lambda
lambda ~ gamma(0.01, 0.01);
// правдоподобие
for (n in 1:N) {
if (delta[n] == 1) {
// для наблюдений с событием
target += log(lambda) - lambda * T[n];
} else {
// для цензурированных наблюдений
target += -lambda * T[n];
}
}
} Модель Вейбулла расширяет экспоненциальную модель, позволяя риску изменяться со временем. Можно реализовать за счет добавления параметра формы:
data {
int N;
vector[N] T;
int delta[N];
}
parameters {
real lambda; // параметр масштаба
real k; // параметр формы
}
model {
lambda ~ gamma(0.01, 0.01);
k ~ gamma(0.01, 0.01);
for (n in 1:N) {
if (delta[n] == 1) {
target += log(k) + (k - 1) * log(T[n]) - pow(T[n] * lambda, k);
} else {
target += -pow(T[n] * lambda, k);
}
}
} Модель пропорциональных рисков Кокса позволяет анализировать влияние объясняющих переменных на риск без необходимости специфицировать базовую функцию риска.
В Stan прямое моделирование модели Кокса может быть сложным, и я думаю, это даже в какой-то мере заслуживает отдельной статьи из-за её частичного правдоподобия.
Для реализации модели Кокса в Stan, можно рассмотреть использование байесовских регрессионных моделей с приорами на коэффициенты регрессии и аппроксимации базовой функции риска, например, с использованием экспоненциального или Вейбулла распределений.
Модели кластеризации
В Stan, существует множество подходов к реализации моделей кластеризации
Модель K-средних может быть реализована через оптимизацию расстояния между точками данных и центрами кластеров.
data {
int N; // колво точек данных
int K; // предполагаемое количество кластеров
vector[N] x; // данные
}
parameters {
vector[K] mu; // предполагаемые центры кластеров
real sigma; // стандартное отклонение
}
model {
sigma ~ normal(0, 1); // априорное распределение для sigma
for (n in 1:N) {
vector[K] log_prob;
for (k in 1:K) {
log_prob[k] = normal_lpdf(x[n] | mu[k], sigma);
}
target += log_sum_exp(log_prob); // логарифм суммы экспонент логарифмических вероятностей
}
} Модель стремится найти такие значения mu и sigma, которые максимизируют правдоподобие наблюдаемых данных, что аналогично минимизации расстояния между точками данных и центрами кластеров в K-средних.
Модель гауссовской смеси позволяет кластеризовать данные, предполагая, что каждый кластер генерируется из гауссовского распределения:
data {
int K; // колво кластеров
int N; // колво наблюдений
vector[2] y[N]; // данные
}
parameters {
simplex[K] theta; // веса смеси
vector[2] mu[K]; // ср. значения кластеров
cov_matrix[2] Sigma[K]; // ковариационные матрицы
}
model {
vector[K] log_theta = log(theta);
for (n in 1:N) {
vector[K] log_phi_n;
for (k in 1:K) {
log_phi_n[k] = multi_normal_log(y[n], mu[k], Sigma[k]);
}
target += log_sum_exp(log_phi_n + log_theta);
}
} LDA является популярной моделью для кластеризации текстовых данных на основе тем:
data {
int W; // колво слов в словаре
int D; // колво документов
int N; // общее количество слов
int word[N]; // идексы слов
int doc[N]; // индексы документов
int K; // количество тем
}
parameters {
simplex[K] theta[D]; // распределение тем в документах
simplex[W] phi[K]; // расределение слов в темах
}
model {
for (n in 1:N) {
vector[K] log_theta = log(theta[doc[n]]);
vector[K] log_phi = log(phi[:, word[n]]);
target += log_sum_exp(log_theta + log_phi);
}
} Немного про визуализацию
После того как модель скомпилирована и выполнено семплирование можно использовать базовые функции R, такие как plot, для визуализации результатов. Например, для простой гистограммы параметра модели:
library(rstan)
# предположим, что `fit` - это объект, возвращенный функцией stan()
samples <- extract(fit)$parameter_name # извлекаем семплы параметра
hist(samples, breaks = 40, main = "Распределение параметра", xlab = "Значение параметра")С помощью bayesplot можно визуализировать выводы MCMC, включая трассировки, графики плотности, интервалы HPD и диагностические графики.
Трассировочные графики показывают, как значения параметра изменялись во время семплирования:
color_scheme_set("brightblue")
mcmc_trace(fit, pars = c("alpha", "beta"))Графики плотности показывают апостериорное распределение параметров:
mcmc_areas(fit, pars = c("alpha", "beta"))Диагностические графики Rhat и ESS помогают оценить качество семплирования:
mcmc_diagnostics(fit)Еще есть shinystan который интерактивный интерфейс для анализа и визуализации результатов моделирования:
# устанавливаем и запускаем shinystan
if (!requireNamespace("shinystan", quietly = TRUE)) {
install.packages("shinystan")
}
library(shinystan)
launch_shinystan(fit)Stan предоставляет интуитивно понятный язык моделирования, поддержку большего спектра статистических моделей и алгоритмы для оценки параметров.
Более подробно с документацией можно ознакомиться здесь.
А про инструменты аналитики эксперты OTUS рассказывают в рамках практических онлайн-курсов. Ознакомиться с каталогом.
