mcmc-loop: Functions to interact with the main loop

Description Usage Arguments Value Advanced usage Examples

Description

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

Usage

1
2
3
4
5

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.

Value

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

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.

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:

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
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

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
#' # 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)

fmcmc documentation built on Jan. 14, 2022, 9:07 a.m.