R包rethinking的介绍

参考网址：https://github.com/rmcelreath/rethinking

这个包附带了一个关于贝叶斯数据分析的课程和书籍:McElreath 2020。《Statistical Rethinking》第二版，CRC出版社。如果您正在使用第一版的书，请参阅该文件底部的注释。
它包含了对后验分布进行快速二次逼近和哈密顿蒙特卡罗(通过RStan - mc-stan.org)的工具。许多包都这样做。此包的签名差异在于，它强制用户将模型指定为显式分布假设的列表。这比典型的基于公式的工具更加繁琐，但是它也更加灵活和强大，并且——最重要的是——对于教学和学习非常有用。当学生必须写出模型的每个细节时，他们实际上是在学习模型。
例如，一个简单的高斯模型可以用这个公式列表来指定:

f <- alist(
    y ~ dnorm( mu , sigma ),
    mu ~ dnorm( 0 , 10 ),
    sigma ~ dexp( 1 )
)

列表中的第一个公式是结果的概率(可能性);第二个是mu的先验;第三个是sigma的先验。

快速安装

您可以在http://xcelab.net/rm/software/中找到带有扩展安装和使用说明的手册
这是一个简单的版本。
您需要首先安装rstan。访问http://mc-stan.org并遵循您的平台的说明。最大的挑战是使一个c++编译器配置为与r的安装一起工作。https://github.com/stan-dev/rstan/wiki/RStan-Getting-Started上的说明非常全面。遵守它们，你就有可能成功。
然后你可以在Rstudio上安装rethinking:

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

如果存在任何问题，那么在尝试安装rstan时可能会出现问题，因此重新考虑包与此关系不大。有关安装rstan的一些提示，请参阅上面链接的手册。但是，请始终在mc-stan.org网站的RStan部分查阅关于RStan的最新信息。

quap的二次逼近

几乎任何普通的广义线性模型都可以用quap来指定。使用二次逼近:

library(rethinking)

f <- alist(
    y ~ dnorm( mu , sigma ),
    mu ~ dnorm( 0 , 10 ),
    sigma ~ dexp( 1 )
)

fit <- quap( 
    f , 
    data=list(y=c(-1,1)) , 
    start=list(mu=0,sigma=1)
)

对象fit保存结果。对于边缘后验分布的总结，使用summary(fit)或precis(fit):

mean   sd  5.5% 94.5%
mu    0.00 0.59 -0.95  0.95
sigma 0.84 0.33  0.31  1.36

它还支持向量化的参数，这对于分类很方便。看到例子吗? quap。
在教科书的第一版中，这个函数被称为map。它仍然可以与别名一起使用。它被重命名了，因为名称映射有误导性。这个函数产生后验分布的二次近似，而不仅仅是最大后验(MAP)估计。

哈密顿蒙特卡罗与ulam(和map2stan)
可以使用以下两种工具之一将相同的公式列表编译成Stan (mc-stan.org)模型:ulam或map2stan。对于简单的模型，它们是相同的。ulam是一个较新的工具，它提供了更大的灵活性，包括显式变量类型和定制发行版。map2stan是该软件包第一版和教科书中的原始工具。未来，ulam将增加新的功能。
Ulam命名来源于Stanisław Ulam, 蒙特卡罗方法的父母之一,也是斯坦的同名项目。它的发音类似[o -lahm]，而不是[YOU-lamm]。
这两个工具采用与quap相同的输入类型:

fit_stan <- ulam( f , data=list(y=c(-1,1)) )

只要安装了rstan，链条就会自动运行。诊断链显示在precis(fit_stan)输出:

      mean   sd  5.5% 94.5% n_eff Rhat
sigma 1.45 0.72  0.67  2.84   145    1
mu    0.12 1.04 -1.46  1.59   163    1

