fmcmc (version 0.5-2)

mcmc-output: Information about the last MCMC call

Description

This environment holds a copy of the last call to MCMC, including the start and end time (to compute total elapsed time) of the call. Since the resulting object of MCMC is an object of class coda::mcmc, this is a way to capture more information in case the user needs it.

Usage

get_(x)
get_logpost()
get_draws()
get_elapsed()
get_initial()
get_fun()
get_nsteps()
get_seed()
get_nchains()
get_burnin()
get_thin()
get_kernel()
get_multicore()
get_conv_checker()
get_cl()
get_progress()
get_chain_id()

Arguments

x

Character scalar. Name of an argument to retrieve. If x was not passed to the last call, the function returns with an error.

Details

The function get_logpost returns the logposterior value at each iteration. The values correspond to a named numeric vector. If nchains > 1 then it will return a list of length nchains with the corresponding logpost values for each chain.

The function get_draws() retrieves the proposed states from the kernel function.

Examples

# Getting the logpost -------------------------------------------------------
set.seed(23133)
x <- rnorm(200)
y <- -4 + x*2 + rnorm(200)
f <- function(p) {
  sum(dnorm(y - p[1] - x*p[2], log = TRUE))
}

# Setting a RAM kernel
kern <- kernel_am(eps = 1e-2)

ans <- MCMC(fun = f, initial = c(0, 1), nsteps = 2000, kernel = kern)
plot(
  # Plotting the logpost from the last run
  -get_logpost(), 
  # Getting the number of chains
  main = paste0("nchains: ", get_nchains()),
  
  # And the elapsed time
  sub  = sprintf("Run time: %.4f(s)", get_elapsed()[3]),
  type = "l",
  log = "y"
) 

# This also works using multiple chains
ans <- MCMC(fun = f, initial = c(0, 0), nsteps=2000, nchains = 2, kernel = kern)

# In this case, just like -ans-, 
draws <- get_draws()

# Plotting proposed points vs accepted
plot(
  draws[[1]], pch = 20,
  col = adjustcolor("gray", alpha = .5),
  main = "Accepted vs proposed states\n(chain 1)"
  )
lines(ans[[1]], pch = 20, col = "tomato", lwd = 2)
legend(
  "topleft", legend = c("Accepted", "Proposed"), pch = c(NA, 20),
  col = c("tomato", "black"), lty = c(1, NA), lwd = c(2, NA)
)