Bernoulli误差与状态空间模型的rstan估计

dxxyhpgq  于 2023-01-03  发布在  其他
关注(0)|答案(1)|浏览(149)

我尝试使用状态空间模型来估计每次判断的准确性。R代码如下。使用MacOS 10.14.6,rstan 2.21.2。

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

我明白了围绕“法官[i] ~伯努利(p[idx_t[i]]);“,可是我一点主意都没有。
输入文件如下所示。

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
oknrviil

oknrviil1#

需要执行以下流程。

dt_idx_time = data.table(
  time_id = sort(unique(data$time)),
  time_id_rev = 1 : length(sort(unique(data$time)))
)
vec_rev_id = dt_idx_time$time_id_rev[match(data$time, dt_idx_time$time_id)]

dt002 = data.table(
  data,
  time_rev = vec_rev_id
)

data_new = dt002

相关问题