diff options
| author | jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb> | 2010-05-11 23:03:37 +0000 |
|---|---|---|
| committer | jranke <jranke@edb9625f-4e0d-4859-8d74-9fd3b1da38cb> | 2010-05-11 23:03:37 +0000 |
| commit | 30cbb4092f6d2d3beff5800603374a0d009ad770 (patch) | |
| tree | ef75421d92823b5b7add1b2d5da9c7499dfd59e5 /R/mkinfit.R | |
Initial upload of the upcoming multicompartment version of kinfit.
Some functionality is still missing (chi2), some may never be implemented
(FOMC model), but in general it is much more powerful than kinfit, owing
to the powerful FME package.
git-svn-id: svn+ssh://svn.r-forge.r-project.org/svnroot/kinfit/pkg/mkin@8 edb9625f-4e0d-4859-8d74-9fd3b1da38cb
Diffstat (limited to 'R/mkinfit.R')
| -rw-r--r-- | R/mkinfit.R | 66 |
1 files changed, 66 insertions, 0 deletions
diff --git a/R/mkinfit.R b/R/mkinfit.R new file mode 100644 index 0000000..9651fd6 --- /dev/null +++ b/R/mkinfit.R @@ -0,0 +1,66 @@ +mkinfit <- function(mkinmod, observed,
+ parms.ini = rep(0.1, length(mkinmod$parms)),
+ state.ini = c(100, rep(0, length(mkinmod$diffs) - 1)),
+ fixed_parms = rep(FALSE, length(mkinmod$parms)),
+ fixed_initials = c(FALSE, rep(TRUE, length(mkinmod$diffs) - 1)),
+ plot = NULL,
+ err = NULL, weight = "none", scaleVar = FALSE,
+ ...)
+{
+ # Name the parameters if they are not named yet
+ if(is.null(names(parms.ini))) names(parms.ini) <- mkinmod$parms
+ # Create a function calculating the differentials specified by the model
+ mkindiff <- function(t, state, parms) {
+ diffs <- vector()
+ for (box in names(mkinmod$diffs))
+ {
+ diffname <- paste("d", box, sep="_")
+ diffs[diffname] <- with(as.list(c(state, parms)),
+ eval(parse(text=mkinmod$diffs[[box]])))
+ }
+ return(list(c(diffs)))
+ }
+
+ # Name the inital parameter values if they are not named yet
+ if(is.null(names(state.ini))) names(state.ini) <- names(mkinmod$diffs)
+
+ # TODO: Collect parameters to be optimised
+ parms.optim <- parms.ini[!fixed_parms]
+ parms.fixed <- parms.ini[fixed_parms]
+
+ state.ini.optim <- state.ini[!fixed_initials]
+ state.ini.optim.boxnames <- names(state.ini.optim)
+ names(state.ini.optim) <- paste(names(state.ini.optim), "0", sep="_")
+ state.ini.fixed <- state.ini[fixed_initials]
+
+ # Define the model cost function
+ cost <- function(P)
+ {
+ if(length(state.ini.optim) > 0) {
+ odeini <- c(P[1:length(state.ini.optim)], state.ini.fixed)
+ names(odeini) <- c(state.ini.optim.boxnames, names(state.ini.fixed))
+ } else odeini <- state.ini.fixed
+
+ odeparms <- c(P[(length(state.ini.optim) + 1):length(P)], parms.fixed)
+ # Solve the ODE
+ out <- ode(
+ y = odeini,
+ times = unique(observed$time),
+ func = mkindiff,
+ parms = odeparms)
+
+ # Output transformation for models with ghost compartments like SFORB
+ out_transformed <- data.frame(time = out[,"time"])
+ for (var in names(mkinmod$map)) {
+ if(length(mkinmod$map[[var]]) == 1) {
+ out_transformed[var] <- out[, var]
+ } else {
+ out_transformed[var] <- rowSums(out[, mkinmod$map[[var]]])
+ }
+ }
+
+ return(modCost(out_transformed, observed, y = "value",
+ err = err, weight = weight, scaleVar = scaleVar))
+ }
+ modFit(cost, c(state.ini.optim, parms.optim), ...)
+}
|