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 Interface Reference

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

Public Member Functions

subroutine mcmc_dp (eval, likelihood, para, rangepar, mcmc_paras, burnin_paras, seed_in, printflag_in, maskpara_in, restart, restart_file, 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 either known or 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.

This sampler does not change the best parameter set, i.e. it cannot be used as an optimiser.
However, the serial and the parallel version give therefore the bitwise same results.
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: Baysian 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)Interface Function which calculates likelihood of given parameter set x
    [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
    Likelihood has to be defined as a function interface
    The maximal number of parameters is 1000.
    Authors
    Juliane Mai
    Date
    Nov 2014

Definition at line 203 of file mo_mcmc.F90.

Member Function/Subroutine Documentation

◆ mcmc_dp()

subroutine mo_mcmc::mcmc::mcmc_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,
logical, intent(in), optional  restart,
character(len = *), intent(in), optional  restart_file,
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 433 of file mo_mcmc.F90.

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