<jupyter_start><jupyter_text>Back-projection for the cases with unknown illness onset<jupyter_code>## adding extra 10 days for the stability of back-projection procedure
temp_lastdays <- matrix(NA, ncol=1, nrow=10)
temp_lastdays %<>% as.data.frame() %>% mutate(labconf = as.Date((max(unknown_Japan$labconf)+1):(max(unknown_Japan$labconf)+10)),
                                              total = 0) %>% dplyr::select(-V1)

rbind(unknown_Japan, temp_lastdays) -> unknown_Japan
unknown_Japan %<>% mutate(time_onset = 1:nrow(unknown_Japan))


## time delay from the illness onset to confirmation
K = nrow(unknown_Japan)

rep_fit=list(shape=1.4875, scale=9.289) 
report_probability = pweibull(1:K, shape=rep_fit$shape, scale=rep_fit$scale) - pweibull(1:K-1, shape=rep_fit$shape, scale=rep_fit$scale)
report_pmf = c(0,report_probability[1:21])


## back-projecton of unknown cases
sts = new("sts", epoch=unknown_Japan$time_onset, observed=unknown_Japan$total)
bpnp.control = list(k = 2, eps = rep(1e-4,2), iter.max=rep(1000,2), 
                    Tmark = nrow(sts), B = -1
                    , alpha = 0.01, verbose = FALSE, lambda0 = NULL, 
                    eq3a.method = c("R","C"))
sts_bp = backprojNP(sts, incu.pmf=report_pmf, control=modifyList(bpnp.control,list(eq3a.method="C")))
unknown_Japan$total_backproj = upperbound(sts_bp)


## normalizing the back-projected cases
unknown_Japan$total_backproj[unknown_Japan$total_backproj<=0.01] <- 0

unknown_Japan %>% mutate(total_normal = round(total_backproj/sum(total_backproj)*sum(total)),
                         time_onset=0:(nrow(unknown_Japan)-1)) %>%
dplyr::select(time_onset, labconf, total_normal) %>% rename(onset = labconf, total = total_normal) -> dt.backproj_onset


## merging the data
merge(master_df_final, dt.backproj_onset, by=c('onset'), all.y=TRUE) %>% 
mutate(total = total.x + total.y) %>% dplyr::select(onset, total) -> df_onset
df_onset[is.na(df_onset)] <- 0<jupyter_output><empty_output><jupyter_text>Likelihood<jupyter_code>date0 = as.Date("2019-12-31")
datel = finalpoint
date2num = function(date){
    as.numeric(date - date0)
}

# infection curve from growth rates
incidence = function(i0 = 10, growthrates=c(0.1,0.2,0.1), timeperiods = c(0,5,15,20), precision = 1){
    tseq = seq(timeperiods[1], timeperiods[length(timeperiods)], precision)
    ires = tseq * 0
    ibegin = i0
    ires[1] = ibegin
    for(period in 1:(length(timeperiods)-1)){
        indexbeginend = c(timeperiods[period],timeperiods[period+1])/precision + 1
        ts = seq(timeperiods[period],timeperiods[period+1], precision)
        ires[seq(indexbeginend[1],indexbeginend[2])] = ires[indexbeginend[1]]*exp(growthrates[period]*(ts-ts[1]))
    }
    return(ires)
}

# infection curve from change points
incidencefrompoints = function(i_t, timeperiods, precision = 1){
    i0 = i_t[1]
    growthrates = diff(log(i_t)) / diff(timeperiods)
    incidence(i0, growthrates, timeperiods, precision)
}

# growth rates from change points
gratefrompoints = function(i_t, timeperiods, precision = 1){
    i0 = i_t[1]
    growthrates = diff(log(i_t)) / diff(timeperiods)
    return(growthrates)
}
# reporting delay
reportingdelay = Weibull(1.4875, 9.289)

# incubation period
incubationperiod = Lnorm(1.525, 0.629)

