The core functions of binaryRL have already been
rewritten in Rcpp. This optimization reduces the runtime to
roughly one-tenth of the previous R-based for-loop implementation. As a
result, running MCMC with LaplacesDemon on
binaryRL is very fast, while still allowing users to easily
define and customize their reinforcement learning models.
Example Data
test_data <- binaryRL::Mason_2024_G2
# Ensure the Subject ID is a numeric value for easier indexing
test_data$Subject <- as.numeric(test_data$Subject)
# Filter the dataset to include only subjects 1 through 5
test_data <- test_data[test_data$Subject <= 5, ]
unique(test_data$Subject)## [1] 1 2 3 4 5Object Function
Argument 1: parm
Similar to other frameworks designed for individual-level analysis,
the primary argument to the LaplacesDemon model function is
a single vector of free parameters. Crucially, this vector does not
represent the parameters for a single “black-box” function call.
Instead, it is a concatenated vector containing the complete set of free
parameters for every subject being modeled. Consequently, its total
length is the number of free parameters per subject multiplied by the
number of subjects.
Argument 2: Data
A key convenience of LaplacesDemon is that the
user-defined Model function can accept a Data argument.
This feature provides a clean and explicit mechanism for passing
supplementary information—such as prior distributions, trial counts, or
policy settings—directly to the model. Consequently, users can avoid
more cumbersome methods, like environment binding, and instead bundle
all necessary data and constants into a single, well-organized list.
Data <- list(
# The prepared dataframe containing the experimental data
df = test_data,
# A vector of unique subject IDs to loop over during modeling
id = unique(test_data$Subject),
# Need to constrain the free parameters with Laplace
priors = NULL,
# The number of free parameters per subject (eta and tau)
n_params = 2,
# The total number of trials per subject in the experiment
n_trials = 360,
# The operational mode for the binaryRL::run_m function, set to "fit"
mode = "fit",
# The learning policy for binaryRL::run_m, set to "off-policy"
policy = "off",
# A vector of descriptive names for each parameter across all subjects
parm.names = LaplacesDemon::as.parm.names(
list(
# Creates names like "[1,1]", "[1,2]", "[2,1]", "[2,2]", etc.
outer(
X = c("eta", "tau"),
Y = 1:length(unique(test_data$Subject)),
FUN = paste, sep = ""
)
)
),
# The names of variables to monitor during the MCMC run (Log Posterior)
mon.names = "LP",
# The total number of subjects
N = length(unique(test_data$Subject))
)Defining Object Function
Create a model function that accepts parm and
Data as its arguments, with a return value that conforms to
the structure required by LaplacesDemon. Within this
function, you need to calculate the sum of the log-likelihoods and the
sum of the log-posteriors across all subjects by using specified prior
distributions.
Model <- function(parm, Data) {
# Initialize accumulators for sum of log-likelihood and log-posterior
sum_LL <- 0
sum_LP <- 0
# Loop through each subject ID provided in the Data list
for (i in as.vector(Data$id)) {
# read unconstrained params
eta_unconstrained <- parm[(i - 1) * 2 + 1]
tau_unconstrained <- parm[(i - 1) * 2 + 2]
# convert to constrained params
eta <- stats::plogis(eta_unconstrained)
tau <- exp(tau_unconstrained)
# pass constrained params to the input vector
#parm[(i - 1) * 2 + 1] <- eta
#parm[(i - 1) * 2 + 2] <- tau
res <- binaryRL::run_m(
data = Data$df,
id = Data$id[i],
eta = eta,
tau = tau,
priors = Data$priors,
n_params = Data$n_params,
n_trials = Data$n_trials,
mode = Data$mode,
policy = Data$policy
)
# Add the log-likelihood of the current subject to the total
sum_LL <- sum_LL + res$ll
# Add the log-posterior of the current subject to the total
LP_eta <- stats::dnorm(x = eta, mean = 0.5, sd = 0.15, log = TRUE)
LP_tau <- stats::dexp(x = tau, rate = 1, log = TRUE)
sum_LP <- sum_LP + LP_eta + LP_tau
}
LP <- sum_LL + sum_LP
# Return a list in the format required by LaplacesDemon
return(
list(
LP = LP,
Dev = -2 * sum_LL,
Monitor = c(LP = LP),
yhat = NULL,
parm = parm
)
)
}Test Object Function
Initial.Values <- stats::runif(Data$N * Data$n_params)
Model(parm = Initial.Values, Data = Data)## $LP
## [1] -1635.234
##
## $Dev
## [1] 3260.12
##
## $Monitor
## LP
## -1635.234
##
## $yhat
## NULL
##
## $parm
## [1] 0.080750138 0.834333037 0.600760886 0.157208442 0.007399441 0.466393497
## [7] 0.497777389 0.289767245 0.732881987 0.772521511Run LaplacesDemon
Fit <- LaplacesDemon::LaplacesDemon(
Model = Model,
Data = Data,
Packages = "binaryRL",
# Set initial values (# of free params * # of subjects)
Initial.Values = stats::runif(length(Data$id) * Data$n_params),
# The relationship between `Algorithm` and `Specs` is special
# Refer to the LaplacesDemon help documents for details
Algorithm = "MWG", Specs = list(B = NULL),
#Algorithm = "NUTS", Specs = list(A = 5, delta = 0.6, epsilon = NULL, Lmax = Inf),
# The total number of iterations to perform in each chain.
Iterations = 10,
# The thinning interval. A value of n keeps one sample every n-th iteration.
Thinning = 1,
# Logs for every ? iterations
Status = 1
)
LaplacesDemon::Consort(Fit)#############################################################
# Consort with Laplace's Demon #
#############################################################
Call:
LaplacesDemon::LaplacesDemon(Model = Model, Data = Data, Initial.Values = stats::runif(length(Data$id) *
Data$n_params), Iterations = 10, Status = 1, Thinning = 1,
Algorithm = "MWG", Specs = list(B = NULL), Packages = "binaryRL")
Acceptance Rate: 0.29091
Algorithm: Metropolis-within-Gibbs
Covariance Matrix: (NOT SHOWN HERE; diagonal shown instead)
[1,1] [2,1] [1,2] [2,2] [1,3] [2,3] [1,4] [2,4] [1,5] [2,5]
0.5670132 0.5670132 0.5670132 0.5670132 0.5670132 0.5670132 0.5670132 0.5670132 0.5670132 0.5670132
Covariance (Diagonal) History: (NOT SHOWN HERE)
Deviance Information Criterion (DIC):
All Stationary
Dbar 3355.960 3355.960
pD 378.877 378.877
DIC 3734.837 3734.837
Initial Values:
[1] 0.09859165 0.69300856 0.75132556 0.99870695 0.34549739 0.29300554 0.66773598 0.59801099 0.70105000 0.10186876
Iterations: 11
Log(Marginal Likelihood): NA
Minutes of run-time: 1.39
Model: (NOT SHOWN HERE)
Monitor: (NOT SHOWN HERE)
Parameters (Number of): 10
Posterior1: (NOT SHOWN HERE)
Posterior2: (NOT SHOWN HERE)
Recommended Burn-In of Thinned Samples: 0
Recommended Burn-In of Un-thinned Samples: 0
Recommended Thinning: 10
Specs: (NOT SHOWN HERE)
Status is displayed every 1 iterations
Summary1: (SHOWN BELOW)
Summary2: (SHOWN BELOW)
Thinned Samples: 11
Thinning: 1
Summary of All Samples
Mean SD MCSE ESS LB Median UB
[1,1] -6.326389e-01 0.00000000 0.00000000 2.220446e-16 -6.326389e-01 -6.326389e-01 -0.6326389
[2,1] 6.930086e-01 0.00000000 0.00000000 2.220446e-16 6.930086e-01 6.930086e-01 0.6930086
[1,2] -7.464188e-03 0.19969803 0.09095104 1.100000e+01 -9.168507e-02 -9.168507e-02 0.4487322
[2,2] 9.987070e-01 0.00000000 0.00000000 2.220446e-16 9.987070e-01 9.987070e-01 0.9987070
[1,3] 6.445531e-01 0.00000000 0.00000000 2.220446e-16 6.445531e-01 6.445531e-01 0.6445531
[2,3] -5.424715e-01 0.31279239 0.14245907 1.100000e+01 -6.743889e-01 -6.743889e-01 0.1720812
[1,4] 6.677360e-01 0.00000000 0.00000000 2.220446e-16 6.677360e-01 6.677360e-01 0.6677360
[2,4] -2.994002e-01 0.30341356 0.22470223 2.220446e-16 -7.339427e-01 -2.994002e-01 0.1351423
[1,5] 7.705329e-01 0.04985889 0.03392331 3.014025e+00 6.831834e-01 8.100874e-01 0.8100874
[2,5] 4.602979e-01 0.59291341 0.43910024 2.220446e-16 -3.888602e-01 4.602979e-01 1.3094560
Deviance 3.355960e+03 27.52734495 14.79642686 3.718372e+00 3.330060e+03 3.349300e+03 3415.7700000
LP -1.683476e+03 13.59142083 7.03639451 3.903581e+00 -1.713494e+03 -1.679777e+03 -1671.8078968
Summary of Stationary Samples
Mean SD MCSE ESS LB Median UB
[1,1] -6.326389e-01 0.00000000 0.00000000 2.220446e-16 -6.326389e-01 -6.326389e-01 -0.6326389
[2,1] 6.930086e-01 0.00000000 0.00000000 2.220446e-16 6.930086e-01 6.930086e-01 0.6930086
[1,2] -7.464188e-03 0.19969803 0.09095104 1.100000e+01 -9.168507e-02 -9.168507e-02 0.4487322
[2,2] 9.987070e-01 0.00000000 0.00000000 2.220446e-16 9.987070e-01 9.987070e-01 0.9987070
[1,3] 6.445531e-01 0.00000000 0.00000000 2.220446e-16 6.445531e-01 6.445531e-01 0.6445531
[2,3] -5.424715e-01 0.31279239 0.14245907 1.100000e+01 -6.743889e-01 -6.743889e-01 0.1720812
[1,4] 6.677360e-01 0.00000000 0.00000000 2.220446e-16 6.677360e-01 6.677360e-01 0.6677360
[2,4] -2.994002e-01 0.30341356 0.22470223 2.220446e-16 -7.339427e-01 -2.994002e-01 0.1351423
[1,5] 7.705329e-01 0.04985889 0.03392331 3.014025e+00 6.831834e-01 8.100874e-01 0.8100874
[2,5] 4.602979e-01 0.59291341 0.43910024 2.220446e-16 -3.888602e-01 4.602979e-01 1.3094560
Deviance 3.355960e+03 27.52734495 14.79642686 3.718372e+00 3.330060e+03 3.349300e+03 3415.7700000
LP -1.683476e+03 13.59142083 7.03639451 3.903581e+00 -1.713494e+03 -1.679777e+03 -1671.8078968
Demonic Suggestion
Due to the combination of the following conditions,
1. Metropolis-within-Gibbs
2. The acceptance rate (0.2909091) is within the interval [0.15,0.5].
3. At least one target MCSE is >= 6.27% of its marginal posterior
standard deviation.
4. At least one target distribution has an effective sample size
(ESS) less than 100. The worst mixing chain is: [1,1] (ESS=2.220446e-16).
5. Each target distribution became stationary by
1 iteration.
Laplace's Demon has not been appeased, and suggests
copy/pasting the following R code into the R console,
and running it.
Initial.Values <- as.initial.values(Fit)
Fit <- LaplacesDemon(Model, Data=Data, Initial.Values,
Covar=Fit$Covar, Iterations=110, Status=7, Thinning=10,
Algorithm="AMWG", Specs=list(B=NULL, n=11, Periodicity=34))
Laplace's Demon is finished consorting.