You can use these functions to read variables, store, and retrieve data during the MCMC process.
ith_step(x)
set_userdata(...)
get_userdata()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:
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.
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
# ////////////////////////////////////////////////////
# 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
# 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())
#> i sum_of_p
#> 1 1 71.00000
#> 2 2 71.00627
#> 3 3 70.97390
#> 4 4 70.97180
#> 5 5 71.02942
#> 6 6 70.95732
# 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
)
#> Warning: While using multiple chains, a single initial point has been passed via `initial`: c(70, -2, 2, 1). The values will be recycled. Ideally you would want to start each chain from different locations.
user_dat <- get_userdata()
lapply(user_dat, head)
#> [[1]]
#> i sum_of_p
#> 1 1 71.00000
#> 2 2 70.95107
#> 3 3 70.66570
#> 4 4 70.65513
#> 5 5 70.87995
#> 6 6 70.61189
#>
#> [[2]]
#> i sum_of_p
#> 1 1 71.00000
#> 2 2 71.25754
#> 3 3 70.62444
#> 4 4 71.38458
#> 5 5 70.87700
#> 6 6 71.01785
#>