rstan MCMC：不同的数据顺序导致不同的结果，为什么？

如何解决rstan MCMC：不同的数据顺序导致不同的结果，为什么？

我是 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)。

即使那样，我也可以想象，由于浮点原因，置换观察结果可能会导致后验分布的绘制实现不同。但这一切都无关紧要（除非您进行的迭代次数太少或收到警告消息），因为即使后验分布的有限实现集是随机不同的数字，后验分布也是相同的。