对于ulam模型，plot显示的信息与precis和traceplot显示的链相同。
extract.samples 返回列表中的样本。extract.prior 从先验中获取先验样本并返回一个列表中的样本。
stanfit对象本身位于@stanfit槽中。用Stan模型做的任何事情都可以直接用那个槽完成。
Stan代码可以通过使用stancode(fit_stan)访问:

data{
    real y[2];
}
parameters{
    real<lower=0> sigma;
    real mu;
}
model{
    sigma ~ exponential( 1 );
    mu ~ normal( 0 , 10 );
    y ~ normal( mu , sigma );
}

注意，ulam并不关心R分布名。你可以使用stan风格的名字:

fit_stan <- ulam(
    alist(
        y ~ normal( mu , sigma ),
        mu ~ normal( 0 , 10 ),
        sigma ~ exponential( 1 )
    ), data=list(y=c(-1,1)) )

后验预测
所有quap、ulam和map2stan对象都可以进行后处理以生成后验预测分布。
link 用于计算任何线性模型的后验分布样本上的值。
用 sim模拟后验预测分布，模拟参数后验分布对样本的结果。sim还可以用来模拟先验。
详情请参见?link和?sim。
postcheck自动计算后验预测(后验?)检查。它只使用link和sim。

多层次模型公式
虽然quap在大多数情况下仅限于固定效果模型，但是ulam可以指定多级模型，甚至非常复杂的模型。例如，一个简单的多种截距模型如下:

# prep data
data( UCBadmit )
UCBadmit$male <- as.integer(UCBadmit$applicant.gender=="male")
UCBadmit$dept <- rep( 1:6 , each=2 )
UCBadmit$applicant.gender <- NULL

# varying intercepts model
m_glmm1 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- a[dept] + b*male,
        a[dept] ~ normal( abar , sigma ),
        abar ~ normal( 0 , 4 ),
        sigma ~ half_normal(0,1),
        b ~ normal(0,1)
    ), data=UCBadmit )

相似多种斜率模型为:

m_glmm2 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- a[dept] + b[dept]*male,
        c( a , b )[dept] ~ multi_normal( c(abar,bbar) , Rho , sigma ),
        abar ~ normal( 0 , 4 ),
        bbar ~ normal(0,1),
        sigma ~ half_normal(0,1),
        Rho ~ lkjcorr(2)
    ),
    data=UCBadmit )

另一种表示变斜率模型的方法是用一个具有不同效果的向量。这可以通过在公式中使用一个显式的向量声明来实现:

m_glmm3 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- v[dept,1] + v[dept,2]*male,
        vector[2]:v[dept] ~ multi_normal( c(abar,bbar) , Rho , sigma ),
        abar ~ normal( 0 , 4 ),
        bbar ~ normal(0,1),
        sigma ~ half_normal(0,1),
        Rho ~ lkjcorr(2)
    ),
    data=UCBadmit )

vector[2]:v[dept] 表示“为每个惟一的dept声明一个长度为2的向量”。为了访问这些向量的元素，线性模型使用括号内的多个索引:[dept,1]。
这一策略可以进一步发展，其手段也可以作为一个载体宣布:

m_glmm4 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- v[dept,1] + v[dept,2]*male,
        vector[2]:v[dept] ~ multi_normal( v_mu , Rho , sigma ),
        vector[2]:v_mu ~ normal(0,1),
        sigma[1] ~ half_normal(0,1),
        sigma[2] ~ half_normal(0,2),
        Rho ~ lkjcorr(2)
    ),
    data=UCBadmit )

完全非中心参数化也可以直接编码:

m_glmm5 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- v_mu[1] + v[dept,1] + (v_mu[2] + v[dept,2])*male,
        matrix[dept,2]: v <- t(diag_pre_multiply( sigma , L_Rho ) * z),
        matrix[2,dept]: z ~ normal( 0 , 1 ),
        vector[2]: v_mu[[1]] ~ normal(0,4),
        vector[2]: v_mu[[2]] ~ normal(0,1),
        vector[2]: sigma ~ half_normal(0,1),
        cholesky_factor_corr[2]: L_Rho ~ lkj_corr_cholesky( 2 )
    ),
    data=UCBadmit )

