使用rstan实现伯努利状态空间模型时遇索引越界错误求助

问题背景 #

在MacOS 10.14.6环境下，使用rstan 2.21.2尝试用状态空间模型估计各时间点的判断准确率，运行代码时出现索引越界错误，定位到Stan代码中judge[i] ~ bernoulli(p[idx_t[i]])行，但不知如何解决。

原R与Stan代码 #

library(rstan)
rstan_options(auto_write=TRUE)
options(mc.cores=parallel::detectCores())

################################
#Stan
stan_model <- "
data {
  int n; //number of data
  int n_t; //number of time index
  int idx_t[n]; //index for time
  int judge[n]; //judgment
}

parameters {
  real<lower = 0, upper = 1> mu[n_t]; // mean of correct judgment
  real<lower = 0> s_mu; 
}

transformed parameters {
  real<lower = 0, upper = 1> p[n_t];
  for(i in 1 : n_t) {
    p[i] = mu[i];
  }
}


model {
  for ( i in 2 : n_t) {
    mu[i] ~ normal(mu[i - 1], s_mu);
  }
  
  for(i in 1 : n) {
    judge[i] ~ bernoulli(p[idx_t[i]]);
  }
}
"
##########################################
##########################################
d <- read.csv('./input100.csv')


d$Time_100 = d$Time * 100
d$Time_100_new = as.integer(d$Time_100)

data = data.frame(d$Time_100_new,d$judgement)
colnames(data) <- c("time","ratio")
#data

qwe = data[order(data$time),]
data = qwe

d_input <- list(n = length(data$ratio), n_t = length(sort(unique(data$time))), idx_t = data$time, judge = data$ratio)
fit <- stan(model_code = stan_model,
               data = d_input)

错误信息 #

SAMPLING FOR MODEL '2cf9ed042422a42c681f58525db6a2eb' NOW (CHAIN 1).
Chain 1: Unrecoverable error evaluating the log probability at the initial value.
Chain 1: Exception: []: accessing element out of range. index 0 out of range; expecting index to be between 1 and 13; index position = 1p  (in 'model11625186b3cc5_2cf9ed042422a42c681f58525db6a2eb' at line 28)

[1] "Error in sampler$call_sampler(args_list[[i]]) : "                                                                                                                                                        
[2] "  Exception: []: accessing element out of range. index 0 out of range; expecting index to be between 1 and 13; index position = 1p  (in 'model11625186b3cc5_2cf9ed042422a42c681f58525db6a2eb' at line 28)"
error occurred during calling the sampler; sampling not done

输入文件内容 #

Time,judgement
0,1
0.04,1
0.07,0
0.08,1
0.08,1
0.08,1
0.08,0
0.08,0
0.08,0
0.08,0
0.09,1
0.09,1
0.09,1
0.09,1
0.09,1
0.09,0
0.09,0
0.09,0
0.09,0
0.1,1
0.1,1
0.1,1
0.1,1
0.1,1
0.1,0
0.1,0
0.1,0
0.1,0
0.11,1
0.11,1
0.11,1
0.11,1
0.11,1
0.11,1
0.11,1
0.11,0
0.11,0
0.11,0
0.11,0
0.11,0
0.11,0
0.11,0
0.12,1
0.12,1
0.12,1
0.12,1
0.12,1
0.12,1
0.12,0
0.12,0
0.12,0
0.12,0
0.13,1
0.13,1
0.13,1
0.13,1
0.13,0
0.13,0
0.13,0
0.13,0
0.13,0
0.13,0
0.14,1
0.14,1
0.14,1
0.14,0
0.14,0
0.14,0
0.14,0
0.14,0
0.15,1
0.15,1
0.15,1
0.15,1
0.15,1
0.15,1
0.15,0
0.15,0
0.15,0
0.15,0
0.16,1
0.16,1
0.16,1
0.16,1
0.16,1
0.16,0
0.16,0
0.16,0
0.16,0
0.17,1
0.17,1
0.17,1
0.17,1
0.17,1
0.17,1
0.17,1
0.17,1
0.17,1
0.17,0

问题原因与解决方法 #

原因分析 #

Stan采用1-based索引（数组从1开始计数），但输入数据中Time为0时，转换后的Time_100_new是0，导致idx_t数组中存在0值。当Stan代码尝试用idx_t[i]作为索引访问p数组时，0超出了p数组的有效索引范围（1到13），触发越界错误。

解决步骤 #

修改R代码中生成时间索引的部分，将所有时间值加1，确保索引从1开始：

# 原代码
d$Time_100_new = as.integer(d$Time_100)
# 修改为
d$Time_100_new = as.integer(d$Time_100) + 1

n_t的计算length(sort(unique(data$time)))会自动更新为14（原13个时间点+1偏移），无需额外修改。

验证修改 #

运行修改后的代码，idx_t中的值将从1开始，对应p数组的有效索引范围，可避免越界错误。

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