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!
No comments:
Post a Comment