Nothing
#' Estimate Death Rate
#'
#' This function provides a first estimate of the death rate, to faciliate
#' convergence of the main algorithm.
#' @param tTarget time vector
#' @param Q target time-histories of the quarantined cases
#' @param D target time-histories of the dead cases
#' @param guess initial guess parameters for kappa
#' @param ftol nls.lm.control object. non-negative numeric. Default is \code{1e-6}
#' @param ptol nls.lm.control object. non-negative numeric. Default is \code{1e-6}
#' @param gtol nls.lm.control object. non-negative numeric. Default is \code{1e-6}
#' @param epsfcn nls.lm.control object. Default is \code{0.001}
#' @param factor nls.lm.control object. Default is \code{100}
#' @param maxfev nls.lm.control object. Default is \code{1000}
#' @param maxiter nls.lm.control object. Default is \code{100}
#' @param nprint nls.lm.control object. Default is \code{1}
#' @param trace Set \code{TRUE} to trace iteration results
#'
#' @importFrom nlsr nlxb
#'
#' @author Selcuk Korkmaz, \email{selcukorkmaz@gmail.com}
#'
#' @return vector of estimation and optimization function for the death rate
#'
#' @seealso \code{\link{SEIQRDP}} \code{\link{fit_SEIQRDP}}
#'
#' @references Peng, L., Yang, W., Zhang, D., Zhuge, C., Hong, L. 2020. “Epidemic analysis of COVID-19 in China by dynamical modeling”, arXiv preprint arXiv:2002.06563.
#' @references \url{https://www.mathworks.com/matlabcentral/fileexchange/74545-generalized-seir-epidemic-model-fitting-and-computation}
getKappaFun <- function(tTarget, Q, D, guess, ftol,
ptol, gtol, epsfcn, factor, maxfev,
maxiter, nprint, trace){
if (max(D)<10){
kappaFun = function(a,t) {a[1] / (exp(a[2]*(t-a[3])))}
}else{
myFun1 = function(a,t) {a[1] / (exp(a[2]*(t-a[3])) + exp(-a[2]*(t-a[3])))}
myFun2 = function(a,t) {a[1]*exp(-(a[2]*(t-a[3]))^2)}
myFun3 = function(a,t) {a[1] + exp(-a[2]*(t+a[3]))}
rate = as.numeric((diff(as.numeric(D))/median(diff(tTarget)))/Q[2:length(Q)])
x = tTarget[2:length(tTarget)]
rate[abs(rate)>3]=NA
if (length(rate[rate==0])/length(rate) < 0.5){
rate[abs(rate)==0] = NA
}
df = cbind.data.frame(tk = x[!is.na(rate)], z = rate[!is.na(rate)])
# ctrl = nls.lm.control(ftol = ftol,
# ptol = ptol, gtol = gtol, diag = list(), epsfcn = epsfcn,
# factor = factor, maxfev = maxfev, maxiter = maxiter, nprint = nprint)
ctrl = list(phi=1, lamda = 0.0001, offset = 100, laminc=10, lamdec = 4)
lsqCurveFit1 <- nlxb(z ~ (a1 / (exp(a2*(tk-a3)) + exp(-a2*(tk-a3)))),
start=c(a1=guess[[8]], a2=guess[[9]], a3=guess[[10]]), data=df, trace= trace,
control = ctrl, lower = c(0,0,0), upper = c(1,1,100))
coeff1 = lsqCurveFit1$coefficients
r1 = lsqCurveFit1$ssquares
ctrl = list(phi=1, lamda = 0.0001, offset = 100, laminc=10, lamdec = 4)
lsqCurveFit2 <- nlxb(z ~ a1*exp(-(a2*(tk-a3))^2),
start=list(a1=guess[[8]], a2=guess[[9]], a3=guess[[10]]), data=df, trace=trace,
control = ctrl, lower = c(0,0,0), upper = c(1,1,100))
coeff2 = lsqCurveFit2$coefficients
r2 = lsqCurveFit2$ssquares
lsqCurveFit3 <- nlxb(z ~ a1 + exp(-a2*(tk+a3)),
start=list(a1=guess[[8]], a2=guess[[9]], a3=guess[[10]]), data=df, trace=trace,
control = ctrl, lower = c(0,0,0), upper = c(1,1,100))
coeff3 = lsqCurveFit3$coefficients
r3 = lsqCurveFit3$ssquares
minR = min(r1,r2,r3)
if (r1==minR){
kappaGuess = coeff1
kappaFun = myFun1
}else if (r2==minR){
kappaGuess = coeff2
kappaFun = myFun2
}else if (r3==minR){
kappaFun = myFun3
kappaGuess = coeff3
}
guess[8:10] = kappaGuess
}
return(list(guess=guess, kappaFun=kappaFun))
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.