Ada>> fmcmc>> 返回
项目作者： USCbiostats

项目描述 ：
A friendly MCMC framework
高级语言： R
项目地址: git://github.com/USCbiostats/fmcmc.git
创建时间： 2017-04-18T18:53:29Z
项目社区：https://github.com/USCbiostats/fmcmc
开源协议：Other
下载


" class="reference-link">fmcmc: A friendly MCMC framework

DOI
R CI
Build status
Coverage Status
Lifecycle: stable
CRAN status
CRAN downloads
status
Integrative Methods of Analysis for Genetic Epidemiology

What

The fmcmc R package provides a lightweight general framework for
implementing Markov Chain Monte Carlo methods based on the
Metropolis-Hastings algorithm. This implementation’s primary purpose
lies in the fact that the user can incorporate the following flexibly:

  1. Automatic convergence checker: The algorithm splits the MCMC
    runs according to the frequency with which it needs to check
    convergence. Users can use either one of the included functions
    (convergence_gelman, convergence_geweke, etc.), or provide their
    own.

  2. Run multiple chains in parallel fashion: Using either a PSOCK
    cluster (default), or providing a personalized cluster object like
    the ones in the parallel R package.

  3. User defined transition kernels: Besides of canonical Gaussian
    Kernel, users can specify their own or use one of the included in
    the package, for example: kernel_adapt, kernel_ram,
    kernel_normal_reflective, kernel_unif, kernel_mirror, or
    kernel_unif_reflective.

All the above without requiring compiled code. For the latest about
fmcmc, checkout the NEWS.md section.

Who is this for?

While a lot of users rely on MCMC tools such as Stan (via the
rstan package) or WinBUGS
(via rstan), in several
settings either these tools are not enough or provide too much for
things that do not need that much. So, this tool is for you if:

  • You have a simple model to estimate with Metropolis-Hastings.

  • You want to run multiple chains of your model using out-of-the-box
    parallel computing.

  • You don’t want (or cannot) rely on external tools (so you need
    good-old base R only for your models).

  • You want to implement a model in which your model parameters are
    either bounded (like a standard error, for example), or are not,
    say, continuous (e.g., a size variable in a Binomial distribution),
    so you need a personalized transition kernel.

In any other case, you may want to take a look at the previously
mentioned R packages, or check out the
mcmc R package, which also
implements the Metropolis-Hastings algorithm (although with not all the
features that this R package has), the
adaptMCMC R package, or
the MCMCpack R package.

Installing

If you want to get the latest bleeding-edge version from Github, you can
use devtools:

  1. devtools::install_github("USCbiostats/fmcmc")

The latest (stable) release is also available on CRAN:

  1. install.packages("fmcmc")

Citation

  1. To cite fmcmc in publications use:
  3.   Vega Yon et al., (2019). fmcmc: A friendly MCMC framework. Journal of
  4.   Open Source Software, 4(39), 1427,
  5.   https://doi.org/10.21105/joss.01427
  7. And the actual R package:
  9. Vega Yon G, Marjoram P (2021). _fmcmc: A friendly MCMC framework_. doi:
  10. 10.5281/zenodo.3378987 (URL: https://doi.org/10.5281/zenodo.3378987), R
  11. package version 0.5-0, <URL: https://github.com/USCbiostats/fmcmc>.
  13. To see these entries in BibTeX format, use 'print(<citation>,
  14. bibtex=TRUE)', 'toBibtex(.)', or set
  15. 'options(citation.bibtex.max=999)'.

Example: Linear regression model

First run

In the following we show how to use the package for estimating
parameters in a linear regression model. First, let’s simulate some data
to use:

  1. set.seed(78845)
  2. n <- 1000
  3. X <- rnorm(n)
  4. y <- 3.0 + 2.0*X + rnorm(n, sd = 4)

In this case, we have three parameters to estimate, the constant (2.0),
the β coefficient (2.0), and the standard deviation parameter of the
error (1.5).

To estimate this model, we can either maximize the log-likelihood
function–which is what is usually done–or we could do it using MCMC. In
either case, we need to specify the log(unnormalized some
times)-likelihood function:

  1. ll <- function(p, X., y.) {
  3.   joint_ll <- dnorm(y. - (p[1] + X.*p[2]), sd = p[3], log = TRUE)
  4.   joint_ll <- sum(joint_ll)
  6.   # If is undefined, then we explicitly return infinte (instead of NaN, for
  7.   # example)
  8.   if (!is.finite(joint_ll))
  9.     return(-Inf)
  11.   joint_ll
  12. }

Notice that the function has more than one argument, in this case, p,
the vector of parameters, X. and y., the data of our model.

Let’s do a first run of the MCMC algorithm using the function of the
same name (first, load the package, of course):

  1. library(fmcmc)
  3. # Running the MCMC (we set the seed first)
  4. set.seed(1215)
  5. ans <- MCMC(
  6.   ll,
  7.   initial = c(0, 0, sd(y)),
  8.   nsteps  = 5000,
  9.   X.      = X,
  10.   y.      = y
  11.   )

