MCMC | R Documentation |

A flexible implementation of the Metropolis-Hastings MCMC algorithm, users can utilize arbitrary transition kernels as well as set-up an automatic stop criterion using a convergence check test.

```
MCMC(
initial,
fun,
nsteps,
...,
seed = NULL,
nchains = 1L,
burnin = 0L,
thin = 1L,
kernel = kernel_normal(),
multicore = FALSE,
conv_checker = NULL,
cl = NULL,
progress = interactive() && !multicore,
chain_id = 1L
)
MCMC_without_conv_checker(
initial,
fun,
nsteps,
...,
nchains = 1L,
burnin = 0L,
thin = 1L,
kernel = kernel_normal(),
multicore = FALSE,
conv_checker = NULL,
cl = NULL,
progress = interactive() && !multicore,
chain_id = 1L
)
MCMC_OUTPUT
```

`initial` |
Either a numeric matrix or vector, or an object of class coda::mcmc or coda::mcmc.list (see details). initial values of the parameters for each chain (See details). |

`fun` |
A function. Returns the log-likelihood. |

`nsteps` |
Integer scalar. Length of each chain. |

`...` |
Further arguments passed to |

`seed` |
If not null, passed to set.seed. |

`nchains` |
Integer scalar. Number of chains to run (in parallel). |

`burnin` |
Integer scalar. Length of burn-in. Passed to
coda::mcmc as |

`thin` |
Integer scalar. Passed to coda::mcmc. |

`kernel` |
An object of class fmcmc_kernel. |

`multicore` |
Logical. If |

`conv_checker` |
A function that receives an object of class coda::mcmc.list,
and returns a logical value with |

`cl` |
A |

`progress` |
Logical scalar. When set to |

