diff options
Diffstat (limited to 'CakeIrlsFit.R')
| -rw-r--r-- | CakeIrlsFit.R | 146 |
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) } |