0.6.2-dev0
FORCES
FORtran lib for Comp. Env. Sys.
Loading...
Searching...
No Matches
Public Member Functions | List of all members
mo_mcmc::mcmc_stddev Interface Reference

This module is Monte Carlo Markov Chain sampling of a posterior parameter distribution (measurement errors are neither known nor modeled). More...

Public Member Functions

subroutine mcmc_stddev_dp (eval, likelihood, para, rangepar, mcmc_paras, burnin_paras, seed_in, printflag_in, maskpara_in, tmp_file, loglike_in, paraselectmode_in, iter_burnin_in, iter_mcmc_in, chains_in, stepsize_in)
 

Detailed Description

This module is Monte Carlo Markov Chain sampling of a posterior parameter distribution (measurement errors are neither known nor modeled).

Sample posterior parameter distribution with Metropolis Algorithm.
This sampling is performed in two steps, i.e. the burn-in phase for adjusting model dependent parameters for the second step which is the proper sampling using the Metropolis Hastings Algorithm.

1. Burn in phase: Find the optimal step size

Purpose:

Find optimal stepsize for each parameter such that the acceptance ratio converges to a value around 0.3.
Important variables:

Variable Description
burnin_iter length of markov chain performed to calculate acceptance ratio
acceptance ratio ratio between accepted jumps and all trials (LEN)
acceptance multiplier stepsize of a parameter is multiplied with this value when jump is accepted (initial : 1.01)
rejection multiplier stepsize of a parameter is multiplied with this value when jump is rejected (initial : 0.99 and
will never be changed)
stepsize a new parameter value is chosen based on a uniform distribution
pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) (initial : stepsize_i = 1.0 for all i)

Algorithm:

  1. start a new markov chain of length burnin_iter with initial parameter set is the OPTIMAL one
    • select a set of parameters to change:
      • accurate --> choose one parameter,
      • comput. efficient --> choose all parameters,
      • moderate accurate & efficient --> choose half of the parameters
    • change parameter(s) based on their stepsize
    • decide whether changed parameter set is accepted or rejected:
      • \( \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \)
      • random number r = Uniform[0,1]
      • odds Ratio > 0 --> positive accept
        odds Ratio > r --> negative accept
        odds Ratio < r --> reject
    • adapt stepsize of parameters changed:
      • accepted step: stepsize_i = stepsize_i * accptance multiplier
      • rejected step: stepsize_i = stepsize_i * rejection multiplier
    • if step is accepted: for all changed parameter(s) change stepsize
  2. calculate acceptance ratio of the Markov Chain
  3. adjust acceptance multiplier acc_mult and store good ratios in history list
    • acceptance ratio < 0.23 --> acc_mult = acc_mult * 0.99
      delete history list
    • acceptance ratio > 0.44 --> acc_mult = acc_mult * 1.01
      delete history list
    • 0.23 < acceptance ratio < 0.44
      add acceptance ratio to history list
  4. check if already 10 values are stored in history list and if they have converged to a value above 0.3
    ( mean above 0.3 and variance less \(\sqrt{1/12*0.05^2}\) = Variance of uniform [acc_ratio +/- 2.5%] )
    • if check is positive abort and save stepsizes
      else goto (1)

2. Monte Carlo Markov Chain: Sample posterioir distribution of parameter

Purpose:

use the previous adapted stepsizes and perform ONE monte carlo markov chain
the accepted parameter sets show the posterior distribution of parameters
Important variables:

Variable Description
iter_mcmc length of the markov chain (>> iter_burnin)
stepsize a new parameter value is chosen based on a uniform distribution
pnew_i = pold_i + Unif(-stepsize_i, stepsize_i) use stepsizes of the burn-in (1)

Algorithm:

  1. select a set of parameters to change
    • accurate --> choose one parameter,
    • comput. efficient --> choose all parameters,
    • moderate accurate & efficient --> choose half of the parameters
  2. change parameter(s) based on their stepsize
  3. decide whether changed parameter set is accepted or rejected:
    • \( \mathrm{odds Ratio} = \frac{\mathrm{likelihood}(p_{new})}{\mathrm{likelihood}(p_{old})} \)
    • random number r = Uniform[0,1]
    • odds Ratio > 0 --> positive accept
      odds Ratio > r --> negative accept
      odds Ratio < r --> reject
  4. if step is accepted: save parameter set
  5. goto (1)

Example

call mcmc( likelihood, para, rangePar, mcmc_paras, burnin_paras, &
seed_in=seeds, printflag_in=printflag, maskpara_in=maskpara, &
tmp_file=tmp_file, loglike_in=loglike, &
paraselectmode_in=paraselectmode, &
iter_burnin_in=iter_burnin, iter_mcmc_in=iter_mcmc, &
chains_in=chains, stepsize_in=stepsize)

See also example in test directory.

