77 lines
2.4 KiB
R
77 lines
2.4 KiB
R
library(RBBGCMuso)
|
|
file.copy("../../start_set/maize.epc","./start_set/maize.epc",overwrite=TRUE)
|
|
setwd("start_set/")
|
|
|
|
rmse <- function(modelled, measured){
|
|
sqrt(mean((modelled-measured)**2))
|
|
}
|
|
|
|
r2 <- function(modelled, measured){
|
|
summary(lm("mod ~ meas", data= data.frame(mod=modelled,meas=measured)))$r.squared
|
|
}
|
|
|
|
modeff <- function(modelled, measured){
|
|
1 - (sum((modelled-measured)**2) / sum((measured - mean(measured))**2))
|
|
}
|
|
|
|
bias <- function(modelled, measured){
|
|
mean(modelled) - mean(measured)
|
|
}
|
|
|
|
get_stats <- function(modelled,measured){
|
|
c(r2=r2(modelled,measured),
|
|
rmse=rmse(modelled,measured),
|
|
bias=bias(modelled,measured),
|
|
modeff=modeff(modelled,measured))
|
|
}
|
|
|
|
get_modelled <- function(obsTable,settings, ...){
|
|
simulation <- runMuso(settings, ...)
|
|
yield <- simulation[,"fruit_DM"]
|
|
modelled <- yield[match(as.Date(obsTable$date),as.Date(names(yield),"%d.%m.%Y"))] * 10
|
|
modelled
|
|
}
|
|
|
|
obsTable <- read.csv2("Martonvasar_maize.obs",stringsAsFactors=FALSE)
|
|
obsTable$mean <- obsTable$mean / 1000 * 0.85
|
|
obsTable$sd <- obsTable$sd / 1000
|
|
measured <- obsTable$mean
|
|
|
|
# apriori
|
|
settings <- setupMuso(iniInput=c("n.ini","n.ini"))
|
|
modelled <- get_modelled(obsTable, settings)
|
|
results <- matrix(ncol=4,nrow=11)
|
|
colnames(results) <- c("r2","rmse","bias","modeff")
|
|
row.names(results) <- 0:10
|
|
results[1,] <- get_stats(modelled,measured)
|
|
# max_likelihood_stats
|
|
for(i in 1:10){
|
|
file.copy(sprintf("../%s/maize_ml.epc",i),"maize.epc",overwrite=TRUE)
|
|
settings <- setupMuso(iniInput=c("n.ini","n.ini"))
|
|
modelled <- get_modelled(obsTable, settings)
|
|
results[i+1,] <- get_stats(modelled, measured)
|
|
}
|
|
# median_stats
|
|
results_med <- results
|
|
|
|
for(i in 1:10){
|
|
file.copy(sprintf("../maize_median_step%02d_corrected.epc",i),"maize.epc",overwrite=TRUE)
|
|
settings <- setupMuso(iniInput=c("n.ini","n.ini"))
|
|
modelled <- get_modelled(obsTable, settings)
|
|
results_med[i+1,] <- get_stats(modelled, measured)
|
|
}
|
|
results_med
|
|
|
|
succes_ratio <- numeric(10)
|
|
for(i in 1:10){
|
|
succes_ratio[i] <- sum(read.csv(sprintf("../%s/result.csv",i),stringsAsFactors=FALSE)$Const[-1])/10000
|
|
}
|
|
names(succes_ratio) <- 1:10
|
|
png("../success_rate.png",height=30,width=30,res=300, units="cm")
|
|
barplot(succes_ratio,ylim=c(0,1),xlab="iteration number",ylab="Succes rate")
|
|
dev.off()
|
|
|
|
postscript("../success_rate.eps")
|
|
barplot(succes_ratio,ylim=c(0,1),xlab="iteration number",ylab="Succes rate")
|
|
dev.off()
|