在上面的例子中，变化效果矩阵v是由一个z-scores z的矩阵和包含在sigma中的协方差结构以及一个Cholesky因子L_Rho构成的。注意双括号符号v_[[1]]允许对向量的每个索引使用不同的先验。

WAIC和LOOCV的对数似然计算
ulam可以有选择地返回点态log-likelihood值。这些是计算WAIC和PSIS-LOO所需要的。log_lik参数打开这个:

m_glmm1 <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- a[dept] + b*male,
        a[dept] ~ normal( abar , sigma ),
        abar ~ normal( 0 , 4 ),
        sigma ~ half_normal(0,1),
        b ~ normal(0,1)
    ), data=UCBadmit , log_lik=TRUE )
WAIC(m_glmm1)

额外的代码已经添加到生成的数量块的Stan模型(看看这个与stancode(m_glmm1)):

generated quantities{
    vector[12] log_lik;
    vector[12] p;
    for ( i in 1:12 ) {
        p[i] = a[dept[i]] + b * male[i];
        p[i] = inv_logit(p[i]);
    }
    for ( i in 1:12 ) log_lik[i] = binomial_lpmf( admit[i] | applications[i] , p[i] );
}

条件语句、自定义分布和混合模型
ulam还支持if-then语句和自定义分配。这对于编码混合模型(如zero-inflated Poisson模型和离散缺失值模型)非常有用。
这是一个zero-inflated Poisson 模型的例子。

# zero-inflated poisson
# gen data first - example from text
prob_drink <- 0.2 # 20% of days
rate_work <- 1    # average 1 manuscript per day
N <- 365
drink <- rbinom( N , 1 , prob_drink )
y <- as.integer( (1-drink)*rpois( N , rate_work ) )
x <- rnorm( N ) # dummy covariate

# now ulam code
m_zip <- ulam(
    alist(
        y|y==0 ~ custom( log_mix( p , 0 , poisson_lpmf(0|lambda) ) ),
        y|y>0 ~ custom( log1m(p) + poisson_lpmf(y|lambda) ),
        logit(p) <- ap,
        log(lambda) <- al + bl*x,
        ap ~ dnorm(0,1),
        al ~ dnorm(0,10),
        bl ~ normal(0,1)
    ) ,
    data=list(y=y,x=x) )

上式前两行对应的Stan代码为:

for ( i in 1:365 ) 
    if ( y[i] > 0 ) target += log1m(p) + poisson_lpmf(y[i] | lambda[i]);
for ( i in 1:365 ) 
    if ( y[i] == 0 ) target += log_mix(p, 0, poisson_lpmf(0 | lambda[i]));

custom 做的是定义 custom target 更新。而“ | ” 操作符使该行成为有条件的。注意，log1m、log_mix和poisson_lpmf是Stan函数。
相同的 custom 分布方法允许在离散缺失值上进行边缘化。让我们引入一些先前 UCBadmit 数据中缺失的值。

UCBadmit$male2 <- UCBadmit$male
UCBadmit$male2[1:2] <- (-1) # missingness code
UCBadmit$male2 <- as.integer(UCBadmit$male2)

现在模型需要检测 male2 何时消失(-1)，然后计算未知状态下的混合物。

m_mix <- ulam(
    alist(
        admit|male2==-1 ~ custom( log_mix( 
            phi_male , 
            binomial_lpmf(admit|applications,p_m1) , 
            binomial_lpmf(admit|applications,p_m0) ) ),
        admit|male2>-1 ~ binomial( applications , p ),
        logit(p) <- a[dept] + b*male2,
        logit(p_m1) <- a[dept] + b*1,
        logit(p_m0) <- a[dept] + b*0,
        male2|male2>-1 ~ bernoulli( phi_male ),
        phi_male ~ beta(2,2),
        a[dept] ~ normal(0,4),
        b ~ normal(0,1)
    ),
    data=UCBadmit )