As the output object is an object of class mcmc from the coda R
package, we can use all the functions from it on our output:

  1. library(coda)
  2. plot(ans)

  1. summary(ans)
  1. ## 
  2. ## Iterations = 1:5000
  3. ## Thinning interval = 1 
  4. ## Number of chains = 1 
  5. ## Sample size per chain = 5000 
  6. ## 
  7. ## 1. Empirical mean and standard deviation for each variable,
  8. ##    plus standard error of the mean:
  9. ## 
  10. ##       Mean      SD Naive SE Time-series SE
  11. ## par1 3.113 0.17593 0.002488       0.024341
  12. ## par2 1.975 0.10647 0.001506       0.022105
  13. ## par3 4.093 0.07843 0.001109       0.005951
  14. ## 
  15. ## 2. Quantiles for each variable:
  16. ## 
  17. ##       2.5%   25%   50%   75% 97.5%
  18. ## par1 2.975 3.029 3.068 3.255 3.354
  19. ## par2 1.749 1.907 1.980 2.020 2.145
  20. ## par3 3.978 4.070 4.101 4.102 4.226

While the summary statistics look very good (we got very close to the
original parameters), the trace of the parameters looks very bad (poor
mixing). We can re-run the algorithm changing the scale parameter in the
kernel_normal function. To do so, we can pass ans as the initial
argument so that the function starts from the last point of that chain:

  1. ans <- MCMC(
  2.   ll,
  3.   initial = ans,
  4.   nsteps  = 5000,
  5.   X.      = X,
  6.   y.      = y,
  7.   kernel  = kernel_normal(scale = .05) # We can set the scale parameter like this
  8.   )
  9. plot(ans)

Much better! Now, what if we use Vihola (2012) Robust Adaptive
Metropolis (which is also implemented in the R package
adaptMCMC)

  1. ans_RAM <- MCMC(
  2.   ll,
  3.   initial = ans,
  4.   nsteps  = 5000,
  5.   X.      = X,
  6.   y.      = y,
  7.   kernel  = kernel_ram() 
  8.   )
  9. plot(ans_RAM)

  1. 1 - rejectionRate(ans_RAM)
  1. ##      par1      par2      par3 
  2. ## 0.3522705 0.3522705 0.3522705

We can also try using Haario et al. (2001) Adaptive Metropolis

  1. ans_AM <- MCMC(
  2.   ll,
  3.   initial = ans,
  4.   nsteps  = 5000,
  5.   X.      = X,
  6.   y.      = y,
  7.   kernel  = kernel_adapt() 
  8.   )
  9. plot(ans_AM)

  1. 1 - rejectionRate(ans_AM)
  1. ##      par1      par2      par3 
  2. ## 0.5457091 0.5457091 0.5457091

Finally, if needed, we can also access information about the last run
using MCMC_OUTPUT. For example, if we wanted to look at the trace of
the logposterior function, we could use the get_logpost() function:

  1. plot(get_logpost(), type = "l")

The set of proposed values is also available using the get_draws()
function:

  1. boxplot(get_draws(), type = "l")

If the previous run featured multiple chains, then get_logpost() would
return a list instead of length get_nchains().

Automatic stop

Now, suppose that the algorithm actually takes a lot of time to actually
reach stationary state, then it would be nice to actually sample from
the posterior distribution once convergence has been reached. In the
following example we use multiple chains and the Gelman-Rubin diagnostic
to check for convergence of the chain:

  1. set.seed(1215) # Same seed as before
  2. ans <- MCMC(
  3.   ll,
  4.   initial = c(0, 0, sd(y)),
  5.   nsteps  = 5000,
  6.   X.      = X,
  7.   y.      = y,
  8.   kernel  = kernel_normal(scale = .05),
  9.   nchains = 2,                           # Multiple chains
  10.   conv_checker = convergence_gelman(200) # Checking for conv. every 200 steps
  11.   )
  1. ## No convergence yet (steps count: 200). Gelman-Rubin's R: 4.5843. Trying with the next bulk.
  3. ## No convergence yet (steps count: 400). Gelman-Rubin's R: 1.1877. Trying with the next bulk.
  5. ## No convergence yet (steps count: 600). Gelman-Rubin's R: 1.4297. Trying with the next bulk.
  7. ## No convergence yet (steps count: 800). Gelman-Rubin's R: 1.1582. Trying with the next bulk.
  9. ## No convergence yet (steps count: 1000). Gelman-Rubin's R: 1.3414. Trying with the next bulk.
  11. ## No convergence yet (steps count: 1200). Gelman-Rubin's R: 1.2727. Trying with the next bulk.
  13. ## No convergence yet (steps count: 1400). Gelman-Rubin's R: 1.4456. Trying with the next bulk.
  15. ## No convergence yet (steps count: 1600). Gelman-Rubin's R: 1.3792. Trying with the next bulk.
  17. ## No convergence yet (steps count: 1800). Gelman-Rubin's R: 1.2069. Trying with the next bulk.
  19. ## No convergence yet (steps count: 2000). Gelman-Rubin's R: 1.1789. Trying with the next bulk.
  21. ## No convergence yet (steps count: 2200). Gelman-Rubin's R: 1.1208. Trying with the next bulk.
  23. ## No convergence yet (steps count: 2400). Gelman-Rubin's R: 1.1196. Trying with the next bulk.
  25. ## Convergence has been reached with 2600 steps. Gelman-Rubin's R: 1.0792. (2600 final count of samples).

