Обзор библиотеки Stan в R

Приветствую!

Stan — это библиотека на C++, предназначенная для байесовского моделирования и вывода. Она использует сэмплер NUTS, чтобы создавать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 состоит из нескольких блоков, каждый из которых выполняет определенную роль в процессе моделирования:

1. Блок данных (data): тут объявляются все входные данные, которые используются в модели. Это могут быть скаляры, векторы, матрицы и так далее. Данные в этом блоке считаются фиксированными и неизменными на протяжении всего процесса моделирования.

2. Блок параметров (parameters): здесь объявляются раметры модели, которые будут оцениваться на основе данных. Это могут быть, например, коэффициенты регрессии, дисперсии, масштабные коэффициенты и так далее. Stan автоматически применяет методы Монте-Карло с цепями Маркова для оценки этих параметров.

3. Блок модели (model): база программы на Stan. Здесь определяется вероятностная модель данных и параметров. Можно юзать встроенные распределения и функции Stan для описания априорных распределений параметров и вероятностей наблюдаемых данных при заданных параметрах.

4. Блоки преобразований (transformed data, transformed parameters): блоки используются для преобразования данных и параметров соответственно. Преобразованные данные могут быть юзабельны для упрощения модели или для предварительной обработки данных, а преобразованные параметры часто используются для выведения вторичных величин, которые вычисляются на основе оцененных параметров.

5. Блок генерации (generated quantities): тут можно вычислять величины, которые зависят от параметров после их оценки. Часто используются для генерации прогнозов и для вычисления статистик постериорного распределения.

Перейдем сразу к созданию, к примеру, модели линейной регрессии:

// определяем модель линейной регрессии data {   int<lower=0> N; // количество наблюдений   vector[N] y; // зависимая переменная   int<lower=0> K; // колво предикторов   matrix[N, K] X; // матрица предикторов }  parameters {   vector[K] beta; // коэф. регрессии   real alpha; // перехват (интерсепт)   real<lower=0> 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<lower=0> N;   vector[N] x;   vector[N] y; } parameters {   real alpha;   real beta;   real<lower=0> 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 <cmath>  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<lower=0> N; // колво наблюдений   vector[N] y; // временной ряд } parameters {   real alpha; // параметр авторегрессии   real beta; // коэффициент уровня   real<lower=0> 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<lower=0> N; // колво наблюдений   int<lower=0> K; // период сезонности   vector[N] y; // временной ряд } parameters {   vector[N] level; // лок. уровень   vector[K] season; // сезонные эффекты   real<lower=0> sigma_level; // шум уровня   real<lower=0> 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<lower=0> N;   vector[N] y; } parameters {   vector[N] mu; // локальный тренд   real<lower=0> 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<lower=0> N_obs; // колво наблюдаемых данных   int<lower=0> N_miss; // колво отсутствующих данных   real y_obs[N_obs]; // наблюд. данные } parameters {   real mu; // ср. значение распределения   real<lower=0> sigma; // дефолтное отклонение распределения   real y_miss[N_miss]; // отсутствующие данные, моделируемые как параметры } model {   y_obs ~ normal(mu, sigma); // модель для наблюдаемых данных   y_miss ~ normal(mu, sigma); // модель для отсутствующих данных }

y_miss моделируется так же, как и наблюдаемые данные, с теми же параметрами mu и sigma.

Использование предиктивных величин для обработки отсутствующих данных. Этот подход заключается в использовании предиктивных величин для заполнения отсутствующих данных:

data {   int<lower=0> N; // общее количество данных (наблюдаемые + отсутствующие)   int<lower=0,upper=1> missing[N]; // индикатор отсутствующих данных   real y[N]; // данные, где отсутствующие значения обозначены как NaN } parameters {   real mu;   real<lower=0> 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<lower=0> N;   real y[N];   real<lower=0> sigma_prior; // частично известное стандартное отклонение } parameters {   real mu;   real<lower=0> sigma; } model {   sigma ~ normal(0, sigma_prior); // использование частично известного параметра   y ~ normal(mu, sigma); }

sigma_prior используется для задания априорного распределения параметра sigma

Модели выживаемости

Экспоненциальная модель предполагает, что интенсивность события (или риск) является постоянной во времени. Это самая простая форма модели выживаемости:

data {   int<lower=0> N; // колво наблюдений   vector[N] T; // тайм до события или цензурирования   int<lower=0,upper=1> delta[N]; // индикатор события (1 если событие произошло, 0 если данные цензурированы) } parameters {   real<lower=0> 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<lower=0> N;   vector[N] T;   int<lower=0,upper=1> delta[N]; } parameters {   real<lower=0> lambda; // параметр масштаба   real<lower=0> 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<lower=0> N; // колво точек данных   int<lower=1> K; // предполагаемое количество кластеров   vector[N] x; // данные } parameters {   vector[K] mu; // предполагаемые центры кластеров   real<lower=0> 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<lower=1> K; // колво кластеров   int<lower=1> 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<lower=1> W; // колво слов в словаре   int<lower=1> D; // колво документов   int<lower=1> N; // общее количество слов   int<lower=1> word[N]; // идексы слов   int<lower=1> doc[N]; // индексы документов   int<lower=1> 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 рассказывают в рамках практических онлайн-курсов. Ознакомиться с каталогом.