请注意将phi_male添加到未知状态的平均值。

连续缺失数据估算
原则上，缺失的实值数据的估算是很容易的:只需用一个参数替换掉每个缺失的值。实际上，这涉及到一堆烦人的统计工作。ulam有一个名为merge_missing的宏来简化它。

UCBadmit$x <- rnorm(12)
UCBadmit$x[1:2] <- NA
m_miss <- ulam(
    alist(
        admit ~ binomial(applications,p),
        logit(p) <- a + b*male + bx*x_merge,
        x_merge ~ normal( 0 , 1 ),
        x_merge <- merge_missing( x , x_impute ),
        a ~ normal(0,4),
        b ~ normal(0,1),
        bx ~ normal(0,1)
    ),
    data=UCBadmit )

merge_missing所做的是找到x中的NA值(无论哪个符号是第一个参数)，构建一个名为x_impute的参数向量(无论第二个参数的名称是什么)，长度合适，然后在合适的位置组合一个包含这两个值的向量x_merge。然后你可以给这个向量分配一个先验，并像往常一样在线性模型中使用它。
在Stan模型运行时，使用自定义函数块完成合并。查看Stan code stancode(m_miss)了解所有的细节。
merge missing是一个宏的例子，它是ulam使用函数名来触发sp的一种方式。在这种情况下，merge_missing在Stan模型中插入一个函数，并构建必要的索引来定位运行时丢失的值。一旦系统完成，宏将在以后得到完整的文档。

高斯过程
一个简单的高斯过程，就像书中第13章中的海洋岛屿的例子，是这样做的:

data(Kline2)
d <- Kline2
data(islandsDistMatrix)
d$society <- 1:10
dat <- list(
    y=d$total_tools,
    society=d$society,
    log_pop = log(d$population),
    Dmat=islandsDistMatrix
)

m_GP1 <- ulam(
    alist(
        y ~ poisson( mu ),
        log(mu) <- a + aj[society] + b*log_pop,
        a ~ normal(0,10),
        b ~ normal(0,1),
        vector[10]: aj ~ multi_normal( 0 , SIGMA ),
        matrix[10,10]: SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),
        etasq ~ exponential(1),
        rhosq ~ exponential(1)
    ),
    data=dat )

这只是一个普通的多种截距模型，但所有的10个截距都来自一个单一的高斯分布。协方差矩阵SIGMA 定义为一般的 L2 -norm。同样，cov_GPL2是一个宏，它在Stan代码中插入一个函数，在模型运行时计算协方差矩阵。
更好看的高斯过程需要不同的参数化。这些也可以建造。这里有一个使用151种灵长类动物和一个系统进化距离矩阵的例子。首先，准备数据:

data(Primates301)
data(Primates301_distance_matrix)
d <- Primates301
d$name <- as.character(d$name)
dstan <- d[ complete.cases( d$social_learning, d$research_effort , d$body , d$brain ) , ]
# prune distance matrix to spp in dstan
spp_obs <- dstan$name
y <- Primates301_distance_matrix
y2 <- y[ spp_obs , spp_obs ]
# scale distances
y3 <- y2/max(y2)

模型是非中心 L2 -norm 高斯过程:

m_GP2 <- ulam(
    alist(
        social_learning ~ poisson( lambda ),
        log(lambda) <- a + g[spp_id] + b_ef*log_research_effort + b_body*log_body + b_eq*log_brain,
        a ~ normal(0,1),
        vector[N_spp]: g <<- L_SIGMA * eta,
        vector[N_spp]: eta ~ normal( 0 , 1 ),
        matrix[N_spp,N_spp]: L_SIGMA <<- cholesky_decompose( SIGMA ),
        matrix[N_spp,N_spp]: SIGMA <- cov_GPL2( Dmat , etasq , rhosq , 0.01 ),
        b_body ~ normal(0,1),
        b_eq ~ normal(0,1),
        b_ef ~ normal(1,1),
        etasq ~ exponential(1),
        rhosq ~ exponential(1)
    ),
    data=list(
        N_spp = nrow(dstan),
        social_learning = dstan$social_learning,
        spp_id = 1:nrow(dstan),
        log_research_effort = log(dstan$research_effort),
        log_body = log(dstan$body),
        log_brain = log(dstan$brain),
        Dmat = y3
    ) , 
    control=list(max_treedepth=15,adapt_delta=0.95) ,
    sample=FALSE )

