mgnk: The multivariate G&K example

In BSL: Bayesian Synthetic Likelihood

The multivariate G&K example

Description

Here we provide the data and tuning parameters required to reproduce the results from the multivariate G & K \insertCiteDrovandi2011BSL example from \insertCiteAn2019;textualBSL.

Usage

data(mgnk)

mgnk_sim(theta_tilde, TT, J, bound)

mgnk_sum(y)

Arguments

theta_tilde	A vector with 15 elements for the proposed model parameters.
TT	The number of observations in the data.
J	The number of variables in the data.
bound	A matrix of boundaries for the uniform prior.
y	A TT \times J matrix of data.

Details

It is not practical to give a reasonable explanation of this example through R documentation given the number of equations involved. We refer the reader to the BSLasso paper \insertCiteAn2019BSL at <doi:10.1080/10618600.2018.1537928> for information on the model and summary statistic used in this example.

An example dataset

We use the foreign currency exchange data available from https://www.rba.gov.au/statistics/historical-data.html as in \insertCiteAn2019;textualBSL.

• data: A 1651 \times 3 matrix of data.

• sim_args: Values of sim_args relevant to this example.

• start: A vector of suitable initial values of the parameters for MCMC.

• cov: The covariance matrix of a multivariate normal random walk proposal distribution used in the MCMC, in the form of a 15 by 15 matrix

Author(s)

Ziwen An, Leah F. South and Christopher Drovandi

References

\insertAllCited

Examples

## Not run: 
require(doParallel) # You can use a different package to set up the parallel backend
require(MASS)
require(elliplot)

# Loading the data for this example
data(mgnk)
model <- newModel(fnSim = mgnk_sim, fnSum = mgnk_sum, simArgs = mgnk$sim_args, theta0 = mgnk$start,
    thetaNames = expression(a[1],b[1],g[1],k[1],a[2],b[2],g[2],k[2],
                            a[3],b[3],g[3],k[3],delta[12],delta[13],delta[23]))

# Performing BSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl <- makeCluster(min(detectCores() - 1,2))
registerDoParallel(cl)
resultMgnkBSL <- bsl(mgnk$data, n = 60, M = 80000, model = model, covRandWalk = mgnk$cov,
    method = "BSL", parallel = FALSE, verbose = 1L, plotOnTheFly = TRUE)
stopCluster(cl)
registerDoSEQ()
show(resultMgnkBSL)
summary(resultMgnkBSL)
plot(resultMgnkBSL, which = 2, thin = 20)

# Performing uBSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl <- makeCluster(min(detectCores() - 1,2))
registerDoParallel(cl)
resultMgnkuBSL <- bsl(mgnk$data, n = 60, M = 80000, model = model, covRandWalk = mgnk$cov,
    method = "uBSL", parallel = FALSE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultMgnkuBSL)
summary(resultMgnkuBSL)
plot(resultMgnkuBSL, which = 2, thin = 20)


# Performing tuning for BSLasso
ssy <- mgnk_sum(mgnk$data)
lambda_all <- list(exp(seq(-2.5,0.5,length.out=20)), exp(seq(-2.5,0.5,length.out=20)),
                   exp(seq(-4,-0.5,length.out=20)), exp(seq(-5,-2,length.out=20)))

# Opening up the parallel pools using doParallel
cl <- makeCluster(min(detectCores() - 1,2))
registerDoParallel(cl)
set.seed(100)
sp_mgnk <- selectPenalty(ssy, n = c(15, 20, 30, 50), lambda = lambda_all, theta = mgnk$start,
    M = 100, sigma = 1.5, model = model, method = "BSL", shrinkage = "glasso", standardise = TRUE,
    parallelSim = TRUE, parallelSimArgs = list(.packages = "MASS", .export = "ninenum"),
    parallelMain = TRUE)
stopCluster(cl)
registerDoSEQ()
sp_mgnk
plot(sp_mgnk)

# Performing BSLasso with a fixed penalty (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl <- makeCluster(min(detectCores() - 1,2))
registerDoParallel(cl)
resultMgnkBSLasso <- bsl(mgnk$data, n = 20, M = 80000, model = model, covRandWalk = mgnk$cov,
    method = "BSL", shrinkage = "glasso", penalty = 0.3, standardise = TRUE, parallel = FALSE,
    verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultMgnkBSLasso)
summary(resultMgnkBSLasso)
plot(resultMgnkBSLasso, which = 2, thin = 20)


# Performing semiBSL (reduce the number of iterations M if desired)
# Opening up the parallel pools using doParallel
cl <- makeCluster(min(detectCores() - 1,2))
registerDoParallel(cl)
resultMgnkSemiBSL <- bsl(mgnk$data, n = 60, M = 80000, model = model, covRandWalk = mgnk$cov,
    method = "semiBSL", parallel = FALSE, verbose = 1L)
stopCluster(cl)
registerDoSEQ()
show(resultMgnkSemiBSL)
summary(resultMgnkSemiBSL)
plot(resultMgnkSemiBSL, which = 2, thin = 20)

# Plotting the results together for comparison
# plot using the R default plot function
oldpar <- par()
par(mar = c(4, 4, 1, 1), oma = c(0, 1, 2, 0))
combinePlotsBSL(list(resultMgnkBSL, resultMgnkuBSL, resultMgnkBSLasso, resultMgnkSemiBSL),
                which = 1, thin = 20, label = c("bsl", "ubsl", "bslasso", "semiBSL"),
                col = c("red", "yellow", "blue", "green"), lty = 2:5, lwd = 1)
mtext("Approximate Univariate Posteriors", outer = TRUE, line = 0.75, cex = 1.2)

# plot using the ggplot2 package
combinePlotsBSL(list(resultMgnkBSL, resultMgnkuBSL, resultMgnkBSLasso, resultMgnkSemiBSL),
    which = 2, thin = 20, label=c("bsl","ubsl","bslasso","semiBSL"),
    options.color=list(values=c("red","yellow","blue","green")),
    options.linetype = list(values = 2:5), options.size = list(values = rep(1, 4)),
    options.theme = list(plot.margin = grid::unit(rep(0.03,4),"npc"),
        axis.title = ggplot2::element_text(size=12), axis.text = ggplot2::element_text(size = 8),
        legend.text = ggplot2::element_text(size = 12)))
par(mar = oldpar$mar, oma = oldpar$oma)

## End(Not run)

BSL documentation built on Nov. 3, 2022, 9:06 a.m.