问题描述
我是 Stan 和 rstan 的新手。
我最近在研究马尔可夫链蒙特卡罗 (MCMC) 时可能会发现一个奇怪的问题。简而言之,例如,数据有 10 个观察值,比如 ID 1 到 10。现在,我通过在原始第一行和第二行之间移动第 10 行来排列它,比如 ID 1、10 和2 到 9。即使我固定相同的随机种子,两种不同的数据场景也会给出不同的估计。
为了以更简单的方式说明问题,我编写了以下 R 脚本。
##TEST 01
# generate data
N <- 100
set.seed(123)
Y <- rnorm(N,1.6,0.2)
stan_code1 <- "
data {
int <lower=0> N; //number of data
real Y[N]; //data in an (C++) array
}
parameters {
real mu; //mean parameter of a normal distribution
real <lower=0> sigma; //standard deviation parameter of a normal distribution
}
model {
//prior distributions for parameters
mu ~ normal(1.7,0.3);
sigma ~ cauchy(0,1);
//likelihood of Y given parameters
for (i in 1:N) {
Y[i] ~ normal(mu,sigma);
}
}
"
# compile model
library(rstan)
model1 <- stan_model(model_code = stan_code1) #usually,take half a minute to run
# pass data to stan and run model
set.seed(123)
fit <- sampling(model1,list(N=N,Y=Y),iter=200,chains=4)
print(fit)
# mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
# mu 1.62 0.00 0.02 1.58 1.61 1.62 1.63 1.66 473 1.00
# sigma 0.18 0.00 0.01 0.16 0.18 0.18 0.19 0.21 141 1.02
# lp__ 117.84 0.07 0.85 115.77 117.37 118.07 118.51 118.78 169 1.01
Yp <- Y[c(1,100,2:99)]
set.seed(123)
fit2 <- sampling(model1,Y=Yp),chains=4)
print(fit2)
# mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
# mu 1.62 0.00 0.02 1.59 1.61 1.62 1.63 1.66 480 0.99
# sigma 0.18 0.00 0.01 0.16 0.18 0.18 0.19 0.21 139 1.02
# lp__ 117.79 0.09 0.95 115.72 117.35 118.05 118.49 118.77 124 1.01
从上面的简单案例中我们可以看出,fit 和 fit2 的两个结果是不同的。 而且,更奇怪的是,如果我在先验之前写可能性(以前,先验写在可能性之前 em>) 在代码文件中,相同的随机种子和相同的数据仍然会给出不同的估计。
##TEST 01'
# generate data
#N <- 100
set.seed(123)
Y <- rnorm(N,0.2)
stan_code11 <- "
data {
int <lower=0> N; //number of data
real Y[N]; //data in an (C++) array
}
parameters {
real mu; //mean parameter of a normal distribution
real <lower=0> sigma; //standard deviation parameter of a normal distribution
}
model {
//likelihood of Y given parameters
for (i in 1:N) {
Y[i] ~ normal(mu,sigma);
}
//prior distributions for parameters
mu ~ normal(1.7,1);
}
"
# compile model
#library(rstan)
model11 <- stan_model(model_code = stan_code11) #usually,take half a minute to run
# pass data to stan and run model
set.seed(123)
fit11 <- sampling(model11,chains=4)
print(fit11)
# mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
# mu 1.62 0.00 0.02 1.58 1.61 1.62 1.63 1.66 455 0.99
# sigma 0.19 0.00 0.01 0.16 0.18 0.18 0.20 0.21 94 1.04
# lp__ 117.68 0.08 0.93 115.24 117.18 117.90 118.45 118.77 149 1.01
##TEST01 was
# mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
# mu 1.62 0.00 0.02 1.58 1.61 1.62 1.63 1.66 473 1.00
# sigma 0.18 0.00 0.01 0.16 0.18 0.18 0.19 0.21 141 1.02
# lp__ 117.84 0.07 0.85 115.77 117.37 118.07 118.51 118.78 169 1.01
解决方法
Stan 不使用与 R 相同的伪随机数生成器。因此,调用 set.seed(123)
只会使 Y
可重复,而不会使 MCMC 采样可重复。为了完成后者,您需要将一个整数作为 seed
参数传递给 rstan 包中的 stan
(或 sampling
)函数,例如
sampling(model11,list(N = N,Y = Y),seed = 1234)
。
即使那样,我也可以想象,由于浮点原因,置换观察结果可能会导致后验分布的绘制实现不同。但这一切都无关紧要(除非您进行的迭代次数太少或收到警告消息),因为即使后验分布的有限实现集是随机不同的数字,后验分布也是相同的。