R Graphical Manual

Browse All

Last data update: 2014.03.03

R: Adaptive Multivariate MCMC sampler

mymcmc

R Documentation

Adaptive Multivariate MCMC sampler

Description

This function generates MCMC-based samples from a (posterior) density f (not necessarily normalized). It uses a Metropolis algorithm in conjunction with a multivariate normal proposal distribution which is updated adaptively by monitoring the correlations of succesive increments of at least 2 pilot chains. The method is described in De Gunst, Dewanji and Luebeck (submitted). The adaptive method is similar to the one proposed in Gelfand and Sahu (JCGS 3:261–276, 1994).

Usage

mymcmc(x, nlogf, m1, m2=m1, m3, scl1=0.5, scl2=2, skip=1, covm=0, nfcn=0, plot=FALSE)

Arguments

`x`	a list with components 'label' (of mode character), 'est' (the parameter vector with the initial guess), 'low' (vector with lower bounds), and 'upp' (vector with upper bounds)
`nlogf`	negative log of the density function (not necessarily normalized) for which samples are to be obtained
`m1`	length of first pilot run (not used when covm supplied)
`m2`	length of second pilot run (not used when covm supplied )
`m3`	length of final run
`scl1`	scale for covariance of mv normal proposal (second pilot run)
`scl2`	scale for covariance of mv normal proposal (final run)
`skip`	number of cycles skipped for graphical output
`covm`	covariance matrix for multivariate normal proposal distribution. If supplied, all pilot runs will be skipped and a run of length m3 will be produced. Useful to continue a simulation from a given point with specified covm
`nfcn`	number of function calls
`plot`	logical variable. If TRUE the chain and the negative log density (nlogf) is plotted. The first m1+m2 cycles are shown in green, other cycles in red

Details

standard output reports a summary of the acceptance fraction, the current values of nlogf and the parameters for every (100*skip) th cycle. Plotted chains show values only for every (skip) th cycle.

Value

list with the following components:

`f`	values of nlogf for the samples obtained
`mcmc`	the chain (samples obtained)
`covm`	current covariance matrix for mv normal proposal distribution

Note

This function is part of the Bhat exploration tool

Author(s)

E. Georg Luebeck (FHCRC)

References

too numerous to be listed here

Examples

# generate some Poisson counts on the fly
  dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
  data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))

# neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
  nlogf <- function (x) { 
  ds <- data[, 1]
  y  <- data[, 2]
  g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
  return(sum(g - y * log(g)))
  }

# start MCMC near mle
  x <- list(label=c("a","b","c"), est=c(20, 0.5, 0.05), low=c(0,0,0), upp=c(100,10,.1))
# samples from posterior density (~exp(-nlogf))) with non-informative
# (random uniform) priors for "a", "b" and "c".
out <- mymcmc(x, nlogf, m1=2000, m2=2000, m3=10000, scl1=0.5, scl2=2, skip=10, plot=TRUE)
# start MCMC from some other point: e.g. try x$est <- c(16,.2,.02)

Results


R version 3.3.1 (2016-06-21) -- "Bug in Your Hair"
Copyright (C) 2016 The R Foundation for Statistical Computing
Platform: x86_64-pc-linux-gnu (64-bit)

R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.

R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.

Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.

> library(Bhat)
> png(filename="/home/ddbj/snapshot/RGM3/R_CC/result/Bhat/mcmc.Rd_%03d_medium.png", width=480, height=480)
> ### Name: mymcmc
> ### Title: Adaptive Multivariate MCMC sampler
> ### Aliases: mymcmc
> ### Keywords: iteration methods optimize
> 
> ### ** Examples
> 
> # generate some Poisson counts on the fly
>   dose <- c(rep(0,50),rep(1,50),rep(5,50),rep(10,50))
>   data <- cbind(dose,rpois(200,20*(1+dose*.5*(1-dose*0.05))))
> 
> # neg. log-likelihood of Poisson model with 'linear-quadratic' mean: 
>   nlogf <- function (x) { 
+   ds <- data[, 1]
+   y  <- data[, 2]
+   g <- x[1] * (1 + ds * x[2] * (1 - x[3] * ds)) 
+   return(sum(g - y * log(g)))
+   }
> 
> # start MCMC near mle
>   x <- list(label=c("a","b","c"), est=c(20, 0.5, 0.05), low=c(0,0,0), upp=c(100,10,.1))
> # samples from posterior density (~exp(-nlogf))) with non-informative
> # (random uniform) priors for "a", "b" and "c".
> out <- mymcmc(x, nlogf, m1=2000, m2=2000, m3=10000, scl1=0.5, scl2=2, skip=10, plot=TRUE)
Mon Jul  4 15:16:25 2016 
trying a proposal COVM using the Hessian 
reducing step size 
first pilot chain: using MLEs and log-transformed covariance 
 
n: 1000 acceptance: 0.443 -log lkh: -25987.048 19.6963 0.57153 0.0556884 
n: 2000 acceptance: 0.43 -log lkh: -25989.011 19.5295 0.600104 0.0571235 

 
second pilot run: 
 
eigen values: 0.1098368 0.0001304661 4.361226e-07 
 
n: 3000 acceptance: 0.41 -log lkh: -25992.487 19.8981 0.509162 0.0498643 
n: 4000 acceptance: 0.381 -log lkh: -25991.694 20.12 0.510698 0.0513577 
n: 5000 acceptance: 0.373 -log lkh: -25990.769 20.3621 0.495959 0.0506152 
n: 6000 acceptance: 0.388 -log lkh: -25985.982 21.4851 0.469771 0.0521133 
n: 7000 acceptance: 0.406 -log lkh: -25991.333 20.3786 0.512277 0.0533012 
n: 8000 acceptance: 0.411 -log lkh: -25990.638 20.8496 0.499264 0.0522433 
n: 9000 acceptance: 0.407 -log lkh: -25991.935 20.0247 0.555182 0.0551934 
n: 10000 acceptance: 0.387 -log lkh: -25990.269 20.1706 0.519683 0.0540228 
n: 11000 acceptance: 0.383 -log lkh: -25992.701 20.5678 0.486743 0.0518812 
n: 12000 acceptance: 0.37 -log lkh: -25987.449 21.2937 0.464376 0.0501704 
> # start MCMC from some other point: e.g. try x$est <- c(16,.2,.02)
> 
> 
> 
> 
> 
> dev.off()
null device 
          1 
>