使用brms做贝叶斯分析遇Stan采样初始化失败错误求助

问题背景 #

使用brms库拟合非线性贝叶斯模型时，出现初始化失败、采样无法执行的错误。

代码实现 #

library(brms)

mymodel <- brms::bf(tau ~ a + b*(gmm^c),
                    a ~ 1,
                    b ~ 1,
                    c ~ 1,
                    nl = TRUE)

myprior <- c(prior(uniform(10^-6, 3000), class = 'b', nlpar = 'a', lb = 10^-6, ub = 3000),
             prior(uniform(10^-6, 1000), class = 'b', nlpar = 'b', lb = 10^-6, ub = 1000),
             prior(uniform(10^-6, 1), class = 'b', nlpar = 'c', lb = 10^-6, ub = 1))

myprior$prior

myfit <- brm(formula = mymodel,
                data = dat,
              family = brmsfamily('Gamma'),
               prior = myprior,
             control = list(adapt_delta = 0.95),
                iter = 10000,
              warmup = 5000,
                thin = 5,
              chains = 1,
               cores = 6,
             refresh = 0)

数据示例 #

> print(dat, n=45)
# A tibble: 14,456 × 4
   fase  ts          gmm   tau
   <fct> <fct>     <dbl> <dbl>
 1 1     30     0.00316   19.3
 2 1     30     0.0153    45.0
 3 1     30     0.0336    55.9
 4 1     30     0.0654    60.5
 5 1     30     0.124     61.7
 6 1     30     0.234     62.4
 7 1     30     0.440     64.0
 8 1     30     0.828     65.2
 9 1     30     1.56      67.3
10 1     30     2.93      71.6
11 1     30     5.51      78.7
12 1     30    10.4       84.6
13 1     30    19.5       91.1
14 1     30    36.8       99.7
15 1     30    69.2      113. 
16 1     30     0.00558   13.8
17 1     30     0.0172    26.4
18 1     30     0.0342    31.5
19 1     30     0.0657    34.3
20 1     30     0.124     34.9
21 1     30     0.234     35.7
22 1     30     0.440     36.7
23 1     30     0.828     38.4
24 1     30     1.56      41.4
25 1     30     2.93      44.0
26 1     30     5.51      49.1
27 1     30    10.4       53.3
28 1     30    19.5       58.2
29 1     30    36.8       65.3
30 1     30    69.2       77.8
31 1     30     0.000369  24.4
32 1     30     0.00252   91.4
33 1     30     0.0194   175. 
34 1     30     0.0599   229  
35 1     30     0.124    240. 
36 1     30     0.232    244. 
37 1     30     0.439    254. 
38 1     30     0.829    256. 
39 1     30     1.56     258. 
40 1     30     2.93     273. 
41 1     30     5.52     289. 
42 1     30    10.4      311. 
43 1     30    19.5      333. 
44 1     30    36.8      346. 
45 1     30    69.2      349. 

错误信息 #

• Error : Initialization failed.
• error occurred during calling the sampler; sampling not done

可能原因分析 #

• 非线性项数值不稳定：模型中的gmm^c项，当gmm取值极小（如0.000369）且c接近0时，计算易出现数值溢出或下溢，导致初始似然计算异常，触发初始化失败。
• 先验与参数范围冲突：c的先验范围包含接近0的值，会放大gmm^c的数值不稳定问题；同时Gamma分布要求响应变量tau严格为正，初始参数组合可能生成非正预测值，违背分布假设。
• 变量尺度差异过大：gmm取值跨度达5个数量级，tau跨度也超过20倍，这种尺度差异会导致模型Hessian矩阵条件数变差，MCMC采样器难以找到有效初始值。
• 单链初始化风险：仅设置1条链，初始化失败时无其他链参考，建议尝试多链（如4链），同时通过control = list(init_r = 0.1)缩小初始值范围，增加初始化成功率。
• 模型参数化方式可优化：将gmm^c重新参数化为exp(c * log(gmm))，能提升数值计算稳定性；或对gmm做对数变换后再建模，降低尺度差异影响。

内容的提问来源于stack exchange，提问作者Gabriel Fonseca