Как перестать бояться и полюбить R. Часть 4. Установка Stan и rethinking

В моей обзорной серии про байесовскую статистику [1] мы говорили о том, что для получения марковских цепочек Монте-Карло MCMC [2] существует два инструмента JAGS и Stan. В тех статьях мы полагались на JAGS, а сегодня мы разберемся с тем, как работать со Stan.

Stan - это современная платформа для статистического моделирования и высокопроизводительных статистических вычислений [3]. Принцип работы Stan такой же как у JAGS, мы должны вручную определять модель используя довольно своеобразную нотацию. Такой подход даёт максимальную свободу, но требует усилий. Вместо этого мы будем работать с библиотекой rethinking [4], которая упрощает взаимодействие со Stan.

Первое, что понадобится сделать – это установить набор инструментов, используемый для сборки пакетов R из исходного кода RTools [5]. Скачиваем инсталлятор, устанавливаем с настройками по умолчанию.

Источник

После установки нам нужно будет внести изменения в системную переменную PATH. Открываем Start >> Edit system environment variables (или Win + R и вводим systempropertiesadvanced).

Добавляем две новые переменные, как это показано на скриншоте:

C:\rtools43\usr\bin
C:\rtools43\mingw64\bin

Настройка переменной PATH на Windows 11.

Если вы следовали инструкциям по умолчанию, то у вас должны быть такие же пути, как у меня, вы всегда можете проверить это самостоятельно в проводнике.

Переходим в R. Наверняка вы помните, как я говорил, что большая часть библиотек размещена на CRAN [6]. Большая – это не все. Многие библиотеки живут за пределами CRAN, к примеру, на GitHub. Для того, чтобы получать к ним доступ нам потребуется установить библиотеку devtools, а после мы установим Stan.

#Устанавливаем и подключаем devtools
install.packages(“devtools”)
library(devtools)

#Скачиваем и устанавливаем Stan
devtools::install_github("stan-dev/cmdstanr")
library(cmdstanr)
cmdstanr::check_cmdstan_toolchain(fix = TRUE)
install_cmdstan()

Ричард МакЭлрит [7] – автор выдающейся книги Statistical Rethinking ("Статистическое переосмысление"), про которую я уже упоминал не раз [8] и буду говорить и дальше, также разработчик библиотеки rethinking [4], которая позволяет удобно работать со Stan. При этом в составе библиотеки более широкий функционал, чем только работа с MCMC.

Установим rethinking и другие библиотеки, которые понадобятся дальше для работы в байесовском фреймворке.

install.packages(c("coda","mvtnorm","loo","dagitty","shape"))
devtools::install_github("rmcelreath/rethinking")

Настало время построить нашу первую байесовскую модель используя rethinking и Stan. Для этого мы воспользуемся функцией ulam, названной в честь Станислава Улама [9], одного из основоположников MCMC.

Предположим, что мы хотим проверить гипотезу о том, какой стаж характерен для сотрудников нашей компании. Пусть наши априорные представления будут равны среднему стажу в 30 месяцев и стандартному отклонению в 5 месяцев.

# Подключаем библиотеку rethinking
library(rethinking)

# Сгенирируем данные с помощью нормального распределения и запишем результат в датафрейм
df <- rnorm(1000, 24, 6) |> as.data.frame()
colnames(df) <- "tenure"

# Определим модель
model <- ulam(
  alist(
    tenure ~ dnorm(mu, sigma), # используем нормальное распределение для расчёта апостериорного распределения с параметрами mu и sigma
    mu ~ dnorm(30, 5), # наше априорное представление о стаже
    sigma ~ dexp(1) # наше априорное представление о вариативности в популяции. Используем экспоненциальное распределение с 1, как неопределённое представление. 
  ), 
  data = df,  # датафрейм с данными
  chains = 4, # количество цепочек MCMC, которые считаются в параллели 
  cores = 4, # количество ядер процессора, которые используем для параллельных вычислений
  iter = 2000, # число итераций по расчёту MCMC
  warmup = 1000 # сколько итераций закладываем на «разогрев» модели
)

Про модели мы ещё будем говорить подробнее в будущем, поэтому не переживайте, если вам не понятно. Расчёт цепочек MCMC потребует какого-то времени, в зависимости от производительности вашей машины.

На двух из трех компьютеров, на которых я устанавливал все требуемые ресурсы, при запуске расчёта модели я получал ошибку про OneDrive. Не знаю, как это лечить правильно, но если вы столкнулись с таким же, то просто удалите ярлык OneDrive из документов.

# Пример ошибки с OneDrive

Exit status: 2
Stderr:
make: *** /c/Users/{myusername}/OneDrive: Is a directory.  Stop.

После окончания расчёта мы можем посмотреть на результат модели.

precis(model)
       mean   sd  5.5% 94.5% n_eff Rhat4
mu    24.51 0.19 24.18 24.82  2589     1
sigma  6.14 0.14  5.93  6.37  2411     1

Итак, модель нам говорит, что апостериорное среднее = 24.5 и стандартное отклонение = 6 с узкими 89% доверительным интервалом (credible interval).

Ссылки