[R-sig-dyn-mod] Qn about FME package
If I try
sensFun(gw_model,parms=pars,sensvar=c("h"))
I get
"Error in func(parms, ...) : argument "rain" is missing, with no default"
Probably because my model function needs more input parameters than the three I'm trying to evaluate..how do I pass only a subset to sensFun ?
If I try model fitting instead, passing the other parameters:
Residuals=function(p)(DataT2$h-gw_model(p,hmin.under.flow, pd.under.flow, rain))
print(system.time(P<-modFit(f=Residuals,p=pars)))
I get: user system elapsed 0.02 0.00 0.02
P$par
hini rf sy
150.000 0.100 0.005
$hessian
hini rf sy
hini 0 0 0
rf 0 0 0
sy 0 0 0
It has not found a new set of parameters.
-----Original Message-----
From: R-sig-dynamic-models [mailto:r-sig-dynamic-models-bounces at r-project.org] On Behalf Of Thomas Petzoldt
Sent: Wednesday, July 20, 2016 12:13 PM
To: Special Interest Group for Dynamic Simulation Models in R
Subject: Re: [R-sig-dyn-mod] Qn about FME package
I guess the error may be in:
Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
See example of ?sensFun, where func is a wrapper around the model that returns the modelled (i.e. dependent or state) variables. In the statement above, it was called with a cost function.
Thomas P.
Am 20.07.2016 um 18:02 schrieb Vishal Mehta:
#########################################
#Jan 2016: Dr. Satkumar Tomer's 1-D model copied ####GW model 1-D GW
model function gw_model <- function(pars, hmin.under.flow,
pd.under.flow, rain){
n <- length(rain)
net_recharge <- rf*rain
h <- array(0,n)
discharge <- array(0,n)
time = array(0,n)
h[1] <- hini
time[1] = 1
for (i in 1:(n-1)){
discharge[i] <- sy*(1-pd.under.flow)*(h[i]-hmin.under.flow+ net_recharge[i]/sy)
if (discharge[i]<0){
discharge[i] <- 0
}
h[i+1] <- h[i] + (net_recharge[i]-discharge[i])/sy
time[i+1] = time[i] + 1
}
return(data.frame(time=time,h=h,discharge=discharge))
}
#synthetic forcing
n = 12*20 # temporal length in months
set.seed(123) # to have the repeatibility in results rain <- rgamma(n,
shape=8, rate = 0.1) # synthetic rainfall
#barplot(rain,xlab="Time",ylab="Rainfall (mm)",main="Monthly
rainfall") rain=rain/1000#mm to m of rain #new
#
#input parameters
rf <- 0.1 # recharge factor, fraction
sy <- 0.005 #specific yield, dimensionless hini <- 150 # initial
condition, metres hmin.under.flow <- 100
pd.under.flow<-0.967
pars=c("hini"=hini,"rf"=rf,"sy"=sy)##this is the subset of parameters
to be SA'd, calibrated, and to inform UA ##
pdf("output.pdf")
output=gw_model(pars,hmin.under.flow,pd.under.flow,rain)
plot(output$h,type='l',xlab='time',ylab='simulated gw level(m)')
dev.off()
##########################
#July 15, 2016
#trial following FME tutorial
#########
library(FME)
#Create dummy observations
#n=240 length of simulated output
#output$h contains the simulated head
#below, observation dataframe is called DataT. It has to(?) contain a
"time" and "h" variable like the simulation output
DataT=data.frame(cbind("time"=seq(1,length(n)),"h"=output$h+rnorm(sd=0
.5,mean=5,n=240),"sd"=4.5))
DataT2=DataT[,c(1,2)]
##
#Create a function to calculate fit
gwcost=function(pars){
out=gw_model(pars,hmin.under.flow, pd.under.flow, rain)
#cost=modCost(model=out,obs=DataT,err="sd")
cost=modCost(model=out,obs=DataT2)
return(cost=cost)
}
###Above seems to be working(maybe not quite?)
gwcost(pars)
gwcost(pars)$model
#[1] 1075 with no weighting
plot(gwcost(pars)$residuals$res)
#
##Sensitivity using sensFun (here onwards results are not following
the JSS article) ##local sensitivity does not seem to be working. does it work only with differential equations?
Sfun=sensFun(gwcost,parms=pars,sensvar=c("h"))
#parameter tinyy to add perturbation? not making any difference
summary(Sfun)
plot(Sfun)
#########
_______________________________________________ R-sig-dynamic-models mailing list R-sig-dynamic-models at r-project.org<mailto:R-sig-dynamic-models at r-project.org> https://stat.ethz.ch/mailman/listinfo/r-sig-dynamic-models