当前位置：文档库 › FMEsteady

FMEsteady

R Package FME:Inverse Modelling,Sensitivity, Monte Carlo–Applied to a Steady-State Model

Karline Soetaert

NIOZ Yerseke

The Netherlands

Abstract

R package FME(Soetaert and Petzoldt2010)contains functions for model calibration, sensitivity,identi?ability,and Monte Carlo analysis of nonlinear models.

This vignette,(vignette("FMEsteady")),applies FME to a partial di?erential equa-tion,solved with a steady-state solver from package rootSolve

A similar vignette(vignette("FMEdyna")),applies the functions to a dynamic simi-

lation model,solved with integration routines from package deSolve

A third vignette(vignette("FMEother")),applies the functions to a simple nonlinear

model

vignette("FMEmcmc")tests the Markov chain Monte Carlo(MCMC)implementation Keywords:steady-state models,di?erential equations,?tting,sensitivity,Monte Carlo,iden-ti?ability,R.

1.A steady-state model of oxygen in a marine sediment This is a simple model of oxygen in a marine(submersed)sediment,di?using along a spatial gradient,with imposed upper boundary concentration oxygen is consumed at maximal?xed rate,and including a monod limitation.

See(Soetaert and Herman2009)for a description of reaction-transport models.

The constitutive equations are:

?O2?t =?

?F lux

?cons·

O2+k s

F lux=?D·

?O2

O2(x=0)=upO2

>par(mfrow=c(2,2))

>require(FME)

First the model parameters are de?ned...

>pars<-c(upO2=360,#concentration at upper boundary,mmolO2/m3 +cons=80,#consumption rate,mmolO2/m3/day

2FME–Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

+ks=1,#O2half-saturation ct,mmolO2/m3

+D=1)#diffusion coefficient,cm2/d

Next the sediment is vertically subdivided into100grid cells,each0.05cm thick.

>n<-100#nr grid points

>dx<-0.05#cm

>dX<-c(dx/2,rep(dx,n-1),dx/2)#dispersion distances;half dx near boundaries >X<-seq(dx/2,len=n,by=dx)#distance from upper interface at middle of box The model function takes as input the parameter values and returns the steady-state condition

of oxygen.Function steady.1D from package rootSolve((Soetaert2009))does this in a very

e?cient way(see(Soetaert and Herman2009)).

>O2fun<-function(pars)

+derivs<-function(t,O2,pars)

+with(as.list(pars),{

+Flux<--D*diff(c(upO2,O2,O2[n]))/dX

+dO2<--diff(Flux)/dx-cons*O2/(O2+ks)

+return(list(dO2,UpFlux=Flux[1],LowFlux=Flux[n+1]))

+})

+#Solve the steady-state conditions of the model

+ox<-steady.1D(y=runif(n),func=derivs,parms=pars,

+nspec=1,positive=TRUE)

+data.frame(X=X,O2=ox$y)

The model is run

>ox<-O2fun(pars)

and the results plotted...

>plot(ox$O2,ox$X,ylim=rev(range(X)),xlab="mmol/m3",

+main="Oxygen",ylab="depth,cm",type="l",lwd=2)

2.Global sensitivity analysis:Sensitivity ranges

The sensitivity of the oxygen pro?le to parameter cons,the consumption rate is estimated.

We assume a normally distributed parameter,with mean=80(parMean),and a variance=100 (parCovar).The model is run100times(num).

Karline Soetaert 3

50100150200250300350

4321

0Oxygen

mmol/m3

d e p t h , c m

Figure 1:The modeled oxygen pro?le -see text for R -code

>print(system.time(

+Sens2<-sensRange(parms =pars,func =O2fun,dist ="norm",+num =100,parMean =c(cons =80),parCovar =100)+))

user system elapsed

1.7

0.0 1.7

The results can be plotted in two ways:>

par(mfrow =c(1,2))

>plot(Sens2,xyswap =TRUE,xlab ="O2",

+ylab ="depth,cm",main ="Sensitivity runs")>plot(summary(Sens2),xyswap =TRUE,xlab ="O2",

+ylab ="depth,cm",main ="Sensitivity ranges")>

par(mfrow =c(1,1))

3.Local sensitivity analysis :Sensitivity functions

Local sensitivity analsysis starts by calculating the sensitivity functions >O2sens <-sensFun(func=O2fun,parms=pars)

The summary of these functions gives information about which parameters have the largest e?ect (univariate sensitivity):>summary(O2sens)

4FME –Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

50100150

200250300350

4321

0Sensitivity runs

d e p t h , c

Figure 2:Results of the sensitivity run -left:all model runs,right:summary -see text for R -code

value scale L1L2Mean Min Max N

upO23603607.00.887.0 1.0e+0013.4176100cons 80808.11.14-8.1-2.2e+01-0.0084100ks 112.20.37 2.2 1.2e-049.6137100D 118.11.148.18.4e-0322.0312

100

In bivariate sensitivity the pair-wise relationship and the correlation is estimated and/or plotted:

>pairs(O2sens)>cor(O2sens[,-(1:2)])

