fmcmc (version 0.5-2)

mcmc-loop: Functions to interact with the main loop

Description

You can use these functions to read variables, store, and retrieve data during the MCMC process.

Usage

ith_step(x)
set_userdata(...)
get_userdata()

Value

The function ith_step() provides access to the following elements:

• i : (int) Step (iteration) number.

• nsteps : (int) Number of steps.

• chain_id : (int) Id of the chain (goes from 1 to -nchains-)

• theta0 : (double vector) Current state of the chain.

• theta1 : (double vector) Proposed state of the chain.

• ans : (double matrix) Set of accepted states (it will be NA for rows >= i).

• draws : (double matrix) Set of proposed states (it will be NA for rows >= i).

• logpost : (double vector) Value of -fun- (it will be NA for elements >= i).

• R : (double vector) Random values from U(0,1). This is used with the Hastings ratio.

• thin : (int) Thinning (applied after the last step).

• burnin : (int) Burn-in (applied after the last step).

• conv_checker : (function) Convergence checker function.

• kernel : (fmcmc_kernel) Kernel object.

• fun : (function) Passed function to MCMC.

• f : (function) Wrapper of -fun-.

• initial : (double vector) Starting point of the chain.

The following objects always have fixed values (see ?ith_step): nchains, cl, multicore

Other available objects: cnames, funargs, MCMC_OUTPUT, passedargs, progress

The function set_userdata() returns invisible(). The only side effect is appending the information by row.

Arguments

x

Name of the element to retrieve. If missing, it will return the entire environment in which the main MCMC loop is running.

...

Named values to be appended to the user data.

Advanced usage

The function ith_step() is a convenience function that provides access to the environment within which the main loop of the MCMC call is being evaluated. This is a wrapper of MCMC_OUTPUT$loop_envir that will either return the value x or, if missing, the entire environment. If ith_step() is called outside of the MCMC call, then it will return with an error.

For example, if you wanted to print information if the current value of logpost is greater than the previous value of logpost, you could define the objective function as follows:

f <- function(p) {
  i            <- ith_step("i")
  logpost_prev <- ith_step("logpost")[i - 1L]
  logpost_curr <- sum(dnorm(y - x*p, log = TRUE))
  
  if (logpost_prev < logpost_curr)
    cat("At a higher point!\n")
    
  return(logpost_curr)
}

In the case of the objects nchains, cl, and multicore, the function will always return the default values 1, NULL, and FALSE, respectively. Thus, the user shouldn't rely on these objects to provide information regarding runs using multiple chains. More examples below.

Examples

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

ans <- MCMC(fun = f, initial = c(0), nsteps=2000)
plot(get_logpost(), type = "l") # Plotting the logpost from the last run


# Printing information every 500 step ---------------------------------------
# for this we use ith_step()

f <- function(p) {

  # Capturing info from within the loop
  i      <- ith_step("i")
  nsteps <- ith_step("nsteps")
  
  if (!(i %% 500)) {
  
    cat(
      "////////////////////////////////////////////////////\n",
      "Step ", i, " of ", nsteps,". Values in the loop:\n",
      "theta0: ", ith_step("theta0"), "\n",
      "theta1: ", ith_step()$theta1, "\n",
      sep = ""
    )
  }
    

  sum(dnorm(y - x*p, log = TRUE))
}

ans0 <- MCMC(fun = f, initial = c(0), nsteps=2000, progress = FALSE, seed = 22)
# ////////////////////////////////////////////////////
# Step 500 of 2000. Values in the loop:
# theta0: 2.025379
# theta1: 1.04524
# ////////////////////////////////////////////////////
# Step 1000 of 2000. Values in the loop:
# theta0: 2.145967
# theta1: 0.2054037
# ////////////////////////////////////////////////////
# Step 1500 of 2000. Values in the loop:
# theta0: 2.211691
# theta1: 2.515361
# ////////////////////////////////////////////////////
# Step 2000 of 2000. Values in the loop:
# theta0: 1.998789
# theta1: 1.33034


# Printing information if the current logpost is greater than max -----------
f <- function(p) {

  i            <- ith_step("i")
  logpost_prev <- max(ith_step("logpost")[1:(i-1)])
  logpost_curr <- sum(dnorm(y - x*p, log = TRUE))
  
  # Only worthwhile after the first step
  if ((i > 1L) && logpost_prev < logpost_curr)
    cat("At a higher point!:", logpost_curr, ", step:", i,"\n")
    
  return(logpost_curr)

}
ans1 <- MCMC(fun = f, initial = c(0), nsteps=1000, progress = FALSE, seed = 22)
# At a higher point!: -357.3584 , step: 2 
# At a higher point!: -272.6816 , step: 6 
# At a higher point!: -270.9969 , step: 7 
# At a higher point!: -269.8128 , step: 24 
# At a higher point!: -269.7435 , step: 46 
# At a higher point!: -269.7422 , step: 543 
# At a higher point!: -269.7421 , step: 788 
# Saving extra information --------------------------------------------------
data("lifeexpect")

# Defining the logposterior
logpost <- function(p) {

  # Reconding the sum of the parameters (just because) 
  # and the number of step.
  set_userdata(i = ith_step("i"), sum_of_p = sum(p))

  with(lifeexpect, {
    sum(dnorm(age - (p[1] + p[2]*smoke + p[3]*female), sd = p[4], log = TRUE))
  })
  
}

# A kernel where sd is positive, the first is average age, so we 
# make it positive too
kern <- kernel_ram(lb = c(10, -20, -20, .0001), eps = .01)
ans <- MCMC(
  initial = c(70, -2, 2, 1), fun = logpost, kernel = kern, nsteps = 1000, seed = 1
  )

# Retrieving the data
head(get_userdata())

# It also works using multiple chains
ans_two <- MCMC(
  initial = c(70, -2, 2, 1), fun = logpost, kernel = kern, nsteps = 1000, seed = 1, nchains = 2
  )
  
user_dat <- get_userdata()
lapply(user_dat, head)