Maybe a silly question. I can call the function using the "odin.dust::odin_dust" statement in the R software, as shown in Figure 1. However, when running "mcstate::particle_filter$new" I encountered an error: "Error in model$public_methods$time_type() : attempt to apply non-function" as shown in Figure 2-3. I also tried using other computers, but the errors seemed to be more frequent, possibly due to the lack of C++ environment configuration, as shown in Figure 4-5.
But if I want to change the initial values of the SIR model or the dynamics of the transmission model (such as SEIR), it seems necessary to write the model myself and call it through "odin.dust::odin_dust". My code is provided below, and I would appreciate your guidance on how to make modifications.
Maybe a silly question. I can call the function using the "odin.dust::odin_dust" statement in the R software, as shown in Figure 1. However, when running "mcstate::particle_filter$new" I encountered an error: "Error in model$public_methods$time_type() : attempt to apply non-function" as shown in Figure 2-3. I also tried using other computers, but the errors seemed to be more frequent, possibly due to the lack of C++ environment configuration, as shown in Figure 4-5.
But if I want to change the initial values of the SIR model or the dynamics of the transmission model (such as SEIR), it seems necessary to write the model myself and call it through "odin.dust::odin_dust". My code is provided below, and I would appreciate your guidance on how to make modifications.
` rm(list = ls()) library(mcstate) library(odin) library("odin.dust")
incidence <- read.csv(system.file("sir_incidence.csv", package = "mcstate")) incidence plot(cases ~ day, incidence, type = "o", xlab = "Day", ylab = "New cases", pch = 19)
history <- readRDS("C:/sir_true_history.rds") history <- t(drop(history)) colnames(history) <- c("S", "I", "R", "cases") history <- data.frame(t = seq_len(nrow(history)) - 1L, history) history matplot(history$t, history[c("S", "I", "R")], type = "l", lty = 1, xlab = "Day", ylab = "Variable")
sir <- odin.dust::odin_dust({ update(S) <- S - n_SI update(I) <- I + n_SI - n_IR update(R) <- R + n_IR
n_IR <- rbinom(I, 1 - exp(-beta * I / (S + I + R))) n_SI <- rbinom(S, 1 - exp(-gamma))
initial(S) <- S_ini initial(I) <- I_ini initial(R) <- 0
S_ini <- user(1000) I_ini <- user(10) beta <- user(0.2) gamma <- user(0.1) })
dt <- 0.25 data <- mcstate::particle_filter_data(incidence, time = "day", rate = 1 / dt)
compare <- function(state, observed, pars = NULL) { if (is.na(observed$cases)) { return(NULL) } exp_noise <- 1e6 incidence_modelled <- state[1, , drop = TRUE] incidence_observed <- observed$cases lambda <- incidence_modelled + rexp(n = length(incidence_modelled), rate = exp_noise) dpois(x = incidence_observed, lambda = lambda, log = TRUE) }
index <- function(info) { list(run = 5L, state = 1:5) }
proposal_kernel <- rbind(c(0.00057, 0.00034), c(0.00034, 0.00026)) pars <- mcstate::pmcmc_parameters$new( list(mcstate::pmcmc_parameter("beta", 0.2, min = 0, max = 1, prior = function(p) log(1e-10)), mcstate::pmcmc_parameter("gamma", 0.1, min = 0, max = 1, prior = function(p) log(1e-10))), proposal = proposal_kernel) pars
n_steps <- 500 n_particles <- 100 n_threads <- dust::dust_openmp_threads()
p1 <- mcstate::particle_filter$new(data, sir, n_particles, compare, index = index, n_threads = n_threads) control1 <- mcstate::pmcmc_control(n_steps, save_trajectories = TRUE, save_state = TRUE, save_restart = 40) res1 <- mcstate::pmcmc(pars, p1, control = control1)
t <- 0:100 matplot(t, t(res1$trajectories$state[2, , ]), type = "l", lty = 1, col = "https://github.com/mrc-ide/odin/issues/11", xlab = "Day", ylab = "Infected") points(I ~ t, history, pch = 19, col = "blue")
matplot(t, t(res1$trajectories$state[5, , ]), type = "l", lty = 1, col = "https://github.com/mrc-ide/odin/issues/5", xlab = "Day", ylab = "Daily incidence") points(cases ~ day, incidence, col = "blue", pch = 19) ` code.txt