0

我的目标是基本上将此代码迁移到 R。所有预处理 wrt 数据集都已经完成，但是现在我被困在编写“模型”文件中。作为第一次尝试，为了清楚起见，我用 R 语言编写了如下所示的代码。

鉴于意大利国家的每日报告数据，我想要做的是运行 MCMC 来估计参数 R_t。已采取的主要步骤是：

  1. 从高斯 RW 分布中采样一个数组参数，即 log(R_t)
Gauss_RandomWalk <- function(N, x0, mu, variance) {   
z <- cumsum(rnorm(n=N, mean=mu, sd=sqrt(variance)))   
t <- 1:N   
x <- (x0 + t*mu + z)   
return(x) 
}

log_R_t <- Gauss_RandomWalk(tot_dates, 0., 0., 0.035**2)
R_t_candidate <- exp(log_R_t)
  1. 计算一些数量，这些数量是该采样参数的函数，即感染数量。这种依赖关系非常简单，因为它是线性代数：
infections <- rep(0. , tot_dates)
infections[1] <- exp(seed)
for (t in 2:tot_dates){
  infections[t] <- sum(R_t_candidate * infections * gt_to_convolution[t-1,])
}
  1. 将我刚刚计算的数组与延迟分布（开始+报告延迟）进行卷积，最后通过曝光变量重新调整它：
test_adjusted_positive <- convolve(infections, delay_distribution_df$density, type = "open")
test_adjusted_positive <- test_adjusted_positive[1:tot_dates]
positive <- round(test_adjusted_positive*exposure)
  1. 通过对上述 log(R_t) 参数进行采样，计算出与观察到某组数据（即每日确诊病例）的概率成正比的可能性，从中计算变量正数。
likelihood <- dnbinom(round(Italian_data$daily_confirmed), mu = positive, size = 1/6)

最后，我们来到我的BUGS模型文件：

model {

  #priors as a Gaussian RW
  log_rt[1] ~ dnorm(0, 0.035)
  log_rt[2] ~ dnorm(0, 0.035)
  for (t in 3:tot_dates) {
    log_rt[t] ~ dnorm(log_rt[t-1] + log_rt[t-2], 0.035)
    R_t_candidate[t] <- exp(log_rt[t])
  }
  
    # data likelihood
    for (t in 2:tot_dates) {
        infections[t] <- sum(R_t_candidate * infections * gt_to_convolution[t-1,])
    }
    
    test_adjusted_positive <- convolve(infections, delay_distribution)
  test_adjusted_positive <- test_adjusted_positive[1:tot_dates]
  positive <- test_adjusted_positive*exposure

  for (t in 2:tot_dates) {
        confirmed[t] ~ dnbinom( obs[t], positive[t], 1/6) 
    }
}

其中gt_to_convolution是一个常数矩阵，tot_dates是一个常数值，exposure是一个常数数组。

尝试通过以下方式编译它时：

data <- NULL
data$obs <- round(Italian_data$daily_confirmed)
data$tot_dates <- n_days
data$delay_distribution <- delay_distribution_df$density
data$exposure <- exposure
data$gt_to_convolution <- gt_to_convolution

inits <- NULL
inits$log_rt  <- rep(0, tot_dates)

library (rjags)
library (coda)

set.seed(1995)

model <- "MyModel.bug"
jm <- jags.model(model , data, inits)

它引发以下引发错误：

Compiling model graph
   Resolving undeclared variables
   Allocating nodes
Deleting model

Error in jags.model(model, data, inits) : RUNTIME ERROR:
Compilation error on line 19.
Possible directed cycle involving test_adjusted_positive

因此，我什至无法对其进行一点调试，即使我很确定总体上存在一些问题，但我无法弄清楚是什么以及为什么。

在这一点上，我认为最好的选择是根据上面的可能性自己实现一个 Metropolis 算法，但显然，我更喜欢使用已经测试过的框架，即 BUGS/JAGS，这就是为什么我我正在寻求帮助。

4

0 回答 0