In kimfinale/particlefilter:
# theta <- data.table::fread("inst/extdata/theta.csv")
# y0 <- data.table::fread("inst/extdata/y0.csv")
# usethis::use_data(theta, overwrite = T)
# usethis::use_data(y0, overwrite = T)
# devtools::document()
# y1 = y0$val
# names(y1) = y0$name
# y0 = y1
Check functions
process_model
nparticle = 100
tend = 20
# beta <- theta[name =="R0", val] * theta[name =="gamma", val]
beta = theta[["R0"]] * theta[["gamma"]]
# volatility for beta (random walk)
beta_vol <- matrix(rnorm(nparticle * tend, mean = 0, sd = theta[["betavol"]]), nrow = tend)
beta_vol[1,] <- exp(beta_vol[1,]) * beta # initial beta
y <- data.frame(matrix(0, nrow = nparticle, ncol = length(y0)))
names(y) <- names(y0)
y[,"S"] <- 1e7
y[,"I"] <- 100
yt <- y1 <- process_model(params = theta, y = y, beta = beta_vol[1,])
for (i in 1:100) {
yt <- process_model(params = theta, y = yt, beta = beta_vol[1,])
}
head(y1)
head(yt)
head(beta_vol[1,])
Rt data generation
nt <- 200 # days
presymp_infect <- TRUE
param <- list()
param$sigma <- 1/4
param$gamma <- 1/4
param$R <- c(rep(1.2, nt/4), rep(1.4, nt/4), rep(0.8, nt/4), rep(1.2, nt/4))
param$presymp_infect <- FALSE
times <- 0:(nt-1) # simulation times
y0 <- c(S = 1e7, E1 = 40, E2 = 40, I = 100, R = 0, CE1 = 0, CE2 = 0, CI = 0) # initial values
library(deSolve)
library(tidyverse)
ode(func = se1e2ir, y = y0, times = times, parms = param) %>%
as.data.frame() -> out
inf_daily <- diff(out$CE1)
onset_daily <- diff(out$CI)
confirm_daily <- diff(out$R)
df <- data.frame(t = 1:length(param$R),
daily_R_true = param$R,
daily_infect = c(0, inf_daily),
daily_onset = c(0, onset_daily),
daily_confirm = c(0, confirm_daily))
df %>%
pivot_longer(cols = -t, names_to = "var") %>%
filter(var != "daily_R_true") %>%
ggplot(aes(t, value, color = var)) +
geom_line() +
labs(x = "time (day)", y = "number of individuals")
# Rt_data <- df
# usethis::use_data(Rt_data, overwrite = T)
particle_filter - sample
pf <- particle_filter(data_type = "infection", npart = 1e2)
pf$lik_overall_average
df <- cbind(Rt_data[,1:2], Rt_pf = pf$trace$beta/theta[["gamma"]])
suppressPackageStartupMessages(library(tidyverse))
df %>% pivot_longer(cols = -t) %>%
ggplot(aes(t, value, color = name)) +
geom_line()
particle_filter - filtered distribution
theta[["betavol"]] = 0.2
pf <- particle_filter(data_type = "symptom onset", npart = 1e3)
pf$lik_overall_average
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)})))
Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]]))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
suppressPackageStartupMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day")
Particle filtering and smoothing
$$X(t_k) | Y(t_1)=y^_1, \dots,\; Y(t_n)=y^_n,\; X_k | Y_1=y_1^, \dots,\; Y_n=y_n^$$
res <- run_model(data_type = "confirmation", rep = 1e2, npart = 1e3, systematic_resampling = T)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)})))
library(data.table)
Rt_mean <- data.frame(mean = rowMeans(res$Rt),
rollmean = frollmean(Rt_quantile[,3], n=7))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
library(tidyverse)
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering with smooothing", y = "R(t)", x = "Day")
repeat particle_filter and sample (backward sampling)
res <- run_model(data_type = "confirmation", rep = 1e3, npart = 1e3,
systematic_resampling = F)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)})))
library(data.table)
Rt_mean <- data.frame(mean = rowMeans(res$Rt),
rollmean = frollmean(Rt_quantile[,3], n=7))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
suppressPackageStartupMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day")
# library(data.table)
# df <- data.frame(pf_Rt_rollmean = frollmean(Rt_quantile[,3], n=7))
# df <- cbind(Rt_data[,1:2], df)
# df %>%
# pivot_longer(cols = - t) %>%
# ggplot(aes(t, value, color = name)) +
# geom_line() +
# geom_point()
Maximum likelihood estimation test
# devtools::load_all(".")
mle <- maxlik(data_type = "confirmation")
df <- cbind(Rt_data[,1:2], (mle$beta_trace/theta[name == "gamma", val]))
suppressPackageStartupMessages(library(tidyverse))
df %>%
pivot_longer(cols = - t) %>%
ggplot(aes(t, value, color = name)) +
geom_line() +
geom_point() +
labs(title = "R(t) via maximum likelihood", y = "R(t)", x = "Day")
RcppArmadillo 3D array test
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/test.cpp")
a = array(1, dim=c(12,12,12))
s = R_C_sum(a)
# > (s = R_C_sum(a))
# [1] 1728
# > 12**3
# [1] 1728
a2 = C_R_times_two(a)
# > class(a2)
# [1] "array"
# > dim(a2)
# [1] 12 12 12
# > a2[1:2,1:2,1:2]
# , , 1
#
# [,1] [,2]
# [1,] 2 2
# [2,] 2 2
#
# , , 2
#
# [,1] [,2]
# [1,] 2 2
# [2,] 2 2
arr = array(0, dim=c(2,2,2))
cube_with_ones(arr)
Rcpp integer casting test
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/test.cpp")
count_dt(tbegin = 10, tend = 12, dt = 0.1)
Rcpp process model test
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/process_model_cpp.cpp")
y0mat = matrix(y0, nrow = 1)
res = process_model_cpp(params = theta, y = y0mat, tbegin = 0, tend = 2, dt = 0.1, beta = 0.3)
beta = 0.5
tend = 200
res = matrix(NA, nrow = tend + 1, ncol = length(y0))
m = matrix(y0, nrow = 1)
res[1, ] = m[1, ]
for (i in 1:tend){
m = process_model_cpp(params = theta, y = m, tbegin = 0, tend = 1, dt = 0.2, beta = beta)
res[i + 1, ] = m[1, ]
}
res = as.data.frame(res)
names(res) = c("S", "E1", "E2", "I", "R", "CE1", "CE2", "CI")
suppressMessages(library(tidyverse))
res %>% mutate(day = (1:n() - 1)) %>%
pivot_longer(cols = - day) %>%
dplyr::filter(name == "I" | name == "R") %>%
ggplot(aes(x = day)) +
geom_line(aes(y = value, color = name))
## outbreak size comparision - to have a sense that numerical solutions are reasonable
R0 = beta / theta[["gamma"]]
outbreak_size <- function (R0, final) {
return (1 - exp( - R0 * final) - final)
}
sols <- rootSolve::multiroot( outbreak_size, c(0, 1), R0=2 )
max(sols$root) # 0 is also a solution
y0[["S"]] * max(sols$root)
res[tend+1, "R"]
Testing particle filter implemented in C++
particle_filter_cpp - filtered distribution
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/particle_filter_cpp.cpp")
devtools::load_all(".")
theta[["betavol"]] = 0.2
dmat = data.matrix(Rt_data)
pf = particle_filter_cpp(theta, y0, dmat, "infection", 1e3L, 200L, 0.1)
pf$lik_overall_average
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)})))
Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]]))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
suppressPackageStartupMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day")
repeat particle_filter and sample (backward sampling)
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/particle_filter_cpp2.cpp")
devtools::load_all(".")
theta[["betavol"]] = 0.2
res <- run_model_cpp(data_type = "infection", rep = 1e2, npart = 1e3)
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)})))
library(data.table)
Rt_mean <- data.frame(mean = rowMeans(res$Rt),
rollmean = frollmean(Rt_quantile[,3], n=7))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
suppressPackageStartupMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day")
performance comparison
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/particle_filter_cpp.cpp")
devtools::load_all(".")
library(microbenchmark)
microbenchmark(run_model_cpp(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1),
run_model(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1), times = 50L)
Parallel
library(parallel)
library(doParallel)
library(foreach)
ncores = detectCores()
cl <- makeCluster(getOption("cl.cores", ncores-2))
doParallel::registerDoParallel(cl)
system.time({
r <- foreach (i=1:1000, .combine = cbind, .packages ="particlefilter", .inorder = F) %dopar% {
run_model(params = theta, y = y0, data = Rt_data, data_type = "infection",
rep = 1, npart = 1e3, tend = 200, dt = 0.1,
systematic_resampling = FALSE)
}})
parallel::stopCluster(cl)
rt = r[9,] # Rt
rtdf = data.frame(matrix(unlist(rt), ncol=length(rt), byrow=F))
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(rtdf, 1, function(x) {quantile(x, pr)})))
library(data.table)
Rt_mean <- data.frame(mean = rowMeans(rtdf),
rollmean = frollmean(Rt_quantile[,3], n=7))
df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2])
suppressPackageStartupMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
geom_point(aes(y = daily_R_true), color = "darkred") +
labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day")
# library(tidyverse)
# rtdf %>%
# mutate(day = 1:n() - 1) %>%
# pivot_longer(cols = - day) %>%
# ggplot(aes(x = day))+
# geom_line(aes(y = value))
# cl <- makeCluster(getOption("cl.cores", 10))
# doParallel::registerDoParallel(cl)
# clusterExport(cl = cl, varlist = c("run_model_cpp", "y0", "theta", "Rt_data", "particle_filter_cpp"))
# # system.time({
# # r <- foreach (i=1:50, .combine=cbind, .packages = "particlefilter") %dopar% {
# # run_model(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1)
# # }})
# system.time({
# r <- foreach (i=1:50, .combine=cbind) %dopar% {
# run_model_cpp(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1)
# }})
# parallel::stopCluster(cl)
KDCA data
Data cleaning
library(readxl)
suppressMessages(library(tidyverse))
d <- read_xlsx("data/covid_kdca.xlsx")
names(d) = c("id", "nationality", "sex", "age", "sido", "sigungu", "date_symptom_onset", "date_diagnosis", "date_report", "route_infection", "country_orgin", "id_infector", "occupation")
d1 <- d %>% filter(route_infection == "국내")%>%
select(date_diagnosis) %>%
group_by(date = date_diagnosis) %>%
summarise(daily_confirm = n())
d2 <- data.frame(date = seq(as.Date("2020-01-03"), as.Date("2021-02-23"), by = "day"))
d2$t <- 1:nrow(d2)
d3 <- left_join(d2, d1, by = "date")
dat <- d3 %>% mutate(daily_confirm = ifelse(is.na(daily_confirm), 0, daily_confirm))
## name standardization
ggplot(dat, aes(t, daily_confirm)) +
geom_line()
Data fitting
theta[["betavol"]] = 0.6
ynew <- c(S = 5*1e7, E1 = 40, E2 = 40, I = 1e2, R = 0, CE1 = 0, CE2 = 0, CI = 0)
pf <- particle_filter(ydata = dat, data_type = "confirmation", npart = 1e4, tend = 418)
pf$lik_overall_average
pr <- c(0.025, 0.25, 0.5, 0.75, 0.975)
Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)})))
Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]]))
df <- cbind(Rt_quantile, Rt_mean, dat[,c("t", "date")])
suppressMessages(library(tidyverse))
ggplot(df, aes(x = t)) +
geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") +
geom_line(aes(y = `50%`), size = 1) +
geom_line(aes(y = mean), color = "steelblue", size = 1) +
labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day") +
geom_hline(yintercept = 1, color = "darkred", linetype = "dotted") +
scale_y_continuous(limits = c(0, 100))
ci <- t(pf$latent_var_filtered[,,8])
daily_confirm <- diff(ci)
for (i in 1:ncol(daily_confirm)){
if (sum(daily_confirm[, i] < 0) > 0) {
cat("part =", i, "\n"); stop("negative")
}
}
plot(daily_confirm[,1], cex = 0.4)
abline(h = 0, col = 2)
points(dat$t, dat$daily_confirm, col = 3, cex = 0.5)
sum(daily_confirm < 0)
ci <- t(apply(pf$latent_var_filtered[,,8], 2, function(x) {quantile(x, pr)}))
daily_confirm <- diff(ci)
Data fitting
library(parallel)
library(doParallel)
library(foreach)
ncores = detectCores()
cl <- makeCluster(getOption("cl.cores", ncores-2))
doParallel::registerDoParallel(cl)
clusterExport(cl = cl, varlist = c("dat"), envir = .GlobalEnv)
system.time({
r <- foreach (i=1:1000, .combine = cbind, .packages ="particlefilter", .inorder = F) %dopar% {
run_model(params = theta, y = y0, data = dat, data_type = "confirmation",
rep = 1, npart = 1e3, tend = 200, dt = 0.1,
systematic_resampling = FALSE)
}})
parallel::stopCluster(cl)
Testing arma::cube can be access through _ just like in Rcpp::NumericMatrix
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/test.cpp")
# a = array(1:12, c(2,3,2))
# a[, , 1]
# (a = export_array(Q = a, slice = 0))
# # (m = export_mat(Q = a, slice = 0))
# (v = export_vec(Q = a, col = 0, slice = 0))
# (i = to_int(x = c(1.0, 2.2)))
# m = matrix(1:12, nrow = 3)
# m[,1]
# (v = export_vec_from_mat(m = m, col = 0))
# (ch = add_x("mr. "))
# all_x("x")
# all_x("y")
int_sample(10)
int_wt_sample(10, wt = c(0.8,0.2,rep(0,8)))
sample_test(10, wt = c(0.8,0.2,rep(0,8)))
testing assign_weights function
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/assign_weights_cpp.cpp")
st_now = data.matrix(data.frame(S=rep(1,10), E1 = rep(1,10), E2 = rep(1,10), I = rep(1,10), R = 1:10, CE1 = 101:110, CI = 101:110))
st_before = data.matrix(data.frame(S=rep(1,10), E1 = rep(1,10), E2=rep(1,10), I = rep(1,10), R = 1:10-1, CE1 = round(seq(1, 100, length.out = 10)), CI = round(seq(1, 100, length.out = 10))))
wt = assign_weights_cpp(st_now, st_before, 53L, data.matrix(Rt_data), "infection")
all.equal(wt[1], dpois(101,100))
testing particle_filter_cpp function
library(Rcpp)
library(RcppArmadillo)
sourceCpp("cpp/particle_filter_cpp.cpp")
Conditional log likelihood
The estimated conditional log likelihood from a fitted model.
The conditional likelihood is defined to be the value of the density
$$ Y(t_k) | Y(t_1),\dots,Y(t_{k-1}) Y_k | Y_1,\dots,Y_{k-1}$$
evaluated at $$Y(t_k) = y^_k Y_k = y_k^$$.
Here, \eqn{Y(t_k)}{Yk} is the observable process, and \eqn{y^_k}{yk} the data, at time \eqn{t_k}.
Thus the conditional log likelihood at time \eqn{t_k} is
\deqn{\ell_k(\theta) = \log f[Y(t_k)=y^_k \vert Y(t_1)=y^1, \dots, Y(t{k-1})=y^_{k-1}],}{ell_k(theta)=log f[Yk = yk | Y1=y1, \dots, Y(k-1)=y(k-1)],}
where \eqn{f} is the probability density above.
pomp package
library(tidyverse)
read_csv(paste0("https://kingaa.github.io/sbied/stochsim/", "Measles_Consett_1948.csv")) %>%
select(week,reports=cases) -> meas
meas %>% as.data.frame() %>% head()
sir_step <- function (S, I, R, N, Beta, mu_IR, delta.t, ...) {
dN_SI <- rbinom(n=1,size=S,prob=1-exp(-Beta*I/N*delta.t))
dN_IR <- rbinom(n=1,size=I,prob=1-exp(-mu_IR*delta.t))
S <- S - dN_SI
I <- I + dN_SI - dN_IR
R <- R + dN_IR
c(S = S, I = I, R = R)
}
sir_rinit <- function (N, eta, ...) {
c(S = round(N*eta), I = 1, R = round(N*(1-eta)))
}
library(pomp)
meas %>%
pomp(times="week",t0=0,
rprocess=euler(sir_step,delta.t=1/7),
rinit=sir_rinit
) -> measSIR
sir_dmeas <- function (reports, H, rho, log, ...) {
dbinom(x=reports, size=H, prob=rho, log=log)
}
sir_rmeas <- function (H, rho, ...) {
c(reports=rbinom(n=1, size=H, prob=rho))
}
measSIR %>%
pomp(
rmeasure=sir_rmeas,
dmeasure=sir_dmeas
) -> measSIR
measSIR %>%
simulate(
params=c(Beta=7.5,mu_IR=0.5,rho=0.5,eta=0.03,N=38000),
nsim=20,format="data.frame",include.data=TRUE
) -> sims
sims %>%
ggplot(aes(x=week,y=reports,group=.id,color=.id=="data"))+
geom_line()+
guides(color=FALSE)
library(pomp)
library(tidyverse)
seir_step <- Csnippet("
double N = S + E1 + E2 + I + R;
double dN_SE1 = beta * S * I / N * dt;
double dN_E1E2 = 2 * sigma * E1 * dt;
double dN_E2I = 2 * sigma * E2 * dt;
double dN_IR = gamma * I * dt;
S += - dN_SE1;
E1 += dN_SE1 - dN_E1E2;
E2 += dN_E1E2 - dN_E2I;
I += dN_E2I - dN_IR;
R += dN_IR;
H += dN_SE1;
")
seir_rinit <- Csnippet("
S = S_0;
E1 = E1_0;
E2 = E2_0;
I = I_0;
R = R_0;
H = H_0;
")
seir_dmeas <- Csnippet("
lik = dpois(daily_infect, H, give_log);
")
seir_rmeas <- Csnippet("
daily_infect = rpois(H);
")
dat <- data.frame(day = Rt_data$t, daily_infect = Rt_data$daily_infect)
dat %>% pomp(
times = "day",
t0 = 1,
rprocess = euler(seir_step, delta.t = 0.2),
rinit = seir_rinit,
rmeasure = seir_rmeas,
dmeasure = seir_dmeas,
accumvars = "H",
statenames = c("S", "E1", "E2", "I", "R", "H"),
paramnames = c("beta", "sigma", "gamma", "S_0", "E1_0", "E2_0", "I_0", "R_0", "H_0")) -> covid_seir
covid_seir %>%
simulate(params = c(beta=1.2*0.25, sigma=0.25, gamma=0.25,
S_0=1e7, E1_0=40, E2_0=40, I_0=100, R_0=0, H_0=0),
nsim = 20, format="data.frame", include.data=TRUE) -> sims
sims %>%
ggplot(aes(x=day,y=daily_infect, group=.id, color=.id=="data"))+
geom_line()+
guides(color=FALSE)
covid_seir %>%
pfilter(Np=1000, paramnames="beta") -> pfrick
library(pomp)
library(tidyverse)
seir_step <- Csnippet("
double N = S + E1 + E2 + I + R;
double dN_SE1 = beta * S * I / N * dt;
double dN_E1E2 = 2 * 0.25 * E1 * dt;
double dN_E2I = 2 * 0.25 * E2 * dt;
double dN_IR = 0.25 * I * dt;
S += - dN_SE1;
E1 += dN_SE1 - dN_E1E2;
E2 += dN_E1E2 - dN_E2I;
I += dN_E2I - dN_IR;
R += dN_IR;
H += dN_SE1;
")
seir_rinit <- Csnippet("
S = 1e7;
E1 = 40;
E2 = 40;
I = 100;
R = 0;
H = 0;
")
seir_dmeas <- Csnippet("
lik = dpois(daily_infect, H, give_log);
")
seir_rmeas <- Csnippet("
daily_infect = rpois(H);
")
dat <- data.frame(day = Rt_data$t, daily_infect = as.integer(Rt_data$daily_infect))
dat %>% pomp(
times = "day",
t0 = 1,
rprocess = euler(seir_step, delta.t = 0.2),
rinit = seir_rinit,
rmeasure = seir_rmeas,
dmeasure = seir_dmeas,
accumvars = "H",
statenames = c("S", "E1", "E2", "I", "R", "H"),
paramnames = c("beta")) -> covid_seir
covid_seir %>%
simulate(params = c(beta=1.2*0.25),
nsim = 20, format="data.frame", include.data=TRUE) -> sims
sims %>%
ggplot(aes(x=day,y=daily_infect, group=.id, color=.id=="data"))+
geom_line()+
guides(color=FALSE)
pf = pfilter(covid_seir, Np=1000, params = c(beta = 0.3),
filter.mean = T, filter.traj = T, verbose = T)
pf@filter.mean %>%
t() %>%
as.data.frame() %>%
mutate(time = 1:n()) %>%
pivot_longer(-time) %>%
filter(name == "I") %>%
ggplot(aes(time, value)) +
geom_line(aes(color = name))
ewmeas %>%
subset(time < 1952) %>%
pomp(
times="time",t0=1948,
rprocess=euler(
Csnippet("
int nrate = 6;
double rate[nrate]; // transition rates
double trans[nrate]; // transition numbers
double dW;
// gamma noise, mean=dt, variance=(sigma^2 dt)
dW = rgammawn(sigma,dt);
// compute the transition rates
rate[0] = mu*pop; // birth into susceptible class
rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection
rate[2] = mu; // death from susceptible class
rate[3] = gamma; // recovery
rate[4] = mu; // death from infectious class
rate[5] = mu; // death from recovered class
// compute the transition numbers
trans[0] = rpois(rate[0]*dt); // births are Poisson
reulermultinom(2,S,&rate[1],dt,&trans[1]);
reulermultinom(2,I,&rate[3],dt,&trans[3]);
reulermultinom(1,R,&rate[5],dt,&trans[5]);
// balance the equations
S += trans[0]-trans[1]-trans[2];
I += trans[1]-trans[3]-trans[4];
R += trans[3]-trans[5];
"),
delta.t=1/52/20
),
rinit=Csnippet("
double m = pop/(S_0+I_0+R_0);
S = nearbyint(m*S_0);
I = nearbyint(m*I_0);
R = nearbyint(m*R_0);
"),
paramnames=c("mu","pop","iota","gamma","Beta","sigma",
"S_0","I_0","R_0"),
statenames=c("S","I","R"),
params=c(mu=1/50,iota=10,pop=50e6,gamma=26,Beta=400,sigma=0.1,
S_0=0.07,I_0=0.001,R_0=0.93)
) -> ew1
kimfinale/particlefilter documentation built on March 6, 2021, 9:49 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
# theta <- data.table::fread("inst/extdata/theta.csv") # y0 <- data.table::fread("inst/extdata/y0.csv") # usethis::use_data(theta, overwrite = T) # usethis::use_data(y0, overwrite = T) # devtools::document() # y1 = y0$val # names(y1) = y0$name # y0 = y1
Check functions
process_model
nparticle = 100 tend = 20 # beta <- theta[name =="R0", val] * theta[name =="gamma", val] beta = theta[["R0"]] * theta[["gamma"]] # volatility for beta (random walk) beta_vol <- matrix(rnorm(nparticle * tend, mean = 0, sd = theta[["betavol"]]), nrow = tend) beta_vol[1,] <- exp(beta_vol[1,]) * beta # initial beta y <- data.frame(matrix(0, nrow = nparticle, ncol = length(y0))) names(y) <- names(y0) y[,"S"] <- 1e7 y[,"I"] <- 100 yt <- y1 <- process_model(params = theta, y = y, beta = beta_vol[1,]) for (i in 1:100) { yt <- process_model(params = theta, y = yt, beta = beta_vol[1,]) } head(y1) head(yt) head(beta_vol[1,])
Rt data generation
nt <- 200 # days presymp_infect <- TRUE param <- list() param$sigma <- 1/4 param$gamma <- 1/4 param$R <- c(rep(1.2, nt/4), rep(1.4, nt/4), rep(0.8, nt/4), rep(1.2, nt/4)) param$presymp_infect <- FALSE times <- 0:(nt-1) # simulation times y0 <- c(S = 1e7, E1 = 40, E2 = 40, I = 100, R = 0, CE1 = 0, CE2 = 0, CI = 0) # initial values library(deSolve) library(tidyverse) ode(func = se1e2ir, y = y0, times = times, parms = param) %>% as.data.frame() -> out inf_daily <- diff(out$CE1) onset_daily <- diff(out$CI) confirm_daily <- diff(out$R) df <- data.frame(t = 1:length(param$R), daily_R_true = param$R, daily_infect = c(0, inf_daily), daily_onset = c(0, onset_daily), daily_confirm = c(0, confirm_daily)) df %>% pivot_longer(cols = -t, names_to = "var") %>% filter(var != "daily_R_true") %>% ggplot(aes(t, value, color = var)) + geom_line() + labs(x = "time (day)", y = "number of individuals") # Rt_data <- df # usethis::use_data(Rt_data, overwrite = T)
particle_filter - sample
pf <- particle_filter(data_type = "infection", npart = 1e2) pf$lik_overall_average df <- cbind(Rt_data[,1:2], Rt_pf = pf$trace$beta/theta[["gamma"]]) suppressPackageStartupMessages(library(tidyverse)) df %>% pivot_longer(cols = -t) %>% ggplot(aes(t, value, color = name)) + geom_line()
particle_filter - filtered distribution
theta[["betavol"]] = 0.2 pf <- particle_filter(data_type = "symptom onset", npart = 1e3) pf$lik_overall_average pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)}))) Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]])) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) suppressPackageStartupMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day")
Particle filtering and smoothing
$$X(t_k) | Y(t_1)=y^_1, \dots,\; Y(t_n)=y^_n,\; X_k | Y_1=y_1^, \dots,\; Y_n=y_n^$$
res <- run_model(data_type = "confirmation", rep = 1e2, npart = 1e3, systematic_resampling = T) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)}))) library(data.table) Rt_mean <- data.frame(mean = rowMeans(res$Rt), rollmean = frollmean(Rt_quantile[,3], n=7)) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) library(tidyverse) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering with smooothing", y = "R(t)", x = "Day")
repeat particle_filter and sample (backward sampling)
res <- run_model(data_type = "confirmation", rep = 1e3, npart = 1e3, systematic_resampling = F) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)}))) library(data.table) Rt_mean <- data.frame(mean = rowMeans(res$Rt), rollmean = frollmean(Rt_quantile[,3], n=7)) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) suppressPackageStartupMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day") # library(data.table) # df <- data.frame(pf_Rt_rollmean = frollmean(Rt_quantile[,3], n=7)) # df <- cbind(Rt_data[,1:2], df) # df %>% # pivot_longer(cols = - t) %>% # ggplot(aes(t, value, color = name)) + # geom_line() + # geom_point()
Maximum likelihood estimation test
# devtools::load_all(".") mle <- maxlik(data_type = "confirmation") df <- cbind(Rt_data[,1:2], (mle$beta_trace/theta[name == "gamma", val])) suppressPackageStartupMessages(library(tidyverse)) df %>% pivot_longer(cols = - t) %>% ggplot(aes(t, value, color = name)) + geom_line() + geom_point() + labs(title = "R(t) via maximum likelihood", y = "R(t)", x = "Day")
RcppArmadillo 3D array test
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/test.cpp") a = array(1, dim=c(12,12,12)) s = R_C_sum(a) # > (s = R_C_sum(a)) # [1] 1728 # > 12**3 # [1] 1728 a2 = C_R_times_two(a) # > class(a2) # [1] "array" # > dim(a2) # [1] 12 12 12 # > a2[1:2,1:2,1:2] # , , 1 # # [,1] [,2] # [1,] 2 2 # [2,] 2 2 # # , , 2 # # [,1] [,2] # [1,] 2 2 # [2,] 2 2 arr = array(0, dim=c(2,2,2)) cube_with_ones(arr)
Rcpp integer casting test
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/test.cpp") count_dt(tbegin = 10, tend = 12, dt = 0.1)
Rcpp process model test
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/process_model_cpp.cpp") y0mat = matrix(y0, nrow = 1) res = process_model_cpp(params = theta, y = y0mat, tbegin = 0, tend = 2, dt = 0.1, beta = 0.3) beta = 0.5 tend = 200 res = matrix(NA, nrow = tend + 1, ncol = length(y0)) m = matrix(y0, nrow = 1) res[1, ] = m[1, ] for (i in 1:tend){ m = process_model_cpp(params = theta, y = m, tbegin = 0, tend = 1, dt = 0.2, beta = beta) res[i + 1, ] = m[1, ] } res = as.data.frame(res) names(res) = c("S", "E1", "E2", "I", "R", "CE1", "CE2", "CI") suppressMessages(library(tidyverse)) res %>% mutate(day = (1:n() - 1)) %>% pivot_longer(cols = - day) %>% dplyr::filter(name == "I" | name == "R") %>% ggplot(aes(x = day)) + geom_line(aes(y = value, color = name)) ## outbreak size comparision - to have a sense that numerical solutions are reasonable R0 = beta / theta[["gamma"]] outbreak_size <- function (R0, final) { return (1 - exp( - R0 * final) - final) } sols <- rootSolve::multiroot( outbreak_size, c(0, 1), R0=2 ) max(sols$root) # 0 is also a solution y0[["S"]] * max(sols$root) res[tend+1, "R"]
Testing particle filter implemented in C++
particle_filter_cpp - filtered distribution
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/particle_filter_cpp.cpp") devtools::load_all(".") theta[["betavol"]] = 0.2 dmat = data.matrix(Rt_data) pf = particle_filter_cpp(theta, y0, dmat, "infection", 1e3L, 200L, 0.1) pf$lik_overall_average pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)}))) Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]])) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) suppressPackageStartupMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day")
repeat particle_filter and sample (backward sampling)
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/particle_filter_cpp2.cpp") devtools::load_all(".") theta[["betavol"]] = 0.2 res <- run_model_cpp(data_type = "infection", rep = 1e2, npart = 1e3) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(res$Rt, 1, function(x) {quantile(x, pr)}))) library(data.table) Rt_mean <- data.frame(mean = rowMeans(res$Rt), rollmean = frollmean(Rt_quantile[,3], n=7)) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) suppressPackageStartupMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day")
performance comparison
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/particle_filter_cpp.cpp") devtools::load_all(".") library(microbenchmark) microbenchmark(run_model_cpp(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1), run_model(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1), times = 50L)
Parallel
library(parallel) library(doParallel) library(foreach) ncores = detectCores() cl <- makeCluster(getOption("cl.cores", ncores-2)) doParallel::registerDoParallel(cl) system.time({ r <- foreach (i=1:1000, .combine = cbind, .packages ="particlefilter", .inorder = F) %dopar% { run_model(params = theta, y = y0, data = Rt_data, data_type = "infection", rep = 1, npart = 1e3, tend = 200, dt = 0.1, systematic_resampling = FALSE) }}) parallel::stopCluster(cl) rt = r[9,] # Rt rtdf = data.frame(matrix(unlist(rt), ncol=length(rt), byrow=F)) pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(rtdf, 1, function(x) {quantile(x, pr)}))) library(data.table) Rt_mean <- data.frame(mean = rowMeans(rtdf), rollmean = frollmean(Rt_quantile[,3], n=7)) df <- cbind(Rt_quantile, Rt_mean, Rt_data[,1:2]) suppressPackageStartupMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + geom_point(aes(y = daily_R_true), color = "darkred") + labs(title = "R(t) via particle filtering with backward sampling", y = "R(t)", x = "Day") # library(tidyverse) # rtdf %>% # mutate(day = 1:n() - 1) %>% # pivot_longer(cols = - day) %>% # ggplot(aes(x = day))+ # geom_line(aes(y = value)) # cl <- makeCluster(getOption("cl.cores", 10)) # doParallel::registerDoParallel(cl) # clusterExport(cl = cl, varlist = c("run_model_cpp", "y0", "theta", "Rt_data", "particle_filter_cpp")) # # system.time({ # # r <- foreach (i=1:50, .combine=cbind, .packages = "particlefilter") %dopar% { # # run_model(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1) # # }}) # system.time({ # r <- foreach (i=1:50, .combine=cbind) %dopar% { # run_model_cpp(data_type = "infection", rep = 1, npart = 1e3, dt = 0.1) # }}) # parallel::stopCluster(cl)
KDCA data
Data cleaning
library(readxl) suppressMessages(library(tidyverse)) d <- read_xlsx("data/covid_kdca.xlsx") names(d) = c("id", "nationality", "sex", "age", "sido", "sigungu", "date_symptom_onset", "date_diagnosis", "date_report", "route_infection", "country_orgin", "id_infector", "occupation") d1 <- d %>% filter(route_infection == "국내")%>% select(date_diagnosis) %>% group_by(date = date_diagnosis) %>% summarise(daily_confirm = n()) d2 <- data.frame(date = seq(as.Date("2020-01-03"), as.Date("2021-02-23"), by = "day")) d2$t <- 1:nrow(d2) d3 <- left_join(d2, d1, by = "date") dat <- d3 %>% mutate(daily_confirm = ifelse(is.na(daily_confirm), 0, daily_confirm)) ## name standardization ggplot(dat, aes(t, daily_confirm)) + geom_line()
Data fitting
theta[["betavol"]] = 0.6 ynew <- c(S = 5*1e7, E1 = 40, E2 = 40, I = 1e2, R = 0, CE1 = 0, CE2 = 0, CI = 0) pf <- particle_filter(ydata = dat, data_type = "confirmation", npart = 1e4, tend = 418) pf$lik_overall_average pr <- c(0.025, 0.25, 0.5, 0.75, 0.975) Rt_quantile <- as.data.frame(t(apply(pf$beta_filtered/theta[["gamma"]], 1, function(x) {quantile(x, pr)}))) Rt_mean <- data.frame(mean = rowMeans(pf$beta_filtered/theta[["gamma"]])) df <- cbind(Rt_quantile, Rt_mean, dat[,c("t", "date")]) suppressMessages(library(tidyverse)) ggplot(df, aes(x = t)) + geom_ribbon(aes(ymin = `2.5%`, ymax = `97.5%`), fill = "grey") + geom_line(aes(y = `50%`), size = 1) + geom_line(aes(y = mean), color = "steelblue", size = 1) + labs(title = "R(t) via particle filtering", y = "R(t)", x = "Day") + geom_hline(yintercept = 1, color = "darkred", linetype = "dotted") + scale_y_continuous(limits = c(0, 100)) ci <- t(pf$latent_var_filtered[,,8]) daily_confirm <- diff(ci) for (i in 1:ncol(daily_confirm)){ if (sum(daily_confirm[, i] < 0) > 0) { cat("part =", i, "\n"); stop("negative") } } plot(daily_confirm[,1], cex = 0.4) abline(h = 0, col = 2) points(dat$t, dat$daily_confirm, col = 3, cex = 0.5) sum(daily_confirm < 0) ci <- t(apply(pf$latent_var_filtered[,,8], 2, function(x) {quantile(x, pr)})) daily_confirm <- diff(ci)
Data fitting
library(parallel) library(doParallel) library(foreach) ncores = detectCores() cl <- makeCluster(getOption("cl.cores", ncores-2)) doParallel::registerDoParallel(cl) clusterExport(cl = cl, varlist = c("dat"), envir = .GlobalEnv) system.time({ r <- foreach (i=1:1000, .combine = cbind, .packages ="particlefilter", .inorder = F) %dopar% { run_model(params = theta, y = y0, data = dat, data_type = "confirmation", rep = 1, npart = 1e3, tend = 200, dt = 0.1, systematic_resampling = FALSE) }}) parallel::stopCluster(cl)
Testing arma::cube can be access through _ just like in Rcpp::NumericMatrix
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/test.cpp") # a = array(1:12, c(2,3,2)) # a[, , 1] # (a = export_array(Q = a, slice = 0)) # # (m = export_mat(Q = a, slice = 0)) # (v = export_vec(Q = a, col = 0, slice = 0)) # (i = to_int(x = c(1.0, 2.2))) # m = matrix(1:12, nrow = 3) # m[,1] # (v = export_vec_from_mat(m = m, col = 0)) # (ch = add_x("mr. ")) # all_x("x") # all_x("y") int_sample(10) int_wt_sample(10, wt = c(0.8,0.2,rep(0,8))) sample_test(10, wt = c(0.8,0.2,rep(0,8)))
testing assign_weights function
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/assign_weights_cpp.cpp") st_now = data.matrix(data.frame(S=rep(1,10), E1 = rep(1,10), E2 = rep(1,10), I = rep(1,10), R = 1:10, CE1 = 101:110, CI = 101:110)) st_before = data.matrix(data.frame(S=rep(1,10), E1 = rep(1,10), E2=rep(1,10), I = rep(1,10), R = 1:10-1, CE1 = round(seq(1, 100, length.out = 10)), CI = round(seq(1, 100, length.out = 10)))) wt = assign_weights_cpp(st_now, st_before, 53L, data.matrix(Rt_data), "infection") all.equal(wt[1], dpois(101,100))
testing particle_filter_cpp function
library(Rcpp) library(RcppArmadillo) sourceCpp("cpp/particle_filter_cpp.cpp")
Conditional log likelihood
The estimated conditional log likelihood from a fitted model. The conditional likelihood is defined to be the value of the density $$ Y(t_k) | Y(t_1),\dots,Y(t_{k-1}) Y_k | Y_1,\dots,Y_{k-1}$$ evaluated at $$Y(t_k) = y^_k Y_k = y_k^$$.
Here, \eqn{Y(t_k)}{Yk} is the observable process, and \eqn{y^_k}{yk} the data, at time \eqn{t_k}.
Thus the conditional log likelihood at time \eqn{t_k} is \deqn{\ell_k(\theta) = \log f[Y(t_k)=y^_k \vert Y(t_1)=y^1, \dots, Y(t{k-1})=y^_{k-1}],}{ell_k(theta)=log f[Yk = yk | Y1=y1, \dots, Y(k-1)=y(k-1)],} where \eqn{f} is the probability density above.
pomp package
library(tidyverse) read_csv(paste0("https://kingaa.github.io/sbied/stochsim/", "Measles_Consett_1948.csv")) %>% select(week,reports=cases) -> meas meas %>% as.data.frame() %>% head() sir_step <- function (S, I, R, N, Beta, mu_IR, delta.t, ...) { dN_SI <- rbinom(n=1,size=S,prob=1-exp(-Beta*I/N*delta.t)) dN_IR <- rbinom(n=1,size=I,prob=1-exp(-mu_IR*delta.t)) S <- S - dN_SI I <- I + dN_SI - dN_IR R <- R + dN_IR c(S = S, I = I, R = R) } sir_rinit <- function (N, eta, ...) { c(S = round(N*eta), I = 1, R = round(N*(1-eta))) } library(pomp) meas %>% pomp(times="week",t0=0, rprocess=euler(sir_step,delta.t=1/7), rinit=sir_rinit ) -> measSIR sir_dmeas <- function (reports, H, rho, log, ...) { dbinom(x=reports, size=H, prob=rho, log=log) } sir_rmeas <- function (H, rho, ...) { c(reports=rbinom(n=1, size=H, prob=rho)) } measSIR %>% pomp( rmeasure=sir_rmeas, dmeasure=sir_dmeas ) -> measSIR measSIR %>% simulate( params=c(Beta=7.5,mu_IR=0.5,rho=0.5,eta=0.03,N=38000), nsim=20,format="data.frame",include.data=TRUE ) -> sims sims %>% ggplot(aes(x=week,y=reports,group=.id,color=.id=="data"))+ geom_line()+ guides(color=FALSE)
library(pomp) library(tidyverse) seir_step <- Csnippet(" double N = S + E1 + E2 + I + R; double dN_SE1 = beta * S * I / N * dt; double dN_E1E2 = 2 * sigma * E1 * dt; double dN_E2I = 2 * sigma * E2 * dt; double dN_IR = gamma * I * dt; S += - dN_SE1; E1 += dN_SE1 - dN_E1E2; E2 += dN_E1E2 - dN_E2I; I += dN_E2I - dN_IR; R += dN_IR; H += dN_SE1; ") seir_rinit <- Csnippet(" S = S_0; E1 = E1_0; E2 = E2_0; I = I_0; R = R_0; H = H_0; ") seir_dmeas <- Csnippet(" lik = dpois(daily_infect, H, give_log); ") seir_rmeas <- Csnippet(" daily_infect = rpois(H); ") dat <- data.frame(day = Rt_data$t, daily_infect = Rt_data$daily_infect) dat %>% pomp( times = "day", t0 = 1, rprocess = euler(seir_step, delta.t = 0.2), rinit = seir_rinit, rmeasure = seir_rmeas, dmeasure = seir_dmeas, accumvars = "H", statenames = c("S", "E1", "E2", "I", "R", "H"), paramnames = c("beta", "sigma", "gamma", "S_0", "E1_0", "E2_0", "I_0", "R_0", "H_0")) -> covid_seir covid_seir %>% simulate(params = c(beta=1.2*0.25, sigma=0.25, gamma=0.25, S_0=1e7, E1_0=40, E2_0=40, I_0=100, R_0=0, H_0=0), nsim = 20, format="data.frame", include.data=TRUE) -> sims sims %>% ggplot(aes(x=day,y=daily_infect, group=.id, color=.id=="data"))+ geom_line()+ guides(color=FALSE) covid_seir %>% pfilter(Np=1000, paramnames="beta") -> pfrick
library(pomp) library(tidyverse) seir_step <- Csnippet(" double N = S + E1 + E2 + I + R; double dN_SE1 = beta * S * I / N * dt; double dN_E1E2 = 2 * 0.25 * E1 * dt; double dN_E2I = 2 * 0.25 * E2 * dt; double dN_IR = 0.25 * I * dt; S += - dN_SE1; E1 += dN_SE1 - dN_E1E2; E2 += dN_E1E2 - dN_E2I; I += dN_E2I - dN_IR; R += dN_IR; H += dN_SE1; ") seir_rinit <- Csnippet(" S = 1e7; E1 = 40; E2 = 40; I = 100; R = 0; H = 0; ") seir_dmeas <- Csnippet(" lik = dpois(daily_infect, H, give_log); ") seir_rmeas <- Csnippet(" daily_infect = rpois(H); ") dat <- data.frame(day = Rt_data$t, daily_infect = as.integer(Rt_data$daily_infect)) dat %>% pomp( times = "day", t0 = 1, rprocess = euler(seir_step, delta.t = 0.2), rinit = seir_rinit, rmeasure = seir_rmeas, dmeasure = seir_dmeas, accumvars = "H", statenames = c("S", "E1", "E2", "I", "R", "H"), paramnames = c("beta")) -> covid_seir covid_seir %>% simulate(params = c(beta=1.2*0.25), nsim = 20, format="data.frame", include.data=TRUE) -> sims sims %>% ggplot(aes(x=day,y=daily_infect, group=.id, color=.id=="data"))+ geom_line()+ guides(color=FALSE) pf = pfilter(covid_seir, Np=1000, params = c(beta = 0.3), filter.mean = T, filter.traj = T, verbose = T) pf@filter.mean %>% t() %>% as.data.frame() %>% mutate(time = 1:n()) %>% pivot_longer(-time) %>% filter(name == "I") %>% ggplot(aes(time, value)) + geom_line(aes(color = name))
ewmeas %>% subset(time < 1952) %>% pomp( times="time",t0=1948, rprocess=euler( Csnippet(" int nrate = 6; double rate[nrate]; // transition rates double trans[nrate]; // transition numbers double dW; // gamma noise, mean=dt, variance=(sigma^2 dt) dW = rgammawn(sigma,dt); // compute the transition rates rate[0] = mu*pop; // birth into susceptible class rate[1] = (iota+Beta*I*dW/dt)/pop; // force of infection rate[2] = mu; // death from susceptible class rate[3] = gamma; // recovery rate[4] = mu; // death from infectious class rate[5] = mu; // death from recovered class // compute the transition numbers trans[0] = rpois(rate[0]*dt); // births are Poisson reulermultinom(2,S,&rate[1],dt,&trans[1]); reulermultinom(2,I,&rate[3],dt,&trans[3]); reulermultinom(1,R,&rate[5],dt,&trans[5]); // balance the equations S += trans[0]-trans[1]-trans[2]; I += trans[1]-trans[3]-trans[4]; R += trans[3]-trans[5]; "), delta.t=1/52/20 ), rinit=Csnippet(" double m = pop/(S_0+I_0+R_0); S = nearbyint(m*S_0); I = nearbyint(m*I_0); R = nearbyint(m*R_0); "), paramnames=c("mu","pop","iota","gamma","Beta","sigma", "S_0","I_0","R_0"), statenames=c("S","I","R"), params=c(mu=1/50,iota=10,pop=50e6,gamma=26,Beta=400,sigma=0.1, S_0=0.07,I_0=0.001,R_0=0.93) ) -> ew1
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.