这个模型采样不快，所以我设置了sample=FALSE。您仍然可以使用stancode(m_GP2)检查Stan代码。
注意协方差 SIGMA 的构建和之前一样，但是我们马上把它分解成一个Cholesky因子然后通过矩阵乘法来构建变化的截距 g 。操作符 <<- 告诉ulam不要循环，而是直接赋值。g <<- L_SIGMA *做了正确的线性代数运算。

过程进行中的工作
ulam 仍在快速发展中。但它将是教科书第二版的核心，允许更精确、更透明的模型定义，减少对实际情况的隐藏。
map2stan语法和特性
旧的map2stan函数对将要看到的公式进行了更强的假设。这允许提供一些额外的自动化，因此它有一些特殊的语法。相反，ulam通过其宏库支持这些特性。

非中心参数化
这是一个非中心参数化，在线性模型之前移动尺度参数的变化效果，这通常是更有效的采样:

f4u <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + zaj[group]*sigma_group[1] + 
         (b + zbj[group]*sigma_group[2])*x,
    c(zaj,zbj)[group] ~ dmvnorm( 0 , Rho_group ),
    a ~ dnorm( 0 , 10 ),
    b ~ dnorm( 0 , 1 ),
    sigma ~ dcauchy( 0 , 1 ),
    sigma_group ~ dcauchy( 0 , 1 ),
    Rho_group ~ dlkjcorr(2)
)

书的第13章提供了很多关于这个问题的细节。
我们可以进一步采用这个策略，并从先验中删除相关矩阵Rho_group。map2stan通过dmvnormNC密度简化了这种形式，它使用相关矩阵的内部Cholesky分解来构建不同的效果。这是之前的变斜率模型，现在用非中心符号表示:

f4nc <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + aj[group] + (b + bj[group])*x,
    c(aj,bj)[group] ~ dmvnormNC( sigma_group , Rho_group ),
    a ~ dnorm( 0 , 10 ),
    b ~ dnorm( 0 , 1 ),
    sigma ~ dcauchy( 0 , 1 ),
    sigma_group ~ dcauchy( 0 , 1 ),
    Rho_group ~ dlkjcorr(2)
)

在内部，使用Cholesky因子L_Rho_group执行抽样。它将出现在返回的示例中，以及由它构造的Rho_group中。

半自动贝叶斯设算（imputation）
编码简单的贝叶斯估算是可能的。例如，让我们模拟一个缺少预测值的简单回归:

N <- 100
N_miss <- 10
x <- rnorm( N )
y <- rnorm( N , 2*x , 1 )
x[ sample(1:N,size=N_miss) ] <- NA

移除10个x值。然后map2stan公式列表定义了一个x的分布:

f5 <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + b*x,
    x ~ dnorm( mu_x, sigma_x ),
    a ~ dnorm( 0 , 100 ),
    b ~ dnorm( 0  , 10 ),
    mu_x ~ dnorm( 0 , 100 ),
    sigma_x ~ dcauchy(0,2),
    sigma ~ dcauchy(0,2)
)
m5 <- map2stan( f5 , data=list(y=y,x=x) )

map2stan所做的是注意缺失值，查看分配给缺失值变量的分布，构建Stan代码，该代码在回归中使用观察到的和估计的x值的混合。有关实现的详细信息，请参见stancode(m5)。