Literature

  1. Gelman et. al (1995) Bayesian Data Analyis. Chapman & Hall.
  2. Gelman et. al (1997) Efficient Metropolis Jumping Rules: Bayesian Statistics 5. pp. 599-607
  3. Tautenhahn et. al (2012) Beyond distance-invariant survival in inverse recruitment modeling: A case study in Siberian Pinus sylvestris forests. Ecological Modelling, 233, 90-103. doi:10.1016/j.ecolmodel.2012.03.009.
    Parameters
    [in]real(dp) :: likelihood(x,sigma,stddev_new,likeli_new)Interface Function which calculates likelihood of given parameter set x and given standard deviation sigma and returns optionally the standard deviation stddev_new of the errors using x and likelihood likeli_new using stddev_new
    [in]real(dp) :: para(:)Inital parameter set (should be GOOD approximation of best parameter set)
    [in]real(dp) :: rangePar(size(para),2)Min/max range of parameters
    [in]integer(i8), optional :: seed_inUser seed to initialise the random number generator
    (default: none --> initialized with timeseed)
    [in]logical, optional :: printflag_inPrint of output on command line
    (default: .False.)
    [in]logical, optional :: maskpara_in(size(para))Parameter will be sampled (.True.) or not (.False.)
    (default: .True.)
    [in]character(len=*), optional :: tmp_filefilename for temporal data saving:
    every iter_mcmc_in iterations parameter sets are appended to this file
    the number of the chain will be prepended to filename
    output format: netcdf
    (default: no file writing)
    [in]logical, optional :: loglike_intrue if loglikelihood function is given instead of likelihood function
    (default: .false.)
    [in]integer(i4), optional :: ParaSelectMode_inHow many parameters will be changed at once?
    • half of the parameter --> 1_i4
    • only one parameter --> 2_i4
    • all parameter --> 3_i4 (default: 2_i4)
    [in]integer(i4), optional :: iter_burnin_inLength of Markov chains of initial burn-in part
    (default: Max(250, 200*count(maskpara)) )
    [in]integer(i4), optional :: iter_mcmc_inLength of Markov chains of proper MCMC part
    (default: 1000 * count(maskpara) )
    [in]integer(i4), optional :: chains_innumber of parallel mcmc chains
    (default: 5_i4)
    [in]real(dp), DIMENSION(size(para,1)), optional :: stepsize_instepsize for each parameter
    if given burn-in is discarded
    (default: none --> adjusted in burn-in)
    [out]real(dp), allocatable :: mcmc_paras(:,:)Parameter sets sampled in proper MCMC part of algorithm
    [out]real(dp), allocatable :: burnin_paras(:,:)Parameter sets sampled during burn-in part of algorithm
    Note
    Restrictions: Likelihood has to be defined as a function interface
    Author
    Maren Goehler
    Date
    Sep. 2012
    • Created using copy of Simulated Annealing:
      • constant temperature T
      • burn-in for stepsize adaption
      • acceptance/ rejection multiplier
    Author
    Juliane Mai
    Date
    Sep. 2012
    • Cleaning code and introduce likelihood:
      • likelihood instead of objective function
      • odds ratio
      • convergence criteria for burn-in
      • different modes of parameter selection
    • OpenMP for chains of MCMC
    • optional file for temporal output
    Nov 2012
    • Temporary file writing as NetCDF
    Aug 2013
    • New likelihood interface to reduce number of function evaluations. Only one seed has to be given.
    Nov 2014
    • add new routine for mcmc:
      • mcmc_stddev, i.e. mcmc with likelihood with "faked sigma".
      • mcmc, i.e. mcmc with "real" likelihood (sigma is given by e.g. error model).
    Apr 2015
    • handling of array-like variables in restart-namelists
    Author
    Mathias Cuntz and Juliane Mai
    Date
    Feb 2015
    • restart

Definition at line 419 of file mo_mcmc.F90.

Member Function/Subroutine Documentation

◆ mcmc_stddev_dp()

subroutine mo_mcmc::mcmc_stddev::mcmc_stddev_dp ( procedure(eval_interface), intent(in), pointer  eval,
procedure(objective_interface), intent(in), pointer  likelihood,
real(dp), dimension(:), intent(in)  para,
real(dp), dimension(:, :), intent(in)  rangepar,
real(dp), dimension(:, :), intent(out), allocatable  mcmc_paras,
real(dp), dimension(:, :), intent(out), allocatable  burnin_paras,
integer(i8), intent(in), optional  seed_in,
logical, intent(in), optional  printflag_in,
logical, dimension(size(para, 1)), intent(in), optional  maskpara_in,
character(len = *), intent(in), optional  tmp_file,
logical, intent(in), optional  loglike_in,
integer(i4), intent(in), optional  paraselectmode_in,
integer(i4), intent(in), optional  iter_burnin_in,
integer(i4), intent(in), optional  iter_mcmc_in,
integer(i4), intent(in), optional  chains_in,
real(dp), dimension(size(para, 1)), intent(in), optional  stepsize_in 
)

Definition at line 1289 of file mo_mcmc.F90.

The documentation for this interface was generated from the following file: