summaryrefslogtreecommitdiff
path: root/CakeIrlsFit.R
diff options
context:
space:
mode:
authorJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:18:07 +0200
committerJohannes Ranke <jranke@uni-bremen.de>2017-10-18 11:18:07 +0200
commit3d6b4b4b8293a4a4ab6f06805e1380600373796c (patch)
tree4ec4dd041427b0154684122f5ae2101c7385f5e2 /CakeIrlsFit.R
parentbe6d42ef636e8e1c9fdcfa6f8738ee10e885d75b (diff)
Version 2.0v2.0
Diffstat (limited to 'CakeIrlsFit.R')
-rw-r--r--CakeIrlsFit.R146
1 files changed, 136 insertions, 10 deletions
diff --git a/CakeIrlsFit.R b/CakeIrlsFit.R
index bf20572..05eff0a 100644
--- a/CakeIrlsFit.R
+++ b/CakeIrlsFit.R
@@ -1,9 +1,9 @@
-#$Id: CakeIrlsFit.R 216 2011-07-05 14:35:03Z nelr $
+#$Id$
#
# The CAKE R modules are based on mkin
# Modifications developed by Tessella Plc for Syngenta: Copyright (C) 2011 Syngenta
-# Authors: Rob Nelson, Richard Smith
-# Tessella Project Reference: 6245
+# Authors: Rob Nelson, Richard Smith, Tamar Christina
+# Tessella Project Reference: 6245, 7247
# This program is free software: you can redistribute it and/or modify
# it under the terms of the GNU General Public License as published by
@@ -22,7 +22,8 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)), lower = 0,
upper = Inf, fixed_parms = NULL, fixed_initials = names(mkinmod$diffs)[-1],
plot = FALSE, quiet = FALSE, err = NULL, weight = "none",
- scaleVar = FALSE, atol=1e-6, control=list(),...)
+ scaleVar = FALSE, atol=1e-6, sannMaxIter = 10000, control=list(),
+ useSolnp = FALSE, method='L-BFGS-B',...)
{
### This is a modification based on the "mkinfit" function.
@@ -35,21 +36,25 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
observed <- subset(observed, name %in% names(mkinmod$map))
ERR <- rep(1,nrow(observed))
observed <- cbind(observed,err=ERR)
+ flag <- 0
+
obs_vars = unique(as.character(observed$name))
if (is.null(names(parms.ini)))
names(parms.ini) <- mkinmod$parms
+
mkindiff <- function(t, state, parms) {
time <- t
diffs <- vector()
for (box in mod_vars) {
diffname <- paste("d", box, sep = "_")
diffs[diffname] <- with(as.list(c(time, state, parms)),
- eval(parse(text = mkinmod$diffs[[box]])))
+ eval(parse(text = mkinmod$diffs[[box]])))
}
return(list(c(diffs)))
}
if (is.null(names(state.ini)))
names(state.ini) <- mod_vars
+
parms.fixed <- parms.ini[fixed_parms]
optim_parms <- setdiff(names(parms.ini), fixed_parms)
parms.optim <- parms.ini[optim_parms]
@@ -57,6 +62,7 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
optim_initials <- setdiff(names(state.ini), fixed_initials)
state.ini.optim <- state.ini[optim_initials]
state.ini.optim.boxnames <- names(state.ini.optim)
+
if (length(state.ini.optim) > 0) {
names(state.ini.optim) <- paste(names(state.ini.optim),
"0", sep = "_")
@@ -67,8 +73,37 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
############### Iteratively Reweighted Least Squares#############
## Start with no weighting
- fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower,
- upper = upper,control=control,...)
+
+ ## a prefitting step since this is usually the most effective method
+ if(useSolnp)
+ {
+ pnames=names(c(state.ini.optim, parms.optim))
+ fn <- function(P){
+ names(P) <- pnames
+ FF<<-costFunctions$cost(P)
+ return(FF$model)}
+ a <- try(fit <- solnp(c(state.ini.optim, parms.optim),fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+ #optimmethod <- method0
+
+ flag <- 1
+ if(class(a) == "try-error")
+ {
+ print('solnp fails, try PORT or other algorithm by users choice, might take longer time. Do something else!')
+ warning('solnp fails, switch to PORT or other algorithm by users choice')
+ ## now using submethod already
+ a <- try(fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method='Port',control=control))
+ flag <- 0
+ if(class(a) == "try-error")
+ {
+ fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower, upper = upper, method=method,control=control)
+ }
+ }
+ }else{
+ fit <- modFit(costFunctions$cost, c(state.ini.optim, parms.optim), lower = lower,
+ upper = upper,control=control,...)
+ }
+
+ ## print(fit$hessian)
if(length(control)==0)
{
@@ -81,27 +116,109 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
control <- list(irls.control=irls.control)
}
}
+
irls.control <- control$irls.control
maxIter <- irls.control$maxIter
tol <- irls.control$tol
####
+
+ if(length(mkinmod$map)==1 || useSolnp){
+ ## there is only one parent just do one iteration:
+ maxIter <- 0
+
+ if(flag==1)## fit from solnp
+ {
+ fit$ssr <- fit$values[length(fit$values)]
+ fit$residuals <-FF$residual$res
+ ## mean square per varaible
+ if (class(FF) == "modCost") {
+ names(fit$residuals) <- FF$residuals$name
+ fit$var_ms <- FF$var$SSR/FF$var$N
+ fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
+ fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+ names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+ names(fit$var_ms) <- FF$var$name
+ } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+ }
+ err1 <- sqrt(fit$var_ms_unweighted)
+ ERR <- err1[as.character(observed$name)]
+ observed$err <-ERR
+ }
+
niter <- 1
## insure one IRLS iteration
diffsigma <- 100
olderr <- rep(1,length(mod_vars))
while(diffsigma>tol & niter<=maxIter)
{
+ if(flag==1 && useSolnp)## fit from solnp
+ {
+ fit$ssr <- fit$values[length(fit$values)]
+ fit$residuals <-FF$residual$res
+ ## mean square per varaible
+ if (class(FF) == "modCost") {
+ names(fit$residuals) <- FF$residuals$name
+ fit$var_ms <- FF$var$SSR/FF$var$N
+ fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
+ fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+ names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+ names(fit$var_ms) <- FF$var$name
+ } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+ }
+
+
err <- sqrt(fit$var_ms_unweighted)
ERR <- err[as.character(observed$name)]
costFunctions$set.error(ERR)
diffsigma <- sum((err-olderr)^2)
cat("IRLS iteration at",niter, "; Diff in error variance ", diffsigma,"\n")
olderr <- err
- fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+
+ if(useSolnp)
+ {
+ flag <- 1
+ a <- try(fit <- solnp(fit$par,fun=fn,LB=lower,UB=upper,control=control),silent=TRUE)
+
+ if(class(a) == "try-error")
+ {
+ flag <- 0
+ print('solnp fails during IRLS iteration, try PORT or other algorithm by users choice.This may takes a while. Do something else!') ## NOTE: because in kingui we switch off the warnings, we need to print out the message instead.
+ warning('solnp fails during IRLS iteration, switch to PORT or other algorithm by users choice')
+
+ fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, method='Port',control=list())
+ }
+ }else{
+ fit <- modFit(costFunctions$cost, fit$par, lower = lower, upper = upper, control=control, ...)
+ }
+
+ ## print(fit$hessian)
+
niter <- niter+1
### If not converged, reweight and fit
}
+
+ if(flag==1 && useSolnp){
+ ## solnp used
+ optimmethod <- 'solnp'
+ fit$ssr <- fit$values[length(fit$values)]
+ fit$residuals <-FF$residual$res
+ ## mean square per varaible
+ if (class(FF) == "modCost") {
+ names(fit$residuals) <- FF$residuals$name
+ fit$var_ms <- FF$var$SSR/FF$var$N
+ fit$var_ms_unscaled <- FF$var$SSR.unscaled/FF$var$N
+ fit$var_ms_unweighted <- FF$var$SSR.unweighted/FF$var$N
+
+ names(fit$var_ms_unweighted) <- names(fit$var_ms_unscaled) <-
+ names(fit$var_ms) <- FF$var$name
+ } else fit$var_ms <- fit$var_ms_unweighted <- fit$var_ms_unscaled <- NA
+ np <- length(c(state.ini.optim, parms.optim))
+ fit$rank <- np
+ fit$df.residual <- length(fit$residuals) - fit$rank
+ }
###########################################
fit$mkindiff <- mkindiff
@@ -125,15 +242,19 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
fit$fixed$type = c(rep("state", length(state.ini.fixed)),
rep("deparm", length(parms.fixed)))
- fit$errmin <- CakeChi2(observed, predicted_long, obs_vars, parms.optim, state.ini.optim)
+ fit$errmin <- CakeChi2(mkinmod, observed, predicted_long, obs_vars, parms.optim, state.ini.optim, state.ini, parms.ini)
parms.all = c(fit$par, parms.fixed)
fit$penalties <- CakePenaltiesLong(parms.all, out_predicted, observed)
fit$distimes <- data.frame(DT50 = rep(NA, length(obs_vars)),
DT90 = rep(NA, length(obs_vars)), row.names = obs_vars)
+ fit$extraDT50<- data.frame(DT50 = rep(NA, 2), row.names = c("k1", "k2"))
+
for (obs_var in obs_vars) {
type = names(mkinmod$map[[obs_var]])[1]
- fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all)
+ fit$distimes[obs_var, ] = CakeDT(type,obs_var,parms.all,sannMaxIter)
}
+
+ fit$extraDT50[ ,c("DT50")] = CakeExtraDT(type,parms.all)
data <- merge(observed, predicted_long, by = c("time", "name"))
@@ -142,6 +263,11 @@ CakeIrlsFit <- function (mkinmod, observed, parms.ini = rep(0.1, length(mkinmod$
data$residual <- data$observed - data$predicted
data$variable <- ordered(data$variable, levels = obs_vars)
fit$data <- data[order(data$variable, data$time), ]
+
+ fit$costFunctions <- costFunctions
+ fit$upper <- upper
+ fit$lower <- lower
+
class(fit) <- c("CakeFit", "mkinfit", "modFit")
return(fit)
}

Contact - Imprint