MCMCBioGro: Simulated Annealing and Markov chain Monte Carlo for...
In BioCro: Simulation of Biomass Crops

Description Usage Arguments Details Value Note Author(s) See Also Examples

This function atempts to implement the simulated annealing method for estimating parameters of a generic C4 crop growth.

MCMCBioGro(niter = 10, niter2 = 10, phen = 6, iCoef = NULL, saTemp = 5, 
    coolSamp = 20, scale = 0.5, WetDat, data, day1 = NULL, dayn = NULL, 
    timestep = 1, lat = 40, iRhizome = 7, irtl = 1e-04, canopyControl = list(), 
    seneControl = list(), photoControl = list(), phenoControl = list(), 
    soilControl = list(), nitroControl = list(), centuryControl = list(), 
    sd = c(0.02, 1e-06))

`niter`	number of iterations for the simulated annealing portion of the optimization.
`niter2`	number of iterations for the Markov chain Monte Carlo portion of the optimization.
`phen`	Phenological stage being optimized.
`iCoef`	initial coefficients for dry biomass partitioning.
`saTemp`	simulated annealing temperature.
`coolSamp`	number of cooling samples.
`scale`	scale parameter to control the standard deviations.
`WetDat`	weather data.
`data`	observed data.
`day1`	first day of the growing season.
`dayn`	last day of the growing season.
`timestep`	Timestep see `BioGro`.
`lat`	latitude.
`iRhizome`	initial rhizome biomass.
`irtl`	See `BioGro`.
`canopyControl`	See `canopyParms`.
`seneControl`	See `seneParms`.
`photoControl`	See `photoParms`.
`phenoControl`	See `phenoParms`.
`soilControl`	See `soilParms`.
`nitroControl`	See `nitroParms`.
`centuryControl`	See `centuryParms`.
`sd`	standard deviations for the parameters to be optimized. The first (0.02) is for the positive dry biomass partitioning coefficients. The second (1e-6) is for the negative dry biomass partitioning coefficients.

This function implements a hybrid algorithm where the first portion is simulated annealing and the second portion is a Markov chain Monte Carlo. The user controls the number of iterations in each portion of the chain with niter and niter2.

An object of class MCMCBioGro consisting of a list with 23 components. The easiest way of accessing the information is with the print and plot methods.

The automatic method for guessing the last day of the growing season differs slightly from that in BioGro. To prevent error due to a shorter simulated growing season than the observed one the method in MCMCBioGro adds one day to the last day of the growing season. Although the upper limit is fixed at 330.

Fernando E. Miguez

See Also as BioGro, OpBioGro and constrOpBioGro.

## Not run: 

data(weather05)

## Some coefficients
pheno.ll <- phenoParms(kLeaf1=0.48,kStem1=0.47,kRoot1=0.05,kRhizome1=-1e-4,
                       kLeaf2=0.14,kStem2=0.65,kRoot2=0.21, kRhizome2=-1e-4,
                       kLeaf3=0.01, kStem3=0.56, kRoot3=0.13, kRhizome3=0.3, 
                       kLeaf4=0.01, kStem4=0.56, kRoot4=0.13, kRhizome4=0.3,
                       kLeaf5=0.01, kStem5=0.56, kRoot5=0.13, kRhizome5=0.3,
                       kLeaf6=0.01, kStem6=0.56, kRoot6=0.13, kRhizome6=0.3)
                                 
system.time(ans <- BioGro(weather05, phenoControl = pheno.ll))

ans.dat <- as.data.frame(unclass(ans)[1:11])
sel.rows <- seq(1,nrow(ans.dat),400)
simDat <- ans.dat[sel.rows,c("ThermalT","Stem","Leaf","Root","Rhizome","Grain","LAI")]
plot(ans,simDat)

## Residual sum of squares before the optimization

ans0 <- BioGro(weather05)
RssBioGro(simDat,ans0)

## This is not working now
## op1.mc <- MCMCBioGro(phen=1, niter=200,niter2=200,
                     WetDat=weather05,
                     data=simDat)


## plot(op1.mc)

## plot(op1.mc, plot.kind="trace", subset = nams 
##				c("kLeaf_1","kStem_1","kRoot_1"))


## End(Not run)