Fitting a dynamic model, and determining the number of parameters that can be fitted.
[This article was first published on gRaphics!, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
Let’s suppose that we have the same dynamic model we presented before – that is, the Lorentz system of differential equations. Remember?Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(deSolve) # require this library | |
solveLorenz <- function(pars, times=seq(0,100,by=0.1)) { | |
derivs <- function(t, state, pars) { # returns rate of change | |
with(as.list(c(state, pars)), { | |
dX <- a*X + Y*Z | |
dY <- b * (Y-Z) | |
dZ <- -X*Y + c*Y - Z | |
return(list(c(dX, dY, dZ))) | |
} | |
) | |
} | |
state <- c(X = 1, Y = 1, Z = 1) | |
## ode solves the model by integration... | |
return(ode(y = state, times = times, func = derivs, parms = pars)) | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
library(FME) | |
Objective <- function(x, parset = names(x)) { | |
guess_pars[parset] <- x | |
tout <- seq(0, 50, by = 0.5) | |
## output times | |
bullet <- solveLorenz(guess_pars, tout) | |
## Model cost | |
return(modCost(obs = target, model = bullet)) | |
} |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
target_pars <- c(a = -9/3, b = -5, c = 30); target_pars | |
target <- solveLorenz(target_pars) | |
guess_pars<-c(a = -6/3, b = -8, c = 20); guess_pars | |
guess <- solveLorenz(guess_pars) | |
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective))) # this works | |
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=100), lower=c(-5,-10,0)))) # this works too | |
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="SANN", control=list(maxit=10000)))) # this works too | |
print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="Nelder-Mead", control=list(maxit=1000)))) # this works too | |
# print(system.time(Fit <- modFit(p = guess_pars, f = Objective, method="Nelder-Mead", control=list(maxit=10000)))) # this works too | |
Fit$par | |
bullet_pars <- Fit$par | |
bullet <- solveLorenz(bullet_pars) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
pXY<-ggplot(as.data.frame(bullet)) +geom_path(aes(X, Y, alpha=Z), col='green') + opts(legend.position = "none") | |
pXY<-pXY +geom_path(data=as.data.frame(target), aes(X, Y, alpha=Z), col='blue') | |
pXY<-pXY +geom_path(data=as.data.frame(guess), aes(X, Y, alpha=Z), col='red') | |
print(pXY) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
sf<-sensFun(func = Objective, parms = bullet_pars, varscale = 1) | |
Coll <- collin(sF) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
> Coll | |
a b c N collinearity | |
1 1 1 0 2 1.1 | |
2 1 0 1 2 2.8 | |
3 0 1 1 2 1.0 | |
4 1 1 1 3 2.9 | |
> |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
require(reshape) | |
# select only the columns relative to parameter on/off, and the collinearity value | |
Coll_m<-melt(t(Coll[,c(1:length(guess_pars),ncol(Coll))])); | |
# split them in two different frames for easier handling | |
Coll_f<-Coll_m[Coll_m$X1=='collinearity',]; | |
Coll_m<-Coll_m[Coll_m$X1!='collinearity',]; | |
# add a categorical label | |
Coll_f$label[Coll_f$value<20]<- sprintf("%3.2g", Coll_f$value[Coll_f$value<20]); Coll_f$label[Coll_f$value>=20] <- ">= 20"; Coll_f$label[Coll_f$value>=100] <- "> 100"; | |
Coll_f$flag<-'ID OK'; Coll_f$flag[Coll_f$value>=20] <- "ID KO"; | |
Coll_f$flag <- factor(Coll_f$flag, levels=c('ID OK', 'ID KO')) |
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
grid.newpage() | |
# create a tile matrix with 0 and 1 depending on which parameter is being fitted. | |
cl1<-ggplot(Coll_m) + geom_tile(aes(X1, X2, fill=value), , col='white') + scale_fill_gradient(low='white', high='black', limits=c(0, 1)) | |
# fill the same matrix with tiles of two different colors | |
cl1<-cl1 + geom_text(aes(X1, X2, label=value, col=value)) + scale_color_gradient(low='black', high='white', limits=c(0, 1)) | |
cl1<-cl1 + opts(title="Parameters' N-ways Identifiability", axis.title.x = theme_blank(), axis.title.y = theme_blank(), axis.text.y = theme_text(colour = NA), axis.ticks = theme_segment(colour = NA), axis.text.x = theme_text(size = 12), legend.position='none') | |
# now, on the right side, fill colors depending on the collinearity values | |
cl2<-ggplot(Coll_f) + geom_tile(aes(X1, X2, fill=flag), col='white') + scale_fill_manual(values=c('green','red')) | |
# finally print clipped collinearity values | |
cl2<-cl2 + geom_text(data=Coll_f, aes(X1, X2, label=label), col='black') + theme_bw() + guides(fill=guide_legend(title=NULL)) + | |
# and all options to make the graph look right | |
opts(title="", axis.title.x = theme_blank(), axis.title.y = theme_blank(), axis.text.y = theme_text(colour = NA), axis.ticks = theme_segment(colour = NA), axis.text.x = theme_text(size = 12), legend.position='none') | |
pushViewport(viewport(layout = grid.layout(6, 7))); | |
print(cl1, vp=vplayout(1:6,1:5)) | |
print(cl2, vp=vplayout(1:6,6:7)) |
This graph suggests, like the table from which is derived, that all parameters can be fit together (with these specific data), and a reliable fit can be found… If only my real-life example were so easy… Homeworks: Try adding some random noise to the XYZ values of the target, and see if the collinearity analysis changes. I’ll try to update this post with a tentative walkthrough in the next few days. For now, enjoy!
To leave a comment for the author, please follow the link and comment on their blog: gRaphics!.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.