[R-sig-dyn-mod] Qn about FME package

[Vishal Mehta]

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

	[[alternative HTML version deleted]]