`chain_id` |
Integer scalar (internal use only). This is an argument
passed to the kernel function and it allows it identify in which of the
chains the process is taking place. This could be relevant for some kernels
(see |

This function implements Markov Chain Monte Carlo (MCMC) using the Metropolis-Hastings ratio with flexible transition kernels. Users can specify either one of the available transition kernels or define one of their own (see kernels). Furthermore, it allows easy parallel implementation running multiple chains in parallel. The function also allows using convergence diagnostics tests to set-up a criterion for automatically stopping the algorithm (see convergence-checker).

The canonical form of the Metropolis Hastings algorithm consists on accepting
a move from state `x`

to state `y`

based on the Hastings ratio `r(x,y)`

:

```
%
r(x,y) = \frac{h(y)q(y,x)}{h(x)q(x,y)},%
```

where `h`

is the unnormalized density of the specified distribution (
the posterior probability), and `q`

has the conditional probability of
moving from state `x`

to `y`

(the proposal density). The move
`x \to y`

is then accepted with probability

```
%
\alpha(x,y) = \min\left(1, r(x,y)\right)%
```

Observe that, in the case that `q()`

is symmetric, meaning `q(x, y) = q(y, x)`

,
the Hastings ration reduces to `h(y)/h(x)`

. Starting version 0.5-0, the value
of the log unnormalized density and the proposed states `y`

can be accessed using
the functions `get_logpost()`

and `get_draws()`

.

We now give details of the various options included in the function.

The function `MCMC_without_conv_checker`

is for internal use
only.

`MCMC`

returns an object of class coda::mcmc from the coda
package. The `mcmc`

object is a matrix with one column per parameter,
and `nsteps`

rows. If `nchains > 1`

, then it returns a coda::mcmc.list.

While the main output of `MCMC`

is the `mcmc`

(`.list`

) object, other information
and intermediate outputs of the process are stored in `MCMC_OUTPUT`

. For information
about how to retrieve/set data, see mcmc-output.

By default, if `initial`

is of class `mcmc`

, `MCMC`

will take the last `nchains`

points from the chain as starting point for the new sequence. If `initial`

is
of class `mcmc.list`

, the number of chains in `initial`

must match the `nchains`

parameter.

If `initial`

is a vector, then it must be of length equal to the number of
parameters used in the model. When using multiple chains, if `initial`

is not
an object of class `mcmc`

or `mcmc.list`

, then it must be a numeric matrix
with as many rows as chains, and as many columns as parameters in the model.

When `nchains > 1`

, the function will run multiple chains. Furthermore,
if `cl`

is not passed, `MCMC`

will create a `PSOCK`

cluster
using `makePSOCKcluster`

with
parallel::detectCores
clusters and attempt to execute using multiple cores. Internally, the function does
the following:

# Creating the cluster ncores <- parallel::detectCores() ncores <- ifelse(nchains < ncores, nchains, ncores) cl <- parallel::makePSOCKcluster(ncores) # Loading the package and setting the seed using clusterRNGStream invisible(parallel::clusterEvalQ(cl, library(fmcmc))) parallel::clusterSetRNGStream(cl, .Random.seed)

When running in parallel, objects that are
used within `fun`

must be passed through `...`

, otherwise the cluster
will return with an error.

The user controls the initial value of the parameters of the MCMC algorithm
using the argument `initial`

. When using multiple chains, i.e., `nchains > 1`

,
the user can specify multiple starting points, which is recommended. In such a
case, each row of `initial`

is use as a starting point for each of the
chains. If `initial`

is a vector and `nchains > 1`

, the value is recycled, so
all chains start from the same point (not recommended, the function throws a
warning message).

By default, no automatic stop is implemented. If one of the functions in
convergence-checker is used, then the MCMC is done by bulks as specified
by the convergence checker function, and thus the algorithm will stop if,
the `conv_checker`

returns `TRUE`

. For more information see convergence-checker.

Brooks, S., Gelman, A., Jones, G. L., & Meng, X. L. (2011). Handbook of Markov Chain Monte Carlo. Handbook of Markov Chain Monte Carlo.

Vega Yon, G., & Marjoram, P. (2019). fmcmc: A friendly MCMC framework. Journal of Open Source Software, 4(39), 1427. \Sexpr[results=rd]{tools:::Rd_expr_doi("10.21105/joss.01427")}

`get_logpost()`

, `get_logpost()`

(mcmc-output) for post execution of `MCMC`

, and
`ith_step()`

for accessing objects within an `MCMC`

call.

```
# Univariate distributed data with multiple parameters ----------------------
# Parameters
set.seed(1231)
n <- 1e3
pars <- c(mean = 2.6, sd = 3)
# Generating data and writing the log likelihood function
D <- rnorm(n, pars[1], pars[2])
fun <- function(x) {
x <- log(dnorm(D, x[1], x[2]))
sum(x)
}
# Calling MCMC, but first, loading the coda R package for
# diagnostics
library(coda)
ans <- MCMC(
fun, initial = c(mu=1, sigma=1), nsteps = 2e3,
kernel = kernel_normal_reflective(scale = .1, ub = 10, lb = 0)
)
# Ploting the output
oldpar <- par(no.readonly = TRUE)
par(mfrow = c(1,2))
boxplot(as.matrix(ans),
main = expression("Posterior distribution of"~mu~and~sigma),
names = expression(mu, sigma), horizontal = TRUE,
col = blues9[c(4,9)],
sub = bquote(mu == .(pars[1])~", and"~sigma == .(pars[2]))
)
abline(v = pars, col = blues9[c(4,9)], lwd = 2, lty = 2)
plot(apply(as.matrix(ans), 1, fun), type = "l",
main = "LogLikelihood",
ylab = expression(L("{"~mu,sigma~"}"~"|"~D))
)
par(oldpar)
# In this example we estimate the parameter for a dataset with ----------------
# With 5,000 draws from a MVN() with parameters M and S.
# Loading the required packages
library(mvtnorm)
library(coda)
# Parameters and data simulation
S <- cbind(c(.8, .2), c(.2, 1))
M <- c(0, 1)
set.seed(123)
D <- rmvnorm(5e3, mean = M, sigma = S)
# Function to pass to MCMC
fun <- function(pars) {
# Putting the parameters in a sensible way
m <- pars[1:2]
s <- cbind( c(pars[3], pars[4]), c(pars[4], pars[5]) )
# Computing the unnormalized log likelihood
sum(log(dmvnorm(D, m, s)))
}
# Calling MCMC
ans <- MCMC(
initial = c(mu0=5, mu1=5, s0=5, s01=0, s2=5),
fun,
kernel = kernel_normal_reflective(
lb = c(-10, -10, .01, -5, .01),
ub = 5,
scale = 0.01
),
nsteps = 1e3,
thin = 10,
burnin = 5e2
)
# Checking out the outcomes
plot(ans)
summary(ans)
# Multiple chains -----------------------------------------------------------
# As we want to run -fun- in multiple cores, we have to
# pass -D- explicitly (unless using Fork Clusters)
# just like specifying that we are calling a function from the
# -mvtnorm- package.
fun <- function(pars, D) {
# Putting the parameters in a sensible way
m <- pars[1:2]
s <- cbind( c(pars[3], pars[4]), c(pars[4], pars[5]) )
# Computing the unnormalized log likelihood
sum(log(mvtnorm::dmvnorm(D, m, s)))
}
# Two chains
ans <- MCMC(
initial = c(mu0=5, mu1=5, s0=5, s01=0, s2=5),
fun,
nchains = 2,
kernel = kernel_normal_reflective(
lb = c(-10, -10, .01, -5, .01),
ub = 5,
scale = 0.01
),
nsteps = 1e3,
thin = 10,
burnin = 5e2,
D = D
)
summary(ans)
```

Embedding an R snippet on your website

Add the following code to your website.

For more information on customizing the embed code, read Embedding Snippets.