Nothing
nls.tv <- function(pars, kappa, active, basisvals, times, data, coefs, lik, proc, start, X.index, control.out, control.in, in.meth){
if(control.out$method == "twoStage"){
delta <- rep(1, length(data))
}
pars.names <- names(pars)
f.conv <- pars.kappa <- c()
maxStep <- 10
lambda.sparse <- control.out$lambda.sparse
for(i in 1:control.out$maxIter){
for(j in 1:maxStep){
linObj <- ProfileSSE.AllPar.tv(pars = c(pars,kappa), times = times, data = data, coefs = coefs, lik = lik, proc = proc, in.meth = in.meth,control.in = control.in)
f.new <- linObj$f
f.new <- sum(f.new^2)
if(control.out$method == "penalized"){
f.new <- f.new + lambda.sparse * sum(abs(kappa[-1] - kappa[-length(kappa)]))
}
if(control.out$echo){
print(x = c(paste("Iter:", i, f.new)))
cat(pars, kappa, "\n")
}
if(i == 1){
break
}else{
if(f.conv[i - 1] - f.new > 0 & f.conv[i - 1] - f.new < control.out$tol){
return(list(pars=pars.old, kappa = kappa.old, coefs = coefs, f = f.new, y = y, Xdf = Xdf, Zdf = Zdf, conv = list(f = f.conv, pars.kappa=pars.kappa, conv.message = "Converged.")))
}
if(f.conv[i - 1] - f.new > 0){
break
}
if(f.conv[i - 1] - f.new < 0 & j == maxStep){
return(list(pars=pars.old, kappa = kappa.old, coefs = coefs, f = f.new, y = y, Xdf = Xdf, Zdf = Zdf, conv = list(f = f.conv, pars.kappa=pars.kappa, conv.message = "Non-dereasing objective.")))
}
pars <- 0.5*(pars - pars.old) + pars.old
kappa <- 0.5*(kappa - kappa.old) + kappa.old
}
}
pars.old <- pars
kappa.old <- kappa
f.conv <- c(f.conv, f.new)
pars.kappa <- rbind(pars.kappa, c(pars, kappa))
Xdf <- - linObj$df[, 1:length(pars), drop = FALSE]
Zdf <- - linObj$df[, (length(pars) +1): dim(linObj$df)[2]]
y <- - linObj$df %*% c(pars, kappa) + linObj$f
coefs <- linObj$coefs
if(control.out$method == "penalized"){
res <- penalized::penalized(response = y, penalized = Zdf, unpenalized = Xdf, lambda2 = lambda.sparse, positive = TRUE, fusedl = TRUE, trace = FALSE)
pars <- res@unpenalized
kappa <- res@penalized
}
if(control.out$method == "nnls"){
res <- nnls::nnls(A = cbind(Xdf, Zdf), b= y)
pars <- res$x[1:length(pars)]
kappa <- res$x[(length(pars) + 1) : length(res$x)]
}
names(pars) <- pars.names
}
return(list(pars=pars.old, kappa = kappa.old, coefs = coefs, f = f.new, y = y, Xdf = Xdf, Zdf = Zdf, conv = list(f = f.conv, pars.kappa=pars.kappa, conv.message = "Maximum iterations reached.")))
}
Profile.LS.tv <- function(fn, data, times, pars, kappa, coefs = NULL, basisvals = NULL,
lambda, fd.obj = NULL, more = NULL, weights = NULL, quadrature = NULL,
in.meth = "nlminb", out.meth = "nls", control.in = list(),
control.out = list(), eps = 1e-06, active = NULL, posproc = FALSE,
poslik = FALSE, discrete = FALSE, names = NULL, sparse = FALSE,
likfn = make.id(), likmore = NULL)
{
if (is.null(active)) {
active = 1:length(c(pars, kappa))
}
profile.obj = LS.setup(pars = c(pars, kappa), coefs = coefs, fn = fn,
basisvals, lambda = lambda, fd.obj, more, data, weights,
times, quadrature, eps = 1e-06, posproc, poslik, discrete,
names, sparse, likfn = make.id(), likmore = NULL)
dims = dim(data)
lik = profile.obj$lik
proc = profile.obj$proc
proc$more$more$nKappa <- length(kappa)
coefs = profile.obj$coefs
data = profile.obj$data
times = profile.obj$times
Ires <- inneropt(data, times, c(pars, kappa), coefs, lik, proc, in.meth, control.in)
## Ires <- IresTmp
ncoefs <- Ires$coefs
apars = c(pars, kappa)[active]
aparamnames = names(apars)
if (is.null(control.out$maxIter)) {
control.out$maxIter = 100
}
if (is.null(control.out$tol)){
control.out$tol = 1e-08
}
res <- nls.tv(pars = pars, kappa = kappa, active = active, basisvals = basisvals, times = times, data = data, coefs = ncoefs, lik = lik, proc = proc, control.out = control.out, control.in = control.in, in.meth = in.meth)
ncoefs <- res$coefs
return(list( data = data,res = res, ncoefs = ncoefs))
}
LS.tv <- function(fn, data, times, pars, kappa, coefs = NULL, basisvals = NULL,
lambda, fd.obj = NULL, more = NULL, weights = NULL, quadrature = NULL,
in.meth = "nlminb", out.meth = "nls", control.in = list(),
control.out = list(), eps = 1e-06, active = NULL, posproc = FALSE,
poslik = FALSE, discrete = FALSE, names = NULL, sparse = FALSE,
likfn = make.id(), likmore = NULL, nnls.res)
{
if (is.null(active)) {
active = 1:length(c(pars, kappa))
}
profile.obj = LS.setup(pars = c(pars, kappa), coefs = coefs, fn = fn,
basisvals, lambda = lambda, fd.obj, more, data, weights,
times, quadrature, eps = 1e-06, posproc, poslik, discrete,
names, sparse, likfn = make.id(), likmore = NULL)
dims = dim(data)
lik = profile.obj$lik
proc = profile.obj$proc
proc$more$more$nKappa <- length(kappa)
coefs = profile.obj$coefs
data = profile.obj$data
times = profile.obj$times
## Ires <- IresTmp
apars = c(pars, kappa)[active]
aparamnames = names(apars)
if (is.null(control.out$maxIter)) {
control.out$maxIter = 100
}
if (is.null(control.out$tol)){
control.out$tol = 1e-08
}
res <- nnls.res
ncoefs <- res$coefs
if(is.null(control.out$pars.c))
control.out$pars.c <- 100
if(control.out$lambda.sparse == -1){
Zdf <- res$Zdf
Xdf <- res$Xdf
y <- res$y
pars <- res$pars
Z0 <- rowSums(Zdf)
beta0 <- sum((y - Xdf %*% pars ) * Z0) / sum(Z0^2)
lambda0 <- max(abs(as.vector(t(y -Xdf %*% pars - Z0 * beta0) %*% sweep(Zdf, 1, mean(Zdf)))))
lambda = exp(seq(log(lambda0), log(lambda0 * 0.01), len = 50))
pars.pen <- kappa.pen <- coefs.pen <- list()
bic <- f <- rep(NA, length(lambda))
for(i in 1:length(lambda)){
lambda.sparse <- lambda[i]
res.sparse <- penalized::penalized(response = y, penalized = Zdf, unpenalized = Xdf, lambda2 = lambda.sparse, fusedl = TRUE, positive = TRUE, trace = FALSE)
pars.pen[[i]] <- res.sparse@unpenalized
kappa.pen[[i]] <- res.sparse@penalized
Ires <- inneropt(data, times, pars = c(pars.pen[[i]],kappa.pen[[i]]), ncoefs, lik, proc, in.meth, control.in)
devals <- as.matrix(lik$bvals%*%Ires$coefs)
f <- as.vector(as.matrix(data - lik$more$fn(times, devals, pars, lik$more$more)))
coefs.pen[[i]] <- Ires$coefs
sd.pen <- sd(f)
ll.pen <- - sum(f^2) / (sd.pen^2) / 2 - length(f) * log(sd.pen)
bic[i] <- -2 * ll.pen + (length(unique(kappa.pen[[i]])) + length(pars.pen[[i]])) * log(length(data))
}
i.select <- which(bic == min(bic))
sel.res <- list(pars.pen = pars.pen[[i.select]], kappa.pen = kappa.pen[[i.select]], bic = bic[i.select], coefs.pen = coefs.pen[[i.select]], lambda = lambda[i.select])
}
return(list(select = sel.res))
}
ProfileSSE.AllPar.tv <- function(pars, times, data, coefs, lik, proc,
in.meth='nlminb', control.in=NULL,
dcdp=NULL, oldpars=NULL, use.nls=TRUE, sgn=1)
{
## Squared Error outer criterion
## coefs = as.vector(coefs) # First run the inner optimization
##
## f1 = SplineCoefsErr.DDE(coefs,times,data,lik,proc,pars, kappa, basisvals = basisvals, fdobj0 = fdobj0)
## The coefficients for delayed times:
##################################################
## Not Sure:
## if(!is.null(dcdp)){
## tcoefs = as.vector(coefs) + dcdp%*%(pars-oldpars);
## f2 = SplineCoefsErr(tcoefs,times,data,lik,proc,pars)
## if(f2 < f1){
## coefs = tcoefs
## f1 = f2
## }
## }
###################################################
## Inner optimization need to be changed as well.
Ires = inneropt(data,times,pars, coefs,lik,proc, in.meth,control.in)
ncoefs = Ires$coefs
## Calculate fitted value after inner optimization:
devals = as.matrix(lik$bvals%*%ncoefs)
colnames(devals) = proc$more$names
## Squared errors: No need to change for DDE
weights = checkweights(lik$more$weights,lik$more$whichobs,data)
f = as.vector(as.matrix(data - lik$more$fn(times, devals, pars, lik$more$more))*sqrt(weights))
isnaf = is.na(f)
f[isnaf] = 0
dlikdp = lik$more$dfdp(times,devals,pars,lik$more$more)
dlikdp = matrix(dlikdp,dim(dlikdp)[1]*dim(dlikdp)[2],dim(dlikdp)[3])
## dlikdp will be zero if the likelihood doesn't directly have ode parameters,
## which is true for least square case.
dlikdx = lik$more$dfdx(times,devals,pars,lik$more$more)
## dlikdx[i,,] is an identity matrix for every i.
dlikdc = c()
for(i in 1:dim(dlikdx)[2]){
tH = c()
for(j in 1:dim(dlikdx)[3]){
## Only dlikdx[,i,i] are non-zero (all 1's)
tH = cbind(tH,as.matrix(diag(dlikdx[,i,j])%*%lik$bvals))
}
dlikdc = rbind(dlikdc,tH)
}
## ??dlikdc: why 0.5 and 0.5 ??
d2Hdc2 = SplineCoefsDC2sparse(ncoefs,times,data,lik,proc,pars)
## need not be changed?
d2Hdcdp = SplineCoefsDCDP(ncoefs, times, data, lik, proc, pars)
## Got warning message:
## In dim(weights[whichrows, ]) == dim(diffs):
## longer object length is not a multiple of shorter object length?
## Use Implicit function theorem:
if(is.matrix(d2Hdc2)){
## When will it not be a matrix? How to use solve in that case?
dcdp = MASS::ginv(d2Hdc2) %*% d2Hdcdp
} else {
dcdp = as.matrix(solve(d2Hdc2,d2Hdcdp))
}
## Chain rule:
df = dlikdc%*%dcdp + dlikdp
df[isnaf,] = 0
colnames(df) = c(proc$more$parnames)
if(!is.null(lik$report)){ print(f) }
f = sgn*f
df = sgn*df
return(list(f = f, df = df, coefs = ncoefs))
}
tvtrans <- function(t,k){
month <- t %% 1
r <- rep(0, length(t))
for(i in 1:length(k)){
r[month >= (i-1)/length(k) & month < i /length(k)] <- k[i]
}
return(r)
}
tvDSIRfn <- list()
tvDSIRfn$fn <- function (t, y, p, more)
{
r = y
nKappa <- more$nKappa
pk <- p[(length(p) - (nKappa - 1)):length(p)]
r[, "S"] = - tvtrans(t, pk) * y[,"I"] * y[, "S"] + 8000 * (sin(t / pi) / 2 + 2)
r[, "I"] = tvtrans(t, pk) * y[,"I"] * y[, "S"] - p["gamma"] * y[, "I"]
return(r)
}
tvDSIRfn$dfdx <- function (t, y, p, more)
{
nKappa <- more$nKappa
r = array(0, c(length(t), ncol(y), ncol(y)))
dimnames(r) = list(NULL, colnames(y), colnames(y))
pk <- p[(length(p) - (nKappa - 1)):length(p)]
r[, "S", "S"] = - tvtrans(t, pk) * y[,"I"]
r[, "I", "S"] = tvtrans(t, pk) * y[,"I"]
r[, "I", "I"] = -p["gamma"]
return(r)
}
tvDSIRfn$dfdp <- function (t, y, p, more)
{
nKappa <- more$nKappa
r = array(0, c(length(t), ncol(y), length(p)))
dimnames(r) = list(NULL, colnames(y), names(p))
r[, "I", "gamma"] = - y[, "I"]
month <- t %% 1
for(i in 1:nKappa){
r[ , "S", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = - (y[,"I"] * y[, "S"])[month >= (i-1)/nKappa & month < i /nKappa]
r[ , "I", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = (y[,"I"] * y[, "S"])[month >= (i-1)/nKappa & month < i /nKappa]
}
return(r)
}
tvDSIRfn$d2fdx2 <- function (t, y, p, more)
{
r = array(0, c(length(t), ncol(y), ncol(y), ncol(y)))
dimnames(r) = list(NULL, colnames(y), colnames(y), colnames(y))
return(r)
}
tvDSIRfn$d2fdxdp <- function (t, y, p, more)
{
r = array(0, c(length(t), ncol(y), ncol(y), length(p)))
dimnames(r) = list(NULL, colnames(y), colnames(y), names(p))
nKappa <- more$nKappa
r[, "I", "I", "gamma"] = -1
month <- t %% 1
for(i in 1:nKappa){
r[ , "S", "S", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = - y[,"I"][month >= (i-1)/nKappa & month < i /nKappa]
r[ , "I", "S", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = y[,"I"][month >= (i-1)/nKappa & month < i /nKappa]
r[ , "S", "I", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = - y[,"S"][month >= (i-1)/nKappa & month < i /nKappa]
r[ , "I", "I", paste("k", i, sep ="")][month >= (i-1)/nKappa & month < i /nKappa] = y[,"S"][month >= (i-1)/nKappa & month < i /nKappa]
}
return(r)
}
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.