# overall delay
overalldelay = incubationperiod
disc_overalldelay = p(overalldelay)(1:100)-p(overalldelay)(1:100-1)
# convolution
observeincidence = function(incidence, delay){
    Ires = incidence * 0
    for(t in 1:(length(Ires)-1)){
        extended_delay=c(delay,numeric(max(0,t-length(delay))))
        Ires[t+1] = sum(incidence[1:t] * extended_delay[t:1])
    }
    Ires
}

# likelihood
nllike = function(param, Data){

    i_t = param
    inc = incidencefrompoints(i_t, Data$timeperiods)
    eI_t = observeincidence(inc, Data$delay)
   
    ll = -sum(dpois(Data$I, eI_t, T))
}
# load data
data <- df_onset
data$onset <- as.Date(data$onset)

# simple nowcasting 
data %>% filter(onset == finaldata) %>% dplyr::select(X) %>% as.numeric() -> Today
delay_fit=list(shape=1.4898, scale=9.1510)
delay_pdf <- function(t){pweibull(Today-t, delay_fit$shape, delay_fit$scale)}

data %<>% mutate(cdf_incu = delay_pdf(X), nowcasted = round(total/cdf_incu)) %>%
filter(onset <= (finalpoint))

casecounts = c(numeric(date2num(as.Date(data$onset[1]))-1), data$nowcasted)<jupyter_output><empty_output><jupyter_text>Estimation<jupyter_code>library(LaplacesDemon)

## likelihood for MCMC
nllike = function(param, Data, q = 1){
    i_t = exp(param)
    inc = incidencefrompoints(i_t, Data$timeperiods)
    eI_t = observeincidence(inc, Data$delay)
    ll = 0
    ll = ll + sum(dpois(Data$I, eI_t, T))
    ll = ll - dnorm(0.1, q, 0.05)
    -ll
}


## MCMC
set.seed(19)
Data = c(list(N = 0, parm.names = c(paste0("i",0:3)), mon.names = "LP"), Data)
MCMC = function(par, Data){
    LL = -nllike(par, Data)
    list(LP = LL, Dev = -2*LL, Monitor = LL, yhat=NULL, parm = par)
}

options(warn=-1)
set.seed(32433)
fit = LaplacesDemon(MCMC, Data, Initial.Values = numeric(4), Covar = NULL, Iterations = 200000, Status = 100000, Thinning = 100)
options(warn=0)
fit_d = deburn(fit, fit$Thinned.Samples/2)
i_t = exp(fit_d$Posterior1)
rs = apply(exp(fit_d$Posterior1),1,gratefrompoints,timeperiods=Data$timeperiods)
rownames(rs) = paste0("r",1:3)
# i_t CrIs
timeperiods = date2num(c(date0, as.Date(c("2020-3-25", "2020-4-7")),datel))
incidences = apply(i_t,1,incidencefrompoints, timeperiods = timeperiods)
infectioncurves = apply(incidences, 1, quantile, probs = c(0.025,0.5,0.975))

onsets = apply(incidences ,2, observeincidence, delay = disc_overalldelay)
onsetcurves = apply(onsets, 1, quantile, probs = c(0.025,0.5,0.975))<jupyter_output><empty_output><jupyter_text>Projected incidence from the last data point to 25 May 2020<jupyter_code>iter = 1000
i_projection1 <- list()
(as.numeric(finalemer - datel)) -> nubmerrow
for (g in 1:iter){i_t[,4][g]*exp((rs[3,][g])*(1:nubmerrow)) -> i_projection1[[g]]}
matrix(unlist(i_projection1),ncol=iter,byrow=FALSE) -> i_projection1_matrix

inci_proj <- matrix(NA, nrow=nubmerrow, ncol=3)

for (j in 1:nubmerrow){
    inci_proj[j,3] <- as.numeric(quantile(i_projection1_matrix[j,1:iter], .975, na.rm = FALSE))
    inci_proj[j,2] <- as.numeric(quantile(i_projection1_matrix[j,1:iter], .5, na.rm = FALSE))
    inci_proj[j,1] <- as.numeric(quantile(i_projection1_matrix[j,1:iter], .025, na.rm = FALSE))}

