Блог

Моя цель состоит в том, чтобы в основном перенести этот код в R. Вся предварительная обработка наборов данных wrt уже выполнена, однако теперь я застрял в написании файла «модель». В качестве первой попытки и для большей ясности я написал код, который показан ниже на языке R.

Что я хочу сделать, так это запустить MCMC, чтобы получить оценку параметра R_t с учетом ежедневных отчетных данных по итальянской стране. Основными шагами, которые были предприняты, являются:

1. Пример параметра массива, а именно журнала(R_t), из гауссовского распределения RW

 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)
 

2. Вычислите некоторые величины, которые являются функцией этих выборочных параметров, а именно количество инфекций. Эта зависимость довольно проста, так как это линейная алгебра:

 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,])
}

 

3. Сверните массив, который я только что вычислил, с распределением задержки (начало задержка сообщения), наконец, масштабируя его по переменной экспозиции:

 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) 
 

4. Вычислите вероятность, которая пропорциональна вероятности того, что наблюдался определенный набор данных (т. Е. ежедневно подтвержденные случаи), путем выборки вышеупомянутого параметра log(R_t), из которого вычисляется переменная positive.

 likelihood <- dnbinom(round(Italian_data$daily_confirmed), mu = positive, size = 1/6) 
 

Наконец, мы подошли к моему файлу модели ОШИБОК:

 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 самостоятельно в соответствии с приведенной выше вероятностью, но, очевидно, я бы гораздо больше предпочел использовать уже протестированную платформу, которая содержит ОШИБКИ/ОШИБКИ, вот почему я прошу о помощи.