upO2cons ks D

upO2 1.0000000-0.97879450.83758060.9787945cons -0.9787945 1.0000000-0.9317287-1.0000000ks 0.8375806-0.9317287 1.00000000.9317287D 0.9787945-1.00000000.9317287 1.0000000

Multivariate sensitivity is done by estimating the collinearity between parameter sets (Brun,Reichert,and Kunsch 2001).>Coll <-collin(O2sens)>Coll

upO2cons ks D N collinearity

110027.72

10102 2.9

Karline Soetaert 5

upO2

?20?10?50

●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

●●●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●05101520

2468

12●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

●●●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●?20?10?5

?0.98

cons

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

0.84?0.93

02468●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●

●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●●2468120510

15200.98?102468

0.93

Figure 3:pairs plot -see text for R -code

100127.7401102 4.4

50101267108864.0

600112 4.471110325.5811013NaN 91011325.5

100111361616597.4

114

NaN

>plot(Coll,log ="y")

4.Fitting the model to the data

Assume both the oxygen ?ux at the upper interface and a vertical pro?le of oxygen has been measured.

These are the data:

>O2dat <-data.frame(x =seq(0.1,3.5,by =0.1),+y =c(279,260,256,220,200,203,189,179,165,140,138,127,116,+109,92,87,78,72,62,55,49,43,35,32,27,20,15,15,10,8,5,3,2,1,0))>O2depth <-cbind(name ="O2",O2dat)#oxygen versus depth >O2flux <-c(UpFlux =170)#measured flux First a function is de?ned that returns only the required model output.>O2fun2<-function(pars)+{

6FME –Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

1e +01

1e +031e +051e +07

Collinearity

Number of parameters

C o l l i n e a r i t y i n d e

Figure 4:collinearity -see text for R -code

+derivs<-function(t,O2,pars)+{+with (as.list(pars),{++Flux <--D*diff(c(upO2,O2,O2[n]))/dX +dO2<--diff(Flux)/dx -cons*O2/(O2+ks)++return(list(dO2,UpFlux =Flux[1],LowFlux =Flux[n+1]))+})+}+

+ox <-steady.1D(y =runif(n),func =derivs,parms =pars,nspec =1,+positive =TRUE,rtol =1e-8,atol =1e-10)+

+list(data.frame(x =X,O2=ox$y),+UpFlux =ox$UpFlux)+}

The function used in the ?tting algorithm returns an instance of type modCost .This is created by calling function modCost twice.First with the modeled oxygen pro?le,then with the modeled ?ux.

