Install package from Github

Only needs to be done the first time :

library(devtools)
install_github("PecanProject/pecan", subdir="all")

Parameter Data Assimilation in PEcAn

Currently, there are four ways of doing Parameter Data Assimilation (PDA) in PEcAn :

Which one to use?

bruteforce : You can choose bruteforce as a method to use PEcAn’s natively implemented Metropolis-Hastings Markov Chain Monte Carlo (MCMC) algorithm. This would perform a Bayesian MCMC on model parameters by proposing one parameter value at a time, and accepting or rejecting it according to the calculated likelihood value. This algorithm also has an adaptation functionality that can be turned on (recommended) and off (see below). As each (i-th) parameter is dependent of the previous (i-1-th) one, this algorithm can only be run sequentially (but different chains can be run in parallel). Therefore, if you have just one parameter to calibrate and a relatively fast model that can be run couple of hundreds (couple of thousands would even be better) of times within an hour, it is possible to use this algorithm.

bruteforce.bs : This algorithm is basically identical to the bruteforce, but rather than proposing parameters one at a time, it proposes new values for all parameters at once (“bs” stands for “block sampling”). If you have more than one parameter to calibrate and a relatively fast model, you can use this algorithm, preferably with the adaptation functionality turned on.

emulator : When a model is slow, it is practically not possible to run it many times in order to explore the parameter space and draw enough samples to converge the target distribution with bruteforce algorithms. Instead, we can run the model for a relatively smaller number of times with parameter values that have been carefully chosen to give good coverage of parameter space. Then we can interpolate the likelihood calculated for each of those runs to get a surface that “emulates” the true likelihood and perform regular MCMC (just like the “bruteforce” approach), except instead of actually running the model on every iteration to get a likelihood, this time we will just get an approximation from the likelihood emulator.

bayesian.tools : There are other MCMC algorithms with different proposal schemes and acceptance criterion than Metropolis-Hastings. BayesianTools is an R-package that includes MCMC and SMC samplers and other tools for Bayesian parameter calibration. If you choose bayesian.tools option, PEcAn framework will hand the PDA calculations over to the BayesianTools package. Although this package includes algorithms that are designed to explore the parameter space more efficiently than the regular MH-MCMC, you would still need a relatively faster model to use these algorithms. The BayesianTools R-package itself is currently under development, once it is fully integrated with PEcAn it will be possible to run some of those algorithms in parallel.

Adding PDA tags to pecan.xml

The easiest way to use PEcAn’s parameter data assimilation module is to add an <assim.batch> block to pecan.xml, load the file with read.settings, and pass the resulting settings object to pda.mcmc(). There are some differences in the settings for using different PDA methods (see below), but here is an example <assim.batch> block :

<assim.batch>
  <iter>10000</iter>
  <method>bruteforce</method>
  <chain>3</chain>
  <param.names>
   <temperate.coniferous>                        <--- YOUR PFT
    <param>veg_respiration_Q10</param>
    <param>stem_respiration_rate</param>
    <param>Vm_low_temp</param>
   </temperate.coniferous>
  </param.names>
  <inputs>
    <file>
      <input.id>1000000384</input.id>
      <path>
        <path>/fs/data1/pecan.data/input/Ameriflux_site_0-796/US-Bar.2005.nc</path>
        <path>/fs/data1/pecan.data/input/Ameriflux_site_0-796/US-Bar.2006.nc</path>
      </path>
      <likelihood>Laplace</likelihood>
      <variable.id>298</variable.id>
      <variable.name>
        <variable.name>LE</variable.name>
        <variable.name>UST</variable.name>
      </variable.name>
    </file>
  </inputs>
  <jump>
    <ar.target>0.3</ar.target>
    <adapt>200</adapt>
    <adj.min>0.1</adj.min>
  </jump>
  <diag.plot.iter>500</diag.plot.iter>
</assim.batch>

Here are details about the settings:

data(trait.dictionary, package = "PEcAn.utils")

head(trait.dictionary[,c("id", "figid")])
##                         id                   figid
## 1           plant_min_temp          Plant Min Temp
## 2                      GDD     Growing Degree Days
## 3                    leafC            Leaf C conc.
## 4                 c2n_leaf                Leaf C:N
## 5 leaf_respiration_rate_m2   Leaf Respiration Rate
## 6  dark_respiration_factor Dark Respiration Factor

Method specific settings

If you want to extend your previous bruteforce and bruteforce.bs MCMC-chains, you can add an <extension> tag within your <assim.batch> chunk after reading your post-PDA settings (which would be saved as “pecan.pda[UNIQUE_ID].xml” in your working directory) and set it to “longer”.

  ...
  <method>bruteforce.bs</method>
  <iter>10000</iter>
  <method>bruteforce</method>
  <chain>3</chain>
  <extension>longer</extension>
  ...

If you are using methods other than bruteforce and bruteforce.bs, some additional tags and different settings may apply.

emulator would use additional tags:

  ...
  <method>emulator</method>
  <n.knot>20</n.knot>
  <GPpckg>GPfit</GPpckg>
  <chain>3</chain>
  <extension>round</extension>
  <knot.par>0.75</knot.par>
  ...

bayesian.tools would look for sampler specific settings that can be passed through the pecan.xml as a block under the <bt.settings> tag. Currently, the available samplers in the BayesianTools package are:

The name of the chosen sampler would be passed under <sampler> tag within the <bt.settings> block :

It is not possible to use some of the samplers in this package for univariate cases. The ones that you can use for univariate cases are: “Metropolis”, “DE”, “DEzs” and “Twalk”.

  ...
  <method>bayesian.tools</method>
  <bt.settings>
    <iter>10000</iter>
    <sampler>Metropolis</sampler>
    <DRlevels>1</DRlevels>
    <optimize>FALSE</optimize>
    <adapt>TRUE</adapt>
    <adaptionInterval>200</adaptionInterval>
    <adaptationNotBefore>500</adaptationNotBefore>
    <consoleUpdates>100</consoleUpdates>
  </bt.settings>
  ...

Some of the samplers are also restartable : “Metropolis”, “DE”, “DEzs”, “DREAM”, “DREAMzs”, “Twalk”