View source: R/derivedEquations.R
sensitivitiesSymb | R Documentation |
Compute sensitivity equations of a function symbolically
sensitivitiesSymb( f, states = names(f), parameters = NULL, inputs = NULL, events = NULL, reduce = FALSE )
f |
named vector of type character, the functions |
states |
Character vector. Sensitivities are computed with respect to initial values of these states |
parameters |
Character vector. Sensitivities are computed with respect to initial values of these parameters |
inputs |
Character vector. Input functions or forcings. They are excluded from the computation of sensitivities. |
events |
data.frame of events with columns "var" (character, the name of the state to be
affected), "time" (numeric or character, time point),
"value" (numeric or character, value), "method" (character, either
"replace" or "add"). See events.
Within |
reduce |
Logical. Attempts to determine vanishing sensitivities, removes their equations and replaces their right-hand side occurences by 0. |
The sensitivity equations are ODEs that are derived from the original ODE f. They describe the sensitivity of the solution curve with respect to parameters like initial values and other parameters contained in f. These equtions are also useful for parameter estimation by the maximum-likelihood method. For consistency with the time-continuous setting provided by adjointSymb, the returned equations contain attributes for the chisquare functional and its gradient.
Named vector of type character with the sensitivity equations. Furthermore,
attributes "chi" (the integrand of the chisquare functional), "grad" (the integrand
of the gradient of the chisquare functional), "forcings" (Character vector of the
additional forcings being necessare to compute chi
and grad
) and "yini" (
The initial values of the sensitivity equations) are returned.
## Not run: ###################################################################### ## Sensitivity analysis of ozone formation ###################################################################### library(deSolve) # O2 + O <-> O3 f <- c( O3 = " build_O3 * O2 * O - decay_O3 * O3", O2 = "-build_O3 * O2 * O + decay_O3 * O3", O = "-build_O3 * O2 * O + decay_O3 * O3" ) # Compute sensitivity equations f_s <- sensitivitiesSymb(f) # Generate ODE function func <- funC(c(f, f_s)) # Initialize times, states, parameters and forcings times <- seq(0, 15, by = .1) yini <- c(O3 = 0, O2 = 3, O = 2, attr(f_s, "yini")) pars <- c(build_O3 = .1, decay_O3 = .01) # Solve ODE out <- odeC(y = yini, times = times, func = func, parms = pars) # Plot solution par(mfcol=c(2,3)) t <- out[,1] M1 <- out[,2:4] M2 <- out[,5:7] M3 <- out[,8:10] M4 <- out[,11:13] M5 <- out[,14:16] M6 <- out[,17:19] matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="solution") legend("topright", legend = c("O3", "O2", "O"), lty=1, col=1:3) matplot(t, M2, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O3)") matplot(t, M3, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O2)") matplot(t, M4, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d O)") matplot(t, M5, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d build_O3)") matplot(t, M6, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="d/(d decay_O3)") ## End(Not run) ## Not run: ###################################################################### ## Estimate parameter values from experimental data ###################################################################### library(deSolve) # O2 + O <-> O3 # diff = O2 - O3 # build_O3 = const. f <- c( O3 = " build_O3 * O2 * O - decay_O3 * O3", O2 = "-build_O3 * O2 * O + decay_O3 * O3", O = "-build_O3 * O2 * O + decay_O3 * O3" ) # Compute sensitivity equations and get attributes f_s <- sensitivitiesSymb(f) chi <- attr(f_s, "chi") grad <- attr(f_s, "grad") forcings <- attr(f_s, "forcings") # Generate ODE function func <- funC(f = c(f, f_s, chi, grad), forcings = forcings, fcontrol = "nospline", modelname = "example3") # Initialize times, states, parameters times <- seq(0, 15, by = .1) yini <- c(O3 = 0, O2 = 2, O = 2.5) yini_s <- attr(f_s, "yini") yini_chi <- c(chi = 0) yini_grad <- rep(0, length(grad)); names(yini_grad) <- names(grad) pars <- c(build_O3 = .2, decay_O3 = .1) # Initialize forcings (the data) data(oxygenData) forcData <- data.frame(time = oxygenData[,1], name = rep( colnames(oxygenData[,-1]), each=dim(oxygenData)[1]), value = as.vector(oxygenData[,-1])) forc <- setForcings(func, forcData) # Solve ODE out <- odeC(y = c(yini, yini_s, yini_chi, yini_grad), times = times, func = func, parms = pars, forcings = forc, method = "lsodes") # Plot solution par(mfcol=c(1,2)) t <- out[,1] M1 <- out[,2:4] M2 <- out[,names(grad)] tD <- oxygenData[,1] M1D <- oxygenData[,2:4] matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states") matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE) legend("topright", legend = names(f), lty=1, col=1:3) matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient") legend("topleft", legend = names(grad), lty=1, col=1:5) # Define objective function obj <- function(p) { out <- odeC(y = c(p[names(f)], yini_s, yini_chi, yini_grad), times = times, func = func, parms = p[names(pars)], forcings = forc, method="lsodes") value <- as.vector(tail(out, 1)[,"chi"]) gradient <- as.vector( tail(out, 1)[,paste("chi", names(p), sep=".")]) hessian <- gradient%*%t(gradient) return(list(value = value, gradient = gradient, hessian = hessian)) } # Fit the data myfit <- optim(par = c(yini, pars), fn = function(p) obj(p)$value, gr = function(p) obj(p)$gradient, method = "L-BFGS-B", lower=0, upper=5) # Model prediction for fit parameters prediction <- odeC(y = c(myfit$par[1:3], yini_s, yini_chi, yini_grad), times = times, func = func, parms = myfit$par[4:5], forcings = forc, method = "lsodes") # Plot solution par(mfcol=c(1,2)) t <- prediction[,1] M1 <- prediction[,2:4] M2 <- prediction[,names(grad)] tD <- oxygenData[,1] M1D <- oxygenData[,2:4] matplot(t, M1, type="l", lty=1, col=1:3, xlab="time", ylab="value", main="states") matplot(tD, M1D, type="b", lty=2, col=1:3, pch=4, add=TRUE) legend("topright", legend = names(f), lty=1, col=1:3) matplot(t, M2, type="l", lty=1, col=1:5, xlab="time", ylab="value", main="gradient") legend("topleft", legend = names(grad), lty=1, col=1:5) ## End(Not run)
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.