二进制离散缺失值的半自动边缘化
缺少值的二进制(0/1)变量是一个特殊的障碍，因为Stan不能对离散参数进行采样。因此，map2stan可以对缺失的二进制值进行平均(边缘化)，而不是输入缺失的二进制值。在上面的例子中，当map2stan检测到一个预测变量中缺失的值时，它会试图找到包含这些值的变量的分布。如果这个变量是二进制的(0/1)，那么它将构建一个混合模型，其中每一项都是log-likelihood，条件是变量采用特定的0/1值组合。
按照上一节的例子，我们可以在一个二元预测器中模拟缺失:

N <- 100
N_miss <- 10
x <- rbinom( N , size=1 , prob=0.5 )
y <- rnorm( N , 2*x , 1 )
x[ sample(1:N,size=N_miss) ] <- NA

模型定义类似于前面的定义，但在为定义x分布的超参数指定约束时也需要注意:

f6 <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + b*x,
    x ~ bernoulli( phi ),
    a ~ dnorm( 0 , 100 ),
    b ~ dnorm( 0  , 10 ),
    phi ~ beta( 1 , 1 ),
    sigma ~ dcauchy(0,2)
)
m6 <- map2stan( f6 , data=list(y=y,x=x) , constraints=list(phi="lower=0,upper=1") )

理论上，该算法适用于任何数量的丢失值的二进制预测器。例如，有两个预测器，每个都有缺失:

N <- 100
N_miss <- 10
x1 <- rbinom( N , size=1 , prob=0.5 )
x2 <- rbinom( N , size=1 , prob=0.1 )
y <- rnorm( N , 2*x1 - x2  , 1 )
x1[ sample(1:N,size=N_miss) ] <- NA
x2[ sample(1:N,size=N_miss) ] <- NA
f7 <- alist(
    y ~ dnorm( mu , sigma ),
    mu <- a + b1*x1 + b2*x2,
    x1 ~ bernoulli( phi1 ),
    x2 ~ bernoulli( phi2 ),
    a ~ dnorm( 0 , 100 ),
    c(b1,b2) ~ dnorm( 0  , 10 ),
    phi1 ~ beta( 1 , 1 ),
    phi2 ~ beta( 1 , 1 ),
    sigma ~ dcauchy(0,2)
)
m7 <- map2stan( f7 , data=list(y=y,x1=x1,x2=x2) , 
      constraints=list(phi1="lower=0,upper=1",phi2="lower=0,upper=1") )

虽然二元预测器的未观察值通常不重要，但它们可以从后验分布计算出来。添加参数do_discrete te_imputation=TRUE指示map2stan自动执行这些计算。例子:

m6 <- map2stan( f6 , data=list(y=y,x=x) , constraints=list(phi="lower=0,upper=1") ,
      do_discrete_imputation=TRUE )
precis( m6 , depth=2 )

输出包含每种情况的样本，假设x取值为1。
该算法通过构造一个混合项列表来计算每个观察到的y值的概率。在最简单的情况下，只有一个预测值缺失，隐含的混合似然包含两个项:

Pr(y[i]) = Pr(x[i]=1)Pr(y[i]|x[i]=1) + Pr(x[i]=0)Pr(y[i]|x[i]=0)

在我们上面的例子模型m6的参数中，这是:

Pr(y[i]) = phi*N(y[i]|a+b,sigma) + (1-phi)*N(y[i]|a,sigma)

现在循环遍历 i 并为每个例子计算上面的值就很简单了。同理，x[i]==1的后验概率为:

Pr(x[i]==1|y[i]) = phi*N(y[i]|a+b,sigma) / Pr(y[i])

当只有一个预测器丢失时，这很简单。如果有两个或更多呢?在这种情况下，失踪的所有可能组合都必须加以说明。例如，假设有两个预测因子，x1和x2，在情形i中均不存在，则隐含的混合概率为:

Pr(y[i]) = Pr(x1=1)Pr(x2=1)*Pr(y[i]|x1=1,x2=1) + Pr(x1=1)Pr(x2=0)Pr(y[i]|x1=1,x2=0) + Pr(x1=0)Pr(x2=1)Pr(y[i]|x1=0,x2=1) + Pr(x1=0)Pr(x2=0)Pr(y[i]|x1=0,x2=0)