>Objective <-function (P)+{

+Pars <-pars

+Pars[names(P)]<-P

+modO2<-O2fun2(Pars)+

+#Model cost:first the oxygen profile

+Cost <-modCost(obs =O2depth,model =modO2[[1]],

Karline Soetaert7 +x="x",y="y")

+#then the flux

+modFl<-c(UpFlux=modO2$UpFlux)

+Cost<-modCost(obs=O2flux,model=modFl,x=NULL,cost=Cost)

+return(Cost)

We?rst estimate the identi?ability of the parameters,given the data:

>print(system.time(

+sF<-sensFun(Objective,parms=pars)

+))

user system elapsed

0.140.000.14

>summary(sF)

value scale L1L2Mean Min Max N

upO23603604.30.97 4.30.506913.336

cons80803.70.99-3.6-15.37220.536

ks110.40.140.4-0.0069 3.136

D113.70.99 3.70.034215.436

>collin(sF)

upO2cons ks D N collinearity

1110028.6

210102 3.1

3100128.7

401102 4.2

50101250.6

600112 4.2

71110314.2

81101350.8

91011314.7

100111350.6

111111451.0

The collinearity of the full set is too high,but as the oxygen di?usion coe?cient is well known,it is left out of the?tting.The combination of the three remaining parameters has a low enough collinearity to enable automatic?tting.The parameters are constrained to be>0 >collin(sF,parset=c("upO2","cons","ks"))

8FME–Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

upO2cons ks D N collinearity

11110314

>print(system.time(

+Fit<-modFit(p=c(upO2=360,cons=80,ks=1),

+f=Objective,lower=c(0,0,0))

+))

user system elapsed

1.250.00 1.26

>(SFit<-summary(Fit))

Parameters:

Estimate Std.Error t value Pr(>|t|)

upO2292.937 2.104139.245<2e-16***

cons49.687 2.36621.001<2e-16***

ks 1.297 1.3620.9530.348

---

Signif.codes:0′S***ˇS0.001′S**ˇS0.01′S*ˇS0.05′S.ˇS0.1′SˇS1

Residual standard error:4.401on33degrees of freedom

Parameter correlation:

upO2cons ks

upO21.00000.57910.2975

cons0.57911.00000.9011

ks0.29750.90111.0000

We next plot the residuals

>plot(Objective(Fit$par),xlab="depth",ylab="",

+main="residual",legpos="top")

and show the best-?t model

>Pars<-pars

>Pars[names(Fit$par)]<-Fit$par

>modO2<-O2fun(Pars)

>plot(O2depth$y,O2depth$x,ylim=rev(range(O2depth$x)),pch=18,

+main="Oxygen-fitted",xlab="mmol/m3",ylab="depth,cm")

>lines(modO2$O2,modO2$X)

5.Running a Markov chain Monte Carlo

We use the parameter covariances of previous?t to update parameters,while the mean squared residual of the?t is use as prior fo the model variance.

Karline Soetaert 9

0.0

0.5 1.0 1.5 2.0 2.5 3.0 3.5

?10

?50

510

residual

depth

Figure 5:residuals -see text for R -code

3.5

3.02.52.01.51.00.5

0.0Oxygen?fitted

d e p t h , c m

Figure 6:Best ?t model -see text for R -code

10FME–Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model >Covar<-SFit$cov.scaled*2.4^2/3

>s2prior<-SFit$modVariance

We run an adaptive Metropolis,making sure that ks does not become negative...

>print(system.time(

+MCMC<-modMCMC(f=Objective,p=Fit$par,jump=Covar,

+niter=1000,ntrydr=2,var0=s2prior,wvar0=1,

+updatecov=100,lower=c(NA,NA,0))

+))

number of accepted runs:699out of1000(69.9%)

user system elapsed

27.780.0228.34

>MCMC$count

dr_steps Alfasteps num_accepted num_covupdate

68420526999

Plotting the results is similar to previous cases.

>plot(MCMC,Full=TRUE)

>hist(MCMC,Full=TRUE)

>pairs(MCMC,Full=TRUE)

or summaries can be created:

>summary(MCMC)

upO2cons ks var_model

mean294.66980153.960056 4.265416618.373713e+02

sd 6.058268 6.329367 4.335104161.127736e+04

min275.59673142.3785600.023612731.982599e+00

max321.49078393.96143037.846720562.610703e+05

q025291.71970449.888962 1.538161031.245861e+01

q050293.60550452.477486 3.159669523.436431e+01

q075296.24619156.415191 5.813362601.034256e+02

>cor(MCMC$pars)

upO2cons ks

upO21.00000000.70644040.3040173

cons0.70644041.00000000.8218140

ks0.30401730.82181401.0000000

Karline Soetaert 11

02006001000280

290

300

310320

upO2

iter

020********

7080

cons

iter

0200

6001000

102030

iter

0200

6001000

10002000300040005000

SSR

iter

0200

6001000

1e +01

1e +03

1e +05

var_model

iter

v a r i a n c

Figure 7:MCMC plot results -see text for R -code

12FME –Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

upO2

2803003200.00

0.05

0.10

0.1

cons ?

406080

0.00

0.04

0.08

0.1

102030

0.00

0.050.10

0.1

SSR

1000

3000

5000

0.0000.0010.0020.0030.00

error var_model

v a r _m o d e l

100000250000

0e +00

2e ?044e ?0

Figure 8:MCMC histogram results -see text for R -code

Karline Soetaert 13

507090100030005000

280300320

0102030

1000

3000

5000

28030032001020300100000250000

0100000250000

Figure 9:MCMC pairs plot -see text for R -code

14FME –Inverse Modelling,Sensitivity,Monte Carlo with a Steady-State Model

050100150200250300

5432

Figure 10:MCMC range plot -see text for R -code

Note:we pass to sensRange the full parameter vector (parms )and the parameters sampled during the MCMC (parInput ).

>plot(summary(sensRange(parms =pars,parInput =MCMC$par,f =O2fun,num =500)),+xyswap =TRUE)

>points(O2depth$y,O2depth$x)

6.Finally

This vignette is made with Sweave (Leisch 2002).

References

Brun R,Reichert P,Kunsch H (2001).“Practical Identi?ability Analysis of Large Environ-mental Simulation Models.”Water Resources Research ,37(4),1015–1030.

Karline Soetaert15 Leisch F(2002).“Dynamic Generation of Statistical Reports Using Literate Data Analysis.”In W H¨a rdle,B R¨o nz(eds.),“COMPSTAT2002–Proceedings in Computational Statistics,”pp.575–580.Physica-Verlag,Heidelberg.

Soetaert K(2009).rootSolve:Nonlinear Root Finding,Equilibrium and Steady-State Analysis of Ordinary Di?erential Equations.R package version1.6,URL http://CRAN. https://www.wendangku.net/doc/5e6540206.html,/package=rootSolve.

Soetaert K,Herman PMJ(2009).A Practical Guide to Ecological https://www.wendangku.net/doc/5e6540206.html,ing R as a Simulation Platform.Springer-Verlag,New York.

Soetaert K,Petzoldt T(2010).“Inverse Modelling,Sensitivity and Monte Carlo Analysis in R Using Package FME.”Journal of Statistical Software,33(3),1–28.URL http://www. https://www.wendangku.net/doc/5e6540206.html,/v33/i03/.

A?liation:

Karline Soetaert

Royal Netherlands Institute of Sea Research(NIOZ)

4401NT Yerseke,Netherlands

E-mail:k.soetaert@nioz.nl

URL:http://www.nioz.nl