Bayesian Inference: Install RTools and package ‘rstan’
Alexander Rodionov
23 Апрель 2017
1 Как инсталлировать библиотеку Stan и начать использовать в R
Для работы Библиотеки Языка моделирования Байесовских выводов Stan необходимо предварительно инсталлировать программу-дополнение R - RTools. Этот набор инструментов содержит компилятор языка C++, который требуется Stan. Хотя Андрей Горюшев описал этот процесс, он далеко не всегда работает. Необходимо, чтобы RTools находился только на диске C: и всенепременно в одноименной папке. Если хотя бы одно из этих условий не будет соблюдено, то уже никогда не вызвать Stan из R посредством пакета rstan 
или попросу строить разнообразные байесовские регрессии функциями пакета rstanarm.
Однако решить проблему конфигурации доступа к C++ из R не сложно. На английском языке об этом можно прочитать на сайте разработки пакета ‘rstan’.
Для вызова функций библиотеки Stan нужно корректно описать пути к компилятора C++ из RTools, которые указаны в переменной окружения операционной системы 'PATH', а затем объяснить R, где помещен этот компилятор.
1.1 Инсталляция Rtools на Windows
Проделайте все шаги инсталляции до пункта Установка RStan (наконец-то!). При этом вы можете устанавливать Rtools куда угодно на локальных дисках вашего компьютера. Например, D:\R\Rtools.
При этом на этапе инсталляции (в окне “Выберите дополнительные задачи”) отметьте галочкой опцию “Edit the system PATH” (нужно, чтобы впоследствии R “увидел” компилятор C++) и введите именно тот путь, где реально будет располагаться Rtools:
Если у вас ранее был установлен Rtools, то старые пути лучше на этом шаге удалить из переменной окружения операционной системы 'PATH'.
Перезапускаем R/RStudio и вбиваем на консоле следующие команды:
Sys.getenv('PATH')
system('g++ -v')
system('where make')
Без перегрузки R/RStudio они сработают не верно. А после вы должны получить нечто подобное:
Консоль в R с путями
1.2 Настройка путей к компилятору C++ в R
Однако при такой попытке обратиться к C++ получается нечто подобное:
По умолчанию R полагает, что компилятор C++ помещен в папке C:/Rtools/mingw_64/bin, а ведь он лежит в D:/R/Rtools/mingw_64/bin. Чтобы переубедить программу необходимо в конфигурационной файле .Rprofile, помещаемом в домашней папке с документами пользователя Windows, добавить строку с описание пути к компилятор C++ от RTools, в нашем случае в переменную BINPREF из окружения R - D:/R/Rtools/mingw_$(WIN)/bin/:
cat("Sys.setenv(BINPREF = \"D:/R/Rtools/mingw_$(WIN)/bin/\")", file = file.path(Sys.getenv("HOME"),
".Rprofile"), sep = "\n", append = TRUE)1.3 Проверка корректной работы компилятора C++ от RTools
После этого следует вновь перегрузить R/RStudio. Теперь он узнает о новом пути к компилятору C++ от RTools. Потом в консоле проверим работу компилятора простым примером функции на C++:
fx <- inline::cxxfunction( signature(x = "integer", y = "numeric" ) , '
return ScalarReal( INTEGER(x)[0] * REAL(y)[0] ) ;
' )
fx( 2L, 5 ) # should be 10## [1] 10Если все нормально, то мы “попадаем” в десятку.
1.4 Инсталляция Stan в R/Rstudio
Теперь можно проинсталлировать пакет rstan с библиотекой Stan:
# check if packages are installed
if ( !is.element("rstan", installed.packages()[,1]) ) {
install.packages("rstan", dependencies = TRUE)
}Можно ожидать, что библиотека Stan со всеми связанными пакетами благополучно проинсталлирована. Настал момент вновь перегрузить R/Rstudio.
1.5 Пакеты R работающие c моделями от Stan
Для работы с Stan могут быть полезны несколько обращающихся к rstan пакетов: bayesplot, beanz, bmlm, brms, ctsem, DeLorean, dfpk, DrBats, hBayesDM, MCMCvis, prophet, rstanarm, rstantools, shinystan. Они существенно облегчают применение этой библиотеки в R/RStudio.
Для построения моделей из них особо популярны пакеты:
rstanarm
- пакет от американских создателей rstanдает доступ к предварительно скомпилированным вrstanрегрессионным моделям для Байесовской оценки. Модели вызываются при помощи традиционного синтаксисаR: формулой с набором данных плюс некоторых дополнительных аргументов для задания одного из 7-ми семейств априорных вероятностных распределений.brms
- относительно молодой конкурент предыдущего пакета, строит на базе rstanболее сложные Байесовские обобщенные (не-)линейные многоуровневые модели. По утверждению его германских создателей, построитель моделей оперирует около 30-тью семействами априорных распределений.prophet
- новинка 2017 г. реализует процедуры Байесовского прогнозирования временных рядов на основе основе Обобщенной Аддитивной Модели (GAM от англ. "Generalized Additive Model", предложенной в 1987), учитывающей (не-)линейные тренды и различные сезональности от американской команды Facebook’s Core Data Science team.
Графический анализ полученных моделей проводят пакетами:
bayesplot
- построение графиков для апостериорного анализа, проверки модели и диагностики MCMC (от англ. "Markov Chain Monte Carlo methods"- методы Монте-Карло с цепями Маркова). Пакет разработан не только для интерфейса сrstan, но и для других пакетов формирующих байесовские выводы. Его загружают для использования моделирующие пакетыrstanarmиbrms.coda
- cодержит функции для суммирования и построения байесовского вывода при моделирования MCMC, а также диагностические тесты сходимости к равновесному распределению цепи Маркова. Его загружает и использует моделирующий пакет brms.shinystan
- интерактивная визуальная и численная диагностика и апостериорный анализ для байесовских моделей, имеющий интерфейс к графикам из пакета bayesplot. Она также вызывается пакетамиrstanarmиbrms.
Кроме того, будет полезен ggmcmc 
- этот графический пакет на основе библиотек ggplot2 
, способен работать не только со Stan, на также с классом mcmc.list из пакета R2jags 
(другая Байесовская библиотека JAGS - Just Another Gibbs Sampler) и классом mcmc из пакета MCMCpack 
(еще одна Байесовская библиотека, но на основе Scythe Statistical Library).
1.6 Проверка работы Stan в R
Не будем обращаться напрямую к библиотеке Stan, обладающей довольно сложным синтаксисом. Для начала воспользуется пакетом rstanarm с более привычным для R синтаксисом моделирования.
# check if packages are installed
if ( !is.element("rstanarm", installed.packages()[,1]) ) {
install.packages("rstanarm", dependencies=TRUE)
}
library("rstan")## Loading required package: ggplot2## Loading required package: StanHeaders## rstan (Version 2.15.1, packaged: 2017-04-19 05:03:57 UTC, GitRev: 2e1f913d3ca3)## For execution on a local, multicore CPU with excess RAM we recommend calling
## rstan_options(auto_write = TRUE)
## options(mc.cores = parallel::detectCores())library("rstanarm")## Loading required package: Rcpp## rstanarm (Version 2.15.2, packaged: 2017-04-19 09:04:44 UTC)## - Do not expect the default priors to remain the same in future rstanarm versions.## Thus, R scripts should specify priors explicitly, even if they are just the defaults.## - For execution on a local, multicore CPU with excess RAM we recommend calling## options(mc.cores = parallel::detectCores())# How many cores can be used in Stan for parallel computing
rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())
# There are cores on this computer
parallel::detectCores()## [1] 8Для параллельных вычислений на центральном процесссоре применяемого компьютера доступно 8 ядер. Построим Байесовскую линейную множественную регрессионную модель на широкоизвестном в R наборе данных mtcars (1974) из пакета datasets.
Для описания сколько американских миль проезжали модели легковых автомобилей 1973-1974 гг. выпуска израсходовав 1 галлон бензина в качестве формулы зададим следующую зависимость:
mpg ~ hp + qsec + factor(am), где
wt- вес автомобиля, тыс. фунтов;qsec- за сколько секунд автомобиль проезжал 1/4 американской мили;am- тип передачи скоростей (0 = автоматическая, 1 = ручная), применяемый в модели как бинарный признак.
Это простая модель - даже утверждение об семействе априорного вероятностного распределения опущено. В линейной регрессионой модели обычно используют нормальное, иначе гаусовское распределение. Для повторного воспроизведения и корректного сравнения зададим одинаковый параметр seed. Он регулирует генератор, встроенный в R, так чтобы случайные числа в этих моделях формировались по одной и той же “последовательности”.
data(mtcars)
## create a bayesian regression model by rstanarm
fit.bayes1 <- stan_lm(mpg ~ wt + qsec + factor(am),
data = mtcars,
prior = NULL,
seed = 12345)
## generate a summary of the results
summary(fit.bayes1, digits = 4)##
## Model Info:
##
## function: stan_lm
## family: gaussian [identity]
## formula: mpg ~ wt + qsec + factor(am)
## algorithm: sampling
## priors: see help('prior_summary')
## sample: 4000 (posterior sample size)
## num obs: 32
##
## Estimates:
## mean sd 2.5% 25% 50% 75%
## (Intercept) 9.9277 7.1671 -4.7176 5.2894 10.0750 14.7682
## wt -3.8010 0.7353 -5.2931 -4.2701 -3.8018 -3.3260
## qsec 1.1905 0.3025 0.6066 0.9860 1.1867 1.3873
## factor(am)1 2.8285 1.4770 -0.1091 1.8846 2.8094 3.8067
## sigma 2.5912 0.3439 2.0139 2.3466 2.5673 2.7980
## log-fit_ratio -0.0035 0.0767 -0.1521 -0.0571 -0.0054 0.0481
## R2 0.8106 0.0514 0.6917 0.7824 0.8181 0.8477
## mean_PPD 20.0847 0.6631 18.7765 19.6446 20.0874 20.5286
## log-posterior -78.3880 2.0795 -83.2444 -79.6067 -78.0732 -76.8688
## 97.5%
## (Intercept) 23.7750
## wt -2.3462
## qsec 1.8138
## factor(am)1 5.7118
## sigma 3.3442
## log-fit_ratio 0.1483
## R2 0.8875
## mean_PPD 21.3503
## log-posterior -75.3358
##
## Diagnostics:
## mcse Rhat n_eff
## (Intercept) 0.1935 1.0005 1371
## wt 0.0157 1.0008 2206
## qsec 0.0104 1.0010 853
## factor(am)1 0.0316 1.0015 2180
## sigma 0.0072 1.0006 2286
## log-fit_ratio 0.0020 1.0018 1436
## R2 0.0012 1.0035 1796
## mean_PPD 0.0105 0.9993 4000
## log-posterior 0.0885 1.0042 553
##
## For each parameter, mcse is Monte Carlo standard error, 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).## mean coefficent R2
r2 <- mean(c(
fit.bayes1$stanfit@sim$samples[[1]]$`R2[1]`,
fit.bayes1$stanfit@sim$samples[[2]]$`R2[1]`,
fit.bayes1$stanfit@sim$samples[[3]]$`R2[1]`,
fit.bayes1$stanfit@sim$samples[[4]]$`R2[1]`
))
cat(paste("\nMean of R2 =", round(r2, digits=4)))##
## Mean of R2 = 0.8106Среднеквадратичное отклонение для байесовской регрессионной модели составляет R2=0.8106. Ее значение достаточно близко единице, что может говорить о высокой значимости модели.
Обладающий расширенными возможностями немецкий пакет brms с библиотекой Stan дает схожие результаты. В этом примере также пропустим большинство параметров, применяемых для построения байесовских моделей:
# check if packages are installed
if ( !is.element("brms", installed.packages()[,1]) ) {
install.packages("brms", dependencies=TRUE)
}
library("brms")## Loading 'brms' package (version 1.6.1). Useful instructions
## can be found by typing help('brms'). A more detailed introduction
## to the package is available through vignette('brms_overview').##
## Attaching package: 'brms'## The following objects are masked from 'package:rstanarm':
##
## exponential, lasso, ngrps## create a bayesian regression model by brms
fit.bayes2 <- brm(mpg ~ wt + qsec + factor(am),
data = mtcars,
prior = NULL,
seed = 12345)## Compiling the C++ model## Start sampling## generate a summary of the results
summary(fit.bayes2)## Family: gaussian(identity)
## Formula: mpg ~ wt + qsec + factor(am)
## Data: mtcars (Number of observations: 32)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
## WAIC: Not computed
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## Intercept 9.37 7.56 -5.86 24.27 1981 1
## wt -3.89 0.78 -5.40 -2.35 2094 1
## qsec 1.23 0.31 0.63 1.85 2384 1
## factoram1 2.98 1.55 -0.07 6.11 2208 1
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat
## sigma 2.58 0.38 1.96 3.45 2850 1
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).Эта модель также построена на четырех цепях Маркова по 2000 итераций, из которых первая тысяча применялась для “разогрева” построителя апостериорного вероятностного распределения.
Однако пакет brms для созданной им регрессионной модели не рассчитывает значения R2, ограничиваясь лишь параметрами семейства вероятностного распределения. Для линейной регрессии, остатки по которой описывают нормальным распределением, это simga. Германские разработчики этого пакета справедливо полагают, что параметра simga, который в этой модели выступает Residual Standard Error достаточно, чтобы рассчитать все остальные коэффициенты, описывающие качество построенной линейной регрессионной модели.
#
# www.learnbymarketing.com/tutorials/explaining-the-lm-summary-in-r/
# en.wikipedia.org/wiki/Anscombe%27s_quartet - four datasets for testing
#
cat("Bayesian Linear Regression Model by 'brms' package:\n")
n = length(mtcars$mpg) # Nums of Observations
k = fit.bayes2$fit@par_dims$b # Nums of Predictors(Beta Coefficients)
SSE = (summary(fit.bayes2)$spec_par[1, 1]**2)*(n-(1+k))
SSyy = sum((mtcars$mpg - mean(mtcars$mpg))**2)
# Residual Standard Error (Like Standard Deviation)
RSE = sqrt(SSE / (n - (1 + k))) #Residual Standard Error
cat(paste("\nResidual standard error:", round(RSE, digits=4), "on",
(n - (1 + k)), "degrees of freedom"))
# Multiple R-Squared: Percent of the variance of Y intact after subtracting the error of the model
r2 <- 1 - (SSE / SSyy)
cat(paste("\nMultiple R-squared: ", round(r2, digits=3)))
# Adjusted R-Squared: Same as multiple R-Squared but takes into account the number of samples and variables you’re using
r2adj <- 1 - (SSE / SSyy) * (n - 1) / (n - (k + 1))
cat(paste(", Adjusted R-squared: ", round(r2adj, digits = 4)))
# F-Statistic: Global test to check if your model has at least one significant variable.
# Ho: All coefficients are zero
# Ha: At least one coefficient is nonzero
# Compare test statistic to F Distribution table
Fstat = ((SSyy - SSE) / k) / (SSE / (n - (k + 1)))
p = 1 - pf(Fstat, df1 = k, df2 = (n - (1 + k))) # p-value for an F-Test
cat(paste("\nF-statistic:", round(Fstat, digits=2), "on", k, "and",
n - (1 + k), "DF, p-value:", formatC(p)))## Bayesian Linear Regression Model by 'brms' package:
##
## Residual standard error: 2.5796 on 28 degrees of freedom
## Multiple R-squared: 0.835, Adjusted R-squared: 0.8168
## F-statistic: 47.07 on 3 and 28 DF, p-value: 4.597e-11Очевидно, что полученная пакетом brms регрессия имеет значимость среднеквадратичного отклонения (‘Adjusted R-Squared’) равного 0.8168 подтверждается высоким значением F-статистики, равным 47.07 и общий уровень значимости (определяемый по этой статистике): ‘p-value’: 4.597e-11, что много меньше 0.001.
Если на основе этого анализа будет составлен отчет, то надо будет указать также значения степеней свободы df: 3 и 28, т.е. по колонкам (предикторам модели) и по строкам (количеству наблюдений минус число предикторов и свободный член уравнения) данных соответственно.
1.7 Графическая диагностика байесовских моделей
На созданные здесь байесовские модели получили графическую диагностику в пакете bayesplot. Ограничимся графиком плотности апостериорного распределения всех параметров модели, полученных в rstanarm.
library("bayesplot")## This is bayesplot version 1.2.0posterior <- as.array(fit.bayes1)
dim(posterior)## [1] 1000 4 7dimnames(posterior)## $iterations
## NULL
##
## $chains
## [1] "chain:1" "chain:2" "chain:3" "chain:4"
##
## $parameters
## [1] "(Intercept)" "wt" "qsec" "factor(am)1"
## [5] "sigma" "log-fit_ratio" "R2"bayesplot::mcmc_areas(
posterior,
prob = 0.8 # 80% Credible Interval (CrI)
)
Очевидно, что рассмотрение всех параметров при одинаковом размахе оси абсцисс малопродуктивно. Поэтому оставляем только предикторы и переменную simga - среднеквадратическое (стандартное) отклонение, описывающее полученное в модели апостериорное нормальное распределение.
Оба моделирующие пакета, использованные в этом примере, самостоятельно обращаются к функциям пакета bayesplot. Это облегчает графическое сравнение полученных моделей, ведь объекты stanreg и brmsfit от этих пакетов имеют различную структуру.
Обратите внимание, что предикторы в модели, созданной пакетом brms получили префиксы 'b_' от греческой буквы “beta”, обозначающей коэффициенты регрессии предикторов.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
## Scale for 'x' is already present. Adding another scale for 'x', which
## will replace the existing scale.
Поскольку на оси абсцисс двух графиков задано одинаковое ограничение xlim(-5.5, 7.5), то хорошо видно, что медианы апостериорного распределения у этих байесовских моделей практически совпадают.
1.8 Описание модели Байесовских выводов на языке моделирования Stan
И наконец, пакет brms способен вызовом функцией извлечь из объекта класса brmsfit объект model - компактное описание модели на языке моделирования Байесовских выводов Stan (см. объект fit.bayes2$model) с комментариями, помогающими его изучению:
# Extract the 'brms' model code in Stan language
# cat (fit.bayes2$model)
stancode(fit.bayes2)## // generated with brms 1.6.1
## functions {
## }
## data {
## int<lower=1> N; // total number of observations
## vector[N] Y; // response variable
## int<lower=1> K; // number of population-level effects
## matrix[N, K] X; // population-level design matrix
## int prior_only; // should the likelihood be ignored?
## }
## transformed data {
## int Kc;
## matrix[N, K - 1] Xc; // centered version of X
## vector[K - 1] means_X; // column means of X before centering
## Kc = K - 1; // the intercept is removed from the design matrix
## for (i in 2:K) {
## means_X[i - 1] = mean(X[, i]);
## Xc[, i - 1] = X[, i] - means_X[i - 1];
## }
## }
## parameters {
## vector[Kc] b; // population-level effects
## real temp_Intercept; // temporary intercept
## real<lower=0> sigma; // residual SD
## }
## transformed parameters {
## }
## model {
## vector[N] mu;
## mu = Xc * b + temp_Intercept;
## // prior specifications
## sigma ~ student_t(3, 0, 10);
## // likelihood contribution
## if (!prior_only) {
## Y ~ normal(mu, sigma);
## }
## }
## generated quantities {
## real b_Intercept; // population-level intercept
## b_Intercept = temp_Intercept - dot_product(means_X, b);
## }В пакете rstanarm для этого придется вывести один из вложенных объектов класса stanreg, порождаемого его регрессионной функцией. Однако, это будет весьма пространное описание: оно в четыре раза превышает аналогичную модель от германского конкурента.
# Extract the 'rstanarm' model code in Stan language
cat(fit.bayes1$stanfit@stanmodel@model_code)##
## // This file is part of rstanarm.
## // Copyright (C) 2015, 2016 2017 Trustees of Columbia University
##
##
## /*
## rstanarm is free software: you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation, either version 3 of the License, or
## (at your option) any later version.
##
## rstanarm is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with rstanarm. If not, see <http://www.gnu.org/licenses/>.
## */
##
## // GLM for a Gaussian outcome with no link function
## functions {
## /**
## * Increments the log-posterior with the logarithm of a multivariate normal
## * likelihood with a scalar standard deviation for all errors
## * Equivalent to normal_lpdf(y | intercept + Q * R * beta, sigma) but faster
## * @param theta vector of coefficients (excluding intercept), equal to R * beta
## * @param b precomputed vector of OLS coefficients (excluding intercept) in Q-space
## * @param intercept scalar (assuming columns of Q have mean zero)
## * @param ybar precomputed sample mean of the outcome
## * @param SSR positive precomputed value of the sum of squared OLS residuals
## * @param sigma positive scalar for the standard deviation of the errors
## * @param N integer equal to the number of observations
## */
## real ll_mvn_ols_qr_lp(vector theta, vector b,
## real intercept, real ybar,
## real SSR, real sigma, int N) {
## target += -0.5 * (dot_self(theta - b) +
## N * square(intercept - ybar) + SSR) /
## square(sigma) -// 0.91... is log(sqrt(2 * pi()))
## N * (log(sigma) + 0.91893853320467267);
## return target();
## }
## }
## data {
## int<lower=0,upper=1> has_intercept; // 0 = no, 1 = yes
## int<lower=0,upper=1> prior_dist_for_intercept; // 0 = none, 1 = normal
## real<lower=0> prior_scale_for_intercept; // 0 = by CLT
## real prior_mean_for_intercept; // expected value for alpha
## int<lower=0,upper=1> prior_dist; // 0 = uniform for R^2, 1 = Beta(K/2,eta)
## int<lower=0,upper=1> prior_PD; // 0 = no, 1 = yes to drawing from the prior
## real<lower=0> eta; // shape hyperparameter
##
## int<lower=1> J; // number of groups
## // the rest of these are indexed by group but should work even if J = 1
## int<lower=1> N[J]; // number of observations
## int<lower=1,upper=min(N)> K; // number of predictors
## vector[K] xbarR_inv[J]; // vector of means of the predictors
## real ybar[J]; // sample mean of outcome
## real center_y; // zero or sample mean of outcome
## real<lower=0> s_Y[J]; // standard deviation of the outcome
## vector[K] Rb[J]; // OLS coefficients
## real<lower=0> SSR[J]; // OLS sum-of-squared residuals
## matrix[K,K] R_inv[J]; // inverse R matrices
## }
## transformed data {
## real half_K = 0.5 * K;
## real sqrt_inv_N[J];
## real sqrt_Nm1[J];
## for (j in 1:J) {
## sqrt_inv_N[j] = sqrt(1.0 / N[j]);
## sqrt_Nm1[j] = sqrt(N[j] - 1.0);
## }
## }
## parameters { // must not call with init="0"
## unit_vector[K] u[J]; // primitives for coefficients
## real z_alpha[J * has_intercept]; // primitives for intercepts
## real<lower=0,upper=1> R2[J]; // proportions of variance explained
## vector[J * (1 - prior_PD)] log_omega; // under/overfitting factors
## }
## transformed parameters {
## real alpha[J * has_intercept]; // uncentered intercepts
## vector[K] theta[J]; // coefficients in Q-space
## real<lower=0> sigma[J]; // error standard deviations
## for (j in 1:J) {
## // marginal standard deviation of outcome for group j
## real Delta_y = prior_PD == 0 ? s_Y[j] * exp(log_omega[j]) : 1;
##
## // coefficients in Q-space
## if (K > 1) theta[j] = u[j] * sqrt(R2[j]) * sqrt_Nm1[j] * Delta_y;
## else theta[j][1] = u[j][1] * sqrt(R2[j]) * sqrt_Nm1[j] * Delta_y;
##
## sigma[j] = Delta_y * sqrt(1 - R2[j]); // standard deviation of errors
##
## if (has_intercept == 1) {
## if (prior_dist_for_intercept == 0) // no information
## alpha[j] = z_alpha[j];
## else if (prior_scale_for_intercept == 0) // central limit theorem
## alpha[j] = z_alpha[j] * Delta_y * sqrt_inv_N[j] + prior_mean_for_intercept;
## else // arbitrary informative prior
## alpha[j] = z_alpha[j] * prior_scale_for_intercept +
## prior_mean_for_intercept;
## }
## }
## }
## model {
## for (j in 1:J) { // likelihood contribution for each group
## if (prior_PD == 0) {
## real dummy; // irrelevant but useful for testing user-defined function
## real shift;
## shift = dot_product(xbarR_inv[j], theta[j]);
## dummy = ll_mvn_ols_qr_lp(theta[j], Rb[j],
## has_intercept == 1 ? alpha[j] + shift : shift,
## ybar[j], SSR[j], sigma[j], N[j]);
## }
## // implicit: u[j] is uniform on the surface of a hypersphere
## }
## if (has_intercept == 1 && prior_dist_for_intercept > 0)
## target += normal_lpdf(z_alpha | 0, 1);
## if (prior_dist == 1) target += beta_lpdf(R2 | half_K, eta);
## // implicit: log_omega is uniform over the real line for all j
## }
## generated quantities {
## real mean_PPD[J];
## vector[K] beta[J];
## for (j in 1:J) {
## real shift;
## shift = dot_product(xbarR_inv[j], theta[j]);
## mean_PPD[j] = normal_rng(has_intercept == 1 ? alpha[j] + shift : shift,
## sigma[j] * sqrt_inv_N[j]);
## beta[j] = R_inv[j] * theta[j];
## }
## }