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

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<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 рассказывают в рамках практических онлайн-курсов. Ознакомиться с каталогом.

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