有四种未观测值的组合，因此混合可能性有四项。当x2被观察到时，我们可以将观察到的值代入上述，那么混合物就很容易简化为我们之前的两项似然:

Pr(y[i]|x2[i]==1) = Pr(x1=1)Pr(x2=1)Pr(y[i]|x1=1,x2=1) + Pr(x1=1)Pr(x2=0)Pr(y[i]|x1=1,x2=1) + Pr(x1=0)Pr(x2=1)Pr(y[i]|x1=0,x2=1) + Pr(x1=0)Pr(x2=0)Pr(y[i]|x1=0,x2=1)
                  = [Pr(x1=1)Pr(x2=1)+Pr(x1=1)Pr(x2=0)]Pr(y[i]|x1=1,x2=1) 
                    + [Pr(x1=0)Pr(x2=1)+Pr(x1=0)Pr(x2=0)]Pr(y[i]|x1=0,x2=1)
                  = Pr(x1=1)Pr(y[i]|x1=1,x2=1) + Pr(x1=0)Pr(y[i]|x1=0,x2=1)

这意味着，如果我们对情况 i 进行循环，并将任何观察到的值插入到一般的混合似然中，我们就可以为每种情况i的特定缺失组合计算相关的混合。这就是map2stan所做的。一般的混合项可以通过算法生成。下面的代码为n个丢失的二进制变量生成一个术语矩阵。

ncombinations <- 2^n
d <- matrix(NA,nrow=ncombinations,ncol=n)
for ( col_var in 1:n ) 
    d[,col_var] <- rep( 0:1 , each=2^(col_var-1) , length.out=ncombinations )

d的行包含terms，列包含变量，每一列的值是每一个变量对应的值。该算法为矩阵中的每一行建立一个线性模型，将混合似然作为这些行的和，并对观测值进行适当的替换。为了精确起见，所有的计算都是在对数刻度上进行的。

高斯过程
基本的高斯过程可以用GPL2分布标签指定。这意味着一个多元高斯与协方差矩阵定义的普通L2范数距离函数:

k(i,j) = eta^2 * exp( -rho^2 * D(i,j)^2 ) + ifelse(i==j,sigma^2,0)

其中D是两两距离的矩阵。例如，在一个空间自相关模型中使用这个约定:

library(rethinking)
data(Kline2)
d <- Kline2
data(islandsDistMatrix)
d$society <- 1:10
mGP <- map2stan(
    alist(
        total_tools ~ dpois( mu ),
        log(mu) <- a + aj[society],
        a ~ dnorm(0,10),
        aj[society] ~ GPL2( Dmat , etasq , rhosq , 0.01 ),
        etasq ~ dcauchy(0,1),
        rhosq ~ dcauchy(0,1)
    ),
    data=list(
        total_tools=d$total_tools,
        society=d$society,
        Dmat=islandsDistMatrix),
    constraints=list(
        etasq="lower=0",
        rhosq="lower=0"
    ),
    warmup=1000 , iter=5000 , chains=4 )

请注意使用constraints列表将custom参数constraints传递给Stan。书中更详细地探讨了这个例子。

信息准则
map和map2stan都提供了DIC和WAIC。大多数情况下是这样的。实际上，这两种工具都足够灵活，您可以指定既不能正确计算DIC也不能正确计算WAIC的模型。但是对于普通的GLMs和GLMMs，它是有效的。看到R帮助了吗?一个方便的比较函数总结了信息标准比较，包括WAIC的标准误差。
ulam支持使用可选的log_lik=TRUE参数进行WAIC计算，该参数返回loo包所需的对数可能性向量。
集成为一个模型集成(eac)计算链路和sim输出

文章仅作参考。如果有一些翻译问题的话，请在下方留言，另外建议还是直接看原网站以及翻阅原书籍。