cbind(infectioncurves, t(inci_proj)) -> infectioncurves_projected<jupyter_output><empty_output><jupyter_text>Projected reported case from the last data point to 25 May 2020<jupyter_code>rbind(incidences, i_projection1_matrix) -> incidences_projected
onsets_projected = apply(incidences_projected ,2, observeincidence, delay = disc_overalldelay)
onsetcurves_projected = apply(onsets_projected, 1, quantile, probs = c(0.025,0.5,0.975))
matplot(t(onsetcurves_projected), type = "l", col = 4, lty = c(2,1,2))<jupyter_output><empty_output><jupyter_text>Computation of NGMs via Julia<jupyter_code>library(JuliaCall)

julia_setup(JULIA_HOME = "C:/Users/##USER##/AppData/Local/Programs/Julia/Julia-1.4.1/bin")
julia_library("Optim, CSV, LinearAlgebra, DelimitedFiles, Distributions")
julia_command("
function KLdiv(dist_p, dist_q, kld = false, upper = 300)
    div = 0
    for i in 0:upper
        if kld
            div += pdf(dist_p,i) * (logpdf(dist_p,i) - logpdf(dist_q, i))
        else
            div += (pdf(dist_p,i) - pdf(dist_q,i))^2
        end
    end
    div
end")
julia_command("MixtureGeomGeom(lambda, mu, mix) = MixtureModel([Geometric(1/(1+lambda)),Geometric(1/(1+mu))], [mix, 1.0-mix])")
julia_command("NBmu(mu, k) = NegativeBinomial(k, k/(k+mu))")
julia_command("rtoR0(r, SI = 4.8, CV = (2.3/4.8)) = (1 + r*SI*CV^2)^(1/(CV^2))")
julia_command("function NGM(Re, lambda, mu, ph)
    mixprop = (Re - mu) / (lambda - mu)
    eqmat = [1 mu-lambda; 0 ph-mixprop]
    rho_b = inv(eqmat) * [mu, (1 - mixprop) * ph]
    rho = rho_b[1]; b = rho_b[2]
    V = [mixprop b; 1-mixprop 1-b]
    V*diagm([Re, rho])*inv(V)
end")
julia_command("function R0toNGM(R0, kldvalues, k = 0.1, phighrisk = 0.99)
    targetNB = NBmu(R0, k)
    function objfunc(params, targetdist = targetNB, kld = true)
        logdiv = log(KLdiv(targetdist, MixtureGeomGeom(params...), kld))
        logdiv
    end
    opt = optimize(objfunc, [1e-4,R0,0.01], [R0,50,0.99], [0.1,5,0.1])
    pars = opt.minimizer
    kldvalues.=opt.minimum
    NGM(R0, pars[1],pars[2], phighrisk)
end")
julia_assign('rsamples',t(rs))
julia_command('kldvalues=[[0.0] for x in rsamples]')
julia_command("R0s = rtoR0.(rsamples);")
julia_command("NGMs = R0toNGM.(R0s,kldvalues);")
julia_command("kldvals=first.(kldvalues)")
kldoutput=julia_eval("kldvals")
# Convert NGMs to column major vectors and save as csv
julia_command("NGMmat=reshape.(NGMs,Ref((4,1)));")
julia_command("NGMoutput = reduce(hcat, [vcat(i...) for i in (eachrow(NGMmat))]);")
NGMoutput = julia_eval("NGMoutput")<jupyter_output><empty_output><jupyter_text>95% CIs of NGMs<jupyter_code>SI = 4.8
CV = 2.3/SI

NGMoutput -> NGMs
NGM1 <- NGMs[1:4,]; NGM2 <- NGMs[5:8,]; NGM3 <- NGMs[9:12,]
upper_r1 <- matrix(NA, nrow=4, ncol=1); median_r1 <- matrix(NA, nrow=4, ncol=1); lower_r1 <- matrix(NA, nrow=4, ncol=1)
upper_r2 <- matrix(NA, nrow=4, ncol=1); median_r2 <- matrix(NA, nrow=4, ncol=1); lower_r2 <- matrix(NA, nrow=4, ncol=1)
upper_r3 <- matrix(NA, nrow=4, ncol=1); median_r3 <- matrix(NA, nrow=4, ncol=1); lower_r3 <- matrix(NA, nrow=4, ncol=1)

for (j in 1:4){
    upper_r1[j,1] <- as.numeric(quantile(NGM1[j,1:iter], .975, na.rm = FALSE))
    median_r1[j,1] <- as.numeric(quantile(NGM1[j,1:iter], .5, na.rm = FALSE))
    lower_r1[j,1] <- as.numeric(quantile(NGM1[j,1:iter], .025, na.rm = FALSE))

    upper_r2[j,1] <- as.numeric(quantile(NGM2[j,1:iter], .975, na.rm = FALSE))
    median_r2[j,1] <- as.numeric(quantile(NGM2[j,1:iter], .5, na.rm = FALSE))
    lower_r2[j,1] <- as.numeric(quantile(NGM2[j,1:iter], .025, na.rm = FALSE))

    upper_r3[j,1] <- as.numeric(quantile(NGM3[j,1:iter], .975, na.rm = FALSE))
    median_r3[j,1] <- as.numeric(quantile(NGM3[j,1:iter], .5, na.rm = FALSE))
    lower_r3[j,1] <- as.numeric(quantile(NGM3[j,1:iter], .025, na.rm = FALSE))}

(cbind(upper_r1, median_r1, lower_r1, 
       upper_r2, median_r2, lower_r2,
       upper_r3, median_r3, lower_r3) -> NGM_CIs)<jupyter_output><empty_output><jupyter_text>Incidence from NGMs<jupyter_code>## decrease from NGM of Period 1
NGM_30 <- rbind(NGM1[1:2,], (NGM1[3:4,]*0.5))
NGM_50 <- rbind(NGM1[1:2,], (NGM1[3:4,]*0.7))
NGM_70 <- rbind(NGM1[1:2,], (NGM1[3:4,]*0.9))
## growth rate for each scenarios

r1 <-  matrix(NA, nrow=iter, ncol=1); r2 <-  matrix(NA, nrow=iter, ncol=1); r3 <-  matrix(NA, nrow=iter, ncol=1);
r4 <-  matrix(NA, nrow=iter, ncol=1); r5 <-  matrix(NA, nrow=iter, ncol=1); r6 <-  matrix(NA, nrow=iter, ncol=1);

for (k in 1:iter){
    r1[k,1] <- ((eigen(matrix(NGM1[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
    r2[k,1] <- ((eigen(matrix(NGM2[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
    r3[k,1] <- ((eigen(matrix(NGM3[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
    
    r4[k,1] <- ((eigen(matrix(NGM_30[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
    r5[k,1] <- ((eigen(matrix(NGM_50[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
    r6[k,1] <- ((eigen(matrix(NGM_70[,k],nrow=2))$value[1])^(CV^2)-1)/(SI*CV^2)
}

## quantile of growth rates
r1_quan <- as.numeric(c(quantile(r1, 0.975), quantile(r1, 0.5), quantile(r1, 0.025)))
r2_quan <- as.numeric(c(quantile(r2, 0.975), quantile(r2, 0.5), quantile(r2, 0.025)))
r3_quan <- as.numeric(c(quantile(r3, 0.975), quantile(r3, 0.5), quantile(r3, 0.025)))

r4_quan <- as.numeric(c(quantile(r4, 0.975), quantile(r4, 0.5), quantile(r4, 0.025)))
r5_quan <- as.numeric(c(quantile(r5, 0.975), quantile(r5, 0.5), quantile(r5, 0.025)))
r6_quan <- as.numeric(c(quantile(r6, 0.975), quantile(r6, 0.5), quantile(r6, 0.025)))
## projection of newly incidence

(Final.T = as.numeric(as.Date("2020-07-30")-finalemer))
(i0 = as.numeric(infectioncurves_projected[2,ncol(infectioncurves_projected)]))

case1 <- matrix(NA, nrow= Final.T, ncol=3); case2 <- matrix(NA, nrow= Final.T, ncol=3); case3 <- matrix(NA, nrow= Final.T, ncol=3)
case4 <- matrix(NA, nrow= Final.T, ncol=3); case5 <- matrix(NA, nrow= Final.T, ncol=3); case6 <- matrix(NA, nrow= Final.T, ncol=3)
for (g in 1:3){
    case1[,g] <- as.matrix(i0*exp(r1_quan[g]*(1:Final.T)))
    case2[,g] <- as.matrix(i0*exp(r2_quan[g]*(1:Final.T)))
    case3[,g] <- as.matrix(i0*exp(r3_quan[g]*(1:Final.T)))
    
    case4[,g] <- as.matrix(i0*exp(r4_quan[g]*(1:Final.T)))
    case5[,g] <- as.matrix(i0*exp(r5_quan[g]*(1:Final.T)))
    case6[,g] <- as.matrix(i0*exp(r6_quan[g]*(1:Final.T)))}
colnames(case1) <- c("upper", "med", "lower"); colnames(case2) <- c("upper", "med", "lower"); colnames(case3) <- c("upper", "med", "lower")
colnames(case4) <- c("upper", "med", "lower"); colnames(case5) <- c("upper", "med", "lower"); colnames(case6) <- c("upper", "med", "lower")

case1 %>% as.data.frame() %>% mutate(time = 1:nrow(case1)) %>% mutate(cal_date = finalemer+time) -> case1_figure
case2 %>% as.data.frame() %>% mutate(time = 1:nrow(case2)) %>% mutate(cal_date = finalemer+time) -> case2_figure
case3 %>% as.data.frame() %>% mutate(time = 1:nrow(case3)) %>% mutate(cal_date = finalemer+time) -> case3_figure

case4 %>% as.data.frame() %>% mutate(time = 1:nrow(case4)) %>% mutate(cal_date = finalemer+time) -> case4_figure
case5 %>% as.data.frame() %>% mutate(time = 1:nrow(case5)) %>% mutate(cal_date = finalemer+time) -> case5_figure
case6 %>% as.data.frame() %>% mutate(time = 1:nrow(case6)) %>% mutate(cal_date = finalemer+time) -> case6_figure

t(infectioncurves_projected) %>% as.data.frame() -> basic_line
basic_line %>% mutate(time = 1:nrow(basic_line)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> basic_line
colnames(basic_line) <- c("upper", "med", "lower", "time", "cal_date")

basic_line %>% tail(1) -> temp_start
rbind(temp_start, case1_figure) -> case1_figure_connect
rbind(temp_start, case2_figure) -> case2_figure_connect
rbind(temp_start, case3_figure) -> case3_figure_connect

rbind(temp_start, case4_figure) -> case4_figure_connect
rbind(temp_start, case5_figure) -> case5_figure_connect
rbind(temp_start, case6_figure) -> case6_figure_connect
options(warn=-1)
options(repr.plot.width=10,repr.plot.height=5)

range = c(as.Date("2020-02-01"), as.Date("2020-07-30"))

case1_figure %>% ggplot() + 
geom_line(data=case1_figure_connect,aes(x=cal_date,y=med,colour="Scenario1 (back to normal)"),size=1) +
geom_ribbon(data=case1_figure_connect,aes(x=cal_date, ymax=case1_figure_connect$upper, ymin=case1_figure_connect$lower),fill="tomato", alpha = 0.4) +

geom_line(data=case4_figure_connect,aes(x=cal_date,y=med,colour="Scenario4 (50% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=case4_figure_connect,aes(x=cal_date, ymax=case4_figure_connect$upper, ymin=case4_figure_connect$lower),fill="salmon3", alpha = 0.4) +

geom_line(data=case5_figure_connect,aes(x=cal_date,y=med,colour="Scenario3 (30% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=case5_figure_connect,aes(x=cal_date, ymax=case5_figure_connect$upper, ymin=case5_figure_connect$lower),fill="violetred3", alpha = 0.4) +

geom_line(data=case6_figure_connect,aes(x=cal_date,y=med,colour="Scenario2 (10% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=case6_figure_connect,aes(x=cal_date, ymax=case6_figure_connect$upper, ymin=case6_figure_connect$lower),fill="seagreen4", alpha = 0.4) +

geom_line(data=basic_line,aes(x=cal_date,y=med),color="gray4",size=1) +
geom_ribbon(data=basic_line,aes(x=cal_date, ymax=basic_line$upper, ymin=basic_line$lower),fill="gray4", alpha = 0.4) +
labs(x="\nDate of infection", y="Number of new infections\n") +
scale_x_date(date_labels="%m/%d",date_breaks  ="3 weeks", limits = range, expand = c(0, 0)) +
theme(text = element_text(size=12, family="sans",color="black"),
      axis.text = element_text(size=10, family="sans",color="black"),
      panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
ggtitle("Tokyo") +
geom_vline(xintercept=as.Date("2020-03-26"), linetype="dashed", color = "tomato", size =0.7) +
geom_vline(xintercept=as.Date("2020-04-08"), linetype="dashed", color = "tomato", size =0.7) +
geom_vline(xintercept=as.Date("2020-05-25"), linetype="dashed", color = "tomato", size =0.7) +
scale_y_continuous(expand = c(0, 0)) + coord_cartesian(ylim = c(0:250)) +
scale_colour_manual("", values = c("Scenario1 (back to normal)"="tomato", 
                                   "Scenario2 (10% decrease in high-risk transmission)"="seagreen4", 
                                   "Scenario3 (30% decrease in high-risk transmission)"="violetred3",
                                   "Scenario4 (50% decrease in high-risk transmission)"="#1380A1"))<jupyter_output><empty_output><jupyter_text>Observed cases from NGMs<jupyter_code>Final.TT <- length(as.vector(c(infectioncurves_projected[2,])))+nrow(case1)

## iterated incidience 
case1 <- matrix(NA, nrow= Final.T, ncol=iter); case2 <- matrix(NA, nrow= Final.T, ncol=iter); case3 <- matrix(NA, nrow= Final.T, ncol=iter)
case4 <- matrix(NA, nrow= Final.T, ncol=iter); case5 <- matrix(NA, nrow= Final.T, ncol=iter); case6 <- matrix(NA, nrow= Final.T, ncol=iter)


for (g in 1:iter){
    case1[,g] <- as.matrix(i0*exp(r1[g]*(1:Final.T)))
    case2[,g] <- as.matrix(i0*exp(r2[g]*(1:Final.T)))
    case3[,g] <- as.matrix(i0*exp(r3[g]*(1:Final.T)))
    
    case4[,g] <- as.matrix(i0*exp(r4[g]*(1:Final.T)))
    case5[,g] <- as.matrix(i0*exp(r5[g]*(1:Final.T)))
    case6[,g] <- as.matrix(i0*exp(r6[g]*(1:Final.T)))
}

## iterated incidience from each scenarios * incubation period

onset1 <- matrix(NA, nrow= Final.TT, ncol=iter)
onset2 <- matrix(NA, nrow= Final.TT, ncol=iter)
onset3 <- matrix(NA, nrow= Final.TT, ncol=iter)

onset4 <- matrix(NA, nrow= Final.TT, ncol=iter)
onset5 <- matrix(NA, nrow= Final.TT, ncol=iter)
onset6 <- matrix(NA, nrow= Final.TT, ncol=iter)

for (g in 1:iter){
    onset1[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case1[,g]))), Data$delay))
    onset2[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case2[,g]))), Data$delay))
    onset3[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case3[,g]))), Data$delay))

    onset4[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case4[,g]))), Data$delay))
    onset5[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case5[,g]))), Data$delay))
    onset6[,g] <- as.matrix(observeincidence(as.vector(c(infectioncurves_projected[2,],as.vector(case6[,g]))), Data$delay))
}

## calculating quantiles

onset_ci1 <- matrix(NA, nrow=nrow(onset1), ncol=3)
onset_ci2 <- matrix(NA, nrow=nrow(onset1), ncol=3)
onset_ci3 <- matrix(NA, nrow=nrow(onset1), ncol=3)

onset_ci4 <- matrix(NA, nrow=nrow(onset1), ncol=3)
onset_ci5 <- matrix(NA, nrow=nrow(onset1), ncol=3)
onset_ci6 <- matrix(NA, nrow=nrow(onset1), ncol=3)

for (j in 1:nrow(onset1)){
    onset_ci1[j,1] <- as.numeric(quantile(onset1[j,1:iter], .975, na.rm = FALSE))
    onset_ci1[j,2] <- as.numeric(quantile(onset1[j,1:iter], .5, na.rm = FALSE))
    onset_ci1[j,3] <- as.numeric(quantile(onset1[j,1:iter], .025, na.rm = FALSE))

    onset_ci2[j,1] <- as.numeric(quantile(onset2[j,1:iter], .975, na.rm = FALSE))
    onset_ci2[j,2] <- as.numeric(quantile(onset2[j,1:iter], .5, na.rm = FALSE))
    onset_ci2[j,3] <- as.numeric(quantile(onset2[j,1:iter], .025, na.rm = FALSE))

    onset_ci3[j,1] <- as.numeric(quantile(onset3[j,1:iter], .975, na.rm = FALSE))
    onset_ci3[j,2] <- as.numeric(quantile(onset3[j,1:iter], .5, na.rm = FALSE))
    onset_ci3[j,3] <- as.numeric(quantile(onset3[j,1:iter], .025, na.rm = FALSE))
    
    onset_ci4[j,1] <- as.numeric(quantile(onset4[j,1:iter], .975, na.rm = FALSE))
    onset_ci4[j,2] <- as.numeric(quantile(onset4[j,1:iter], .5, na.rm = FALSE))
    onset_ci4[j,3] <- as.numeric(quantile(onset4[j,1:iter], .025, na.rm = FALSE))
    
    onset_ci5[j,1] <- as.numeric(quantile(onset5[j,1:iter], .975, na.rm = FALSE))
    onset_ci5[j,2] <- as.numeric(quantile(onset5[j,1:iter], .5, na.rm = FALSE))
    onset_ci5[j,3] <- as.numeric(quantile(onset5[j,1:iter], .025, na.rm = FALSE))
    
    onset_ci6[j,1] <- as.numeric(quantile(onset6[j,1:iter], .975, na.rm = FALSE))
    onset_ci6[j,2] <- as.numeric(quantile(onset6[j,1:iter], .5, na.rm = FALSE))
    onset_ci6[j,3] <- as.numeric(quantile(onset6[j,1:iter], .025, na.rm = FALSE))
}
colnames(onset_ci1) <- c("upper", "med", "lower")
colnames(onset_ci2) <- c("upper", "med", "lower")
colnames(onset_ci3) <- c("upper", "med", "lower")

colnames(onset_ci4) <- c("upper", "med", "lower")
colnames(onset_ci5) <- c("upper", "med", "lower")
colnames(onset_ci6) <- c("upper", "med", "lower")

onset_ci1 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci1_figure
onset_ci2 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci2_figure
onset_ci3 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci3_figure

onset_ci4 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci4_figure
onset_ci5 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci5_figure
onset_ci6 %>% as.data.frame() %>% mutate(time = 1:nrow(onset_ci1)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> onset_ci6_figure

t(onsetcurves_projected) %>% as.data.frame() -> basic_line2
basic_line2 %>% mutate(time = 1:nrow(basic_line)) %>% mutate(cal_date = as.Date("2019-12-30")+time) -> basic_line2
colnames(basic_line2) <- c("upper", "med", "lower", "time", "cal_date")
options(warn=-1)

data %>% ggplot()+ 
geom_bar(aes(x=onset, y=nowcasted), stat='identity', fill="#FAAB18",size=0.7) +
geom_line(data=onset_ci1_figure,aes(x=cal_date,y=med,colour ="Scenario1 (back to normal)"),size=1) +
geom_ribbon(data=onset_ci1_figure,aes(x=cal_date, ymax=onset_ci1_figure$upper, ymin=onset_ci1_figure$lower),fill="tomato", alpha = 0.4) +

geom_line(data=onset_ci4_figure,aes(x=cal_date,y=med,color="Scenario4 (50% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=onset_ci4_figure,aes(x=cal_date, ymax=onset_ci4_figure$upper, ymin=onset_ci4_figure$lower),fill="#1380A1", alpha = 0.4) +

geom_line(data=onset_ci5_figure,aes(x=cal_date,y=med,colour ="Scenario3 (30% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=onset_ci5_figure,aes(x=cal_date, ymax=onset_ci5_figure$upper, ymin=onset_ci5_figure$lower),fill="violetred3", alpha = 0.4) +

geom_line(data=onset_ci6_figure,aes(x=cal_date,y=med,colour ="Scenario2 (10% decrease in high-risk transmission)"),size=1) +
geom_ribbon(data=onset_ci6_figure,aes(x=cal_date, ymax=onset_ci6_figure$upper, ymin=onset_ci6_figure$lower),fill="seagreen4", alpha = 0.4) +

geom_line(data=basic_line2,aes(x=cal_date,y=med),color="gray4",size=1) +
geom_ribbon(data=basic_line2,aes(x=cal_date, ymax=basic_line2$upper, ymin=basic_line2$lower),fill="gray4", alpha = 0.4) +
labs(x="\nDate of illness onset", y="Number of new cases\n") +
scale_x_date(date_labels="%m/%d",date_breaks  ="3 weeks",limits = range, expand = c(0, 0)) +
theme(text = element_text(size=12, family="sans",color="black"),
      axis.text = element_text(size=10, family="sans",color="black"),
      panel.grid.major = element_blank(), panel.grid.minor = element_blank()) +
ggtitle("Tokyo") +
geom_vline(xintercept=as.Date("2020-03-26"), linetype="dashed", color = "tomato", size =1) +
geom_vline(xintercept=as.Date("2020-04-08"), linetype="dashed", color = "tomato", size =1) +
geom_vline(xintercept=as.Date("2020-05-25"), linetype="dashed", color = "tomato", size =0.7) +
scale_y_continuous(expand = c(0, 0)) + coord_cartesian(ylim = c(0:250)) +
scale_colour_manual("", values = c("Scenario1 (back to normal)"="tomato", 
                                   "Scenario2 (10% decrease in high-risk transmission)"="seagreen4", 
                                   "Scenario3 (30% decrease in high-risk transmission)"="violetred3",
                                   "Scenario4 (50% decrease in high-risk transmission)"="#1380A1"))<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.19405180161777924, 'fraction_numerical': 0.051924176811831035, 'mean_word_length': 5.493103448275862, 'pattern_counts': {'":': 0, '<': 170, '<?xml version=': 0, '>': 130, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '9', 'dir': 'jupyter-structured-clean-dedup', 'id': '289111', 'license': 'MIT', 'n_tokens_mistral': 11048, 'n_tokens_neox': 10110, 'n_words': 1564, 'path': 'projection.ipynb', 'repo_name': 'SungmokJung/Projection_Japan_COVID19', 'size': '24478'}