As a difference from the previous case, now we didn’t have to wait until
the 5,000 completed, but the algorithm stopped for us, allowing us to
start generating the desired sample much quicker.

Kernels: Making sure that we get positive values

For this final example, we will use the kernel argument and provide
what corresponds to a transition kernel which makes proposals within
certain boundaries, in particular, we want the algorithm to propose only
positive values for the sd parameter, which we now must be positive.

Moreover, since we know that we will only get positive values, we can go
further and modify ll skipping the check for finite values:

  1. ll <- function(p, X., y.) {
  3.   sum(dnorm(y. - (p[1] + X.*p[2]), sd = p[3], log = TRUE))
  5. }

Much simpler function! Let’s do the call of the MCMC function specifying
the right transition kernel to increase the acceptance rate. In this
example, we will set the max of all parameters to be 5.0, and the min to
be -5.0 for the constant and 0 for the beta coefficient and the variance
parameter, all this using the kernel_normal_reflective (which
implements a normal kernel with boundaries) function:

  1. set.seed(1215) # Same seed as before
  2. ans <- MCMC(
  3.   ll,
  4.   initial = c(0, 0, sd(y)),
  5.   nsteps  = 5000,
  6.   X.      = X,
  7.   y.      = y,
  8.   kernel  = kernel_normal_reflective(
  9.     ub    = 5.0,               # All parameters have the same upper bound
  10.     lb    = c(-5.0, 0.0, 0.0), # But lower bound is specified per parameter
  11.     scale = 0.05               # This is the same scale as before
  12.     ),
  13.   nchains = 2,                           
  14.   conv_checker = convergence_gelman(200)
  15.   )
  1. ## No convergence yet (steps count: 200). Gelman-Rubin's R: 3.7891. Trying with the next bulk.
  3. ## No convergence yet (steps count: 400). Gelman-Rubin's R: 1.1257. Trying with the next bulk.
  5. ## No convergence yet (steps count: 600). Gelman-Rubin's R: 1.4696. Trying with the next bulk.
  7. ## No convergence yet (steps count: 800). Gelman-Rubin's R: 1.1313. Trying with the next bulk.
  9. ## No convergence yet (steps count: 1000). Gelman-Rubin's R: 1.4384. Trying with the next bulk.
  11. ## No convergence yet (steps count: 1200). Gelman-Rubin's R: 1.3696. Trying with the next bulk.
  13. ## No convergence yet (steps count: 1400). Gelman-Rubin's R: 1.5243. Trying with the next bulk.
  15. ## No convergence yet (steps count: 1600). Gelman-Rubin's R: 1.3720. Trying with the next bulk.
  17. ## No convergence yet (steps count: 1800). Gelman-Rubin's R: 1.1722. Trying with the next bulk.
  19. ## No convergence yet (steps count: 2000). Gelman-Rubin's R: 1.1492. Trying with the next bulk.
  21. ## No convergence yet (steps count: 2200). Gelman-Rubin's R: 1.1004. Trying with the next bulk.
  23. ## No convergence yet (steps count: 2400). Gelman-Rubin's R: 1.1161. Trying with the next bulk.
  25. ## Convergence has been reached with 2600 steps. Gelman-Rubin's R: 1.0815. (2600 final count of samples).

Again, as the proposal kernel has lower and upper bounds, then we are
guaranteed that all proposed states are within the support of the
parameter space.

Other tools

fmcmc is just one way to work with Markov Chain Monte Carlo. Besides
stan and WinBUGS, there are other ways to do MCMC in R:
mcmc,
HybridMC,
adaptMCMC, and
elhmc among others (take a
look at the CRAN Task View on Bayesian
Inference.)

Contributing to fmcmc

We welcome contributions to fmcmc. Whether reporting a bug, starting a
discussion by asking a question, or proposing/requesting a new feature,
please go by creating a new issue
here so that we can talk
about it.

Please note that the ‘fmcmc’ project is released with a Contributor
Code of Conduct. By contributing to this project,
you agree to abide by its terms.

Funding

Supported by National Cancer Institute Grant #1P01CA196596.