From ad707902f9cac15ec7f781c79bd3c2a88ec23ae4 Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Tue, 30 Jan 2024 14:52:51 +1100 Subject: [PATCH] Added a var.mcmc.list() pseudo-method and cleaned up documentation for the others. --- DESCRIPTION | 4 ++-- NAMESPACE | 1 + R/mcmc-utils.R | 40 +++++++++++++++++++++++++++++++++------- man/mcmc-utilities.Rd | 32 +++++++++++++++++++++++++------- 4 files changed, 61 insertions(+), 16 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 2e8f86b..a039370 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: statnet.common -Version: 4.10.0-439 -Date: 2023-12-18 +Version: 4.10.0-441 +Date: 2024-01-30 Title: Common R Scripts and Utilities Used by the Statnet Project Software Authors@R: c( person(c("Pavel", "N."), "Krivitsky", role=c("aut","cre"), email="pavel@statnet.org", comment=c(ORCID="0000-0002-9101-3362", affiliation="University of New South Wales")), diff --git a/NAMESPACE b/NAMESPACE index ffb7833..2f76992 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -138,6 +138,7 @@ export(ult) export(unused_dots_warning) export(unwhich) export(update_snctrl) +export(var.mcmc.list) export(vector.namesmatch) export(xAxT) export(xTAx) diff --git a/R/mcmc-utils.R b/R/mcmc-utils.R index b0308c2..d40075b 100644 --- a/R/mcmc-utils.R +++ b/R/mcmc-utils.R @@ -10,7 +10,7 @@ #' @name mcmc-utilities #' @title Utility operations for [`mcmc.list`] objects #' -#' @description \code{colMeans.mcmc.list} is a "method" for (non-generic) [`colMeans`] applicable to [`mcmc.list`] objects. +#' @description \code{colMeans.mcmc.list} is a "method" for (non-generic) [colMeans()] applicable to [`mcmc.list`] objects. #' #' @param x a \code{\link{mcmc.list}} object. #' @param \dots additional arguments to \code{\link{colMeans}} or @@ -18,7 +18,7 @@ #' @return \code{colMeans.mcmc} returns a vector with length equal to #' the number of mcmc chains in \code{x} with the mean value for #' each chain. -#' @seealso [`colMeans`], [`mcmc.list`] +#' @seealso [colMeans()], [`mcmc.list`] #' @examples #' data(line, package="coda") #' summary(line) # coda @@ -29,15 +29,41 @@ #' @export colMeans.mcmc.list colMeans.mcmc.list<-function(x,...) colMeans(as.matrix(x),...) +#' @rdname mcmc-utilities +#' +#' @description \code{var.mcmc.list} is a "method" for (non-generic) +#' [var()] applicable to [`mcmc.list`] objects. Since MCMC chains +#' are assumed to all be sampling from the same underlying +#' distribution, pooled mean is used. This implementation should be +#' equivalent (within numerical error) to `var(as.matrix(x))` while +#' avoiding constructing the large matrix. +#' +#' @seealso [var()] +#' @examples +#' data(line, package="coda") +#' var(as.matrix(line)) # coda +#' var.mcmc.list(line) # "Method" +#' \dontshow{ +#' stopifnot(isTRUE(all.equal(var.mcmc.list(line), var(as.matrix(line))))) +#' } +#' @export var.mcmc.list +var.mcmc.list <- function(x, ...){ + nchain <- length(x) + niter <- NROW(x[[1]]) + SSW <- Reduce(`+`, lapply(x, cov)) * (niter-1) + SSB <- cov(t(sapply(x, colMeans))) * niter + (SSW+SSB) / (niter*nchain-1) +} + #' @rdname mcmc-utilities #' #' @description \code{sweep.mcmc.list} is a "method" for (non-generic) -#' [`sweep`] applicable to [`mcmc.list`] objects. +#' [sweep()] applicable to [`mcmc.list`] objects. #' -#' @param STATS,FUN,check.margin See help for [`sweep`]. +#' @param STATS,FUN,check.margin See help for [sweep()]. #' @return \code{sweep.mcmc.list} returns an appropriately modified #' version of \code{x} -#' @seealso [`sweep`] +#' @seealso [sweep()] #' @examples #' data(line, package="coda") #' colMeans.mcmc.list(line)-1:3 @@ -56,12 +82,12 @@ sweep.mcmc.list<-function(x, STATS, FUN="-", check.margin=TRUE, ...){ #' @rdname mcmc-utilities #' #' @description \code{lapply.mcmc.list} is a "method" for (non-generic) -#' [`lapply`] applicable to [`mcmc.list`] objects. +#' [lapply()] applicable to [`mcmc.list`] objects. #' #' @param X An [`mcmc.list`] object. #' @return `lapply.mcmc.list` returns an [`mcmc.list`] each of #' whose chains had been passed through `FUN`. -#' @seealso [`lapply`] +#' @seealso [lapply()] #' @examples #' data(line, package="coda") #' colMeans.mcmc.list(line)[c(2,3,1)] diff --git a/man/mcmc-utilities.Rd b/man/mcmc-utilities.Rd index 2e11603..55deec3 100644 --- a/man/mcmc-utilities.Rd +++ b/man/mcmc-utilities.Rd @@ -3,12 +3,15 @@ \name{mcmc-utilities} \alias{mcmc-utilities} \alias{colMeans.mcmc.list} +\alias{var.mcmc.list} \alias{sweep.mcmc.list} \alias{lapply.mcmc.list} \title{Utility operations for \code{\link{mcmc.list}} objects} \usage{ colMeans.mcmc.list(x, ...) +var.mcmc.list(x, ...) + sweep.mcmc.list(x, STATS, FUN = "-", check.margin = TRUE, ...) lapply.mcmc.list(X, FUN, ...) @@ -19,7 +22,7 @@ lapply.mcmc.list(X, FUN, ...) \item{\dots}{additional arguments to \code{\link{colMeans}} or \code{\link{sweep}}.} -\item{STATS, FUN, check.margin}{See help for \code{\link{sweep}}.} +\item{STATS, FUN, check.margin}{See help for \code{\link[=sweep]{sweep()}}.} \item{X}{An \code{\link{mcmc.list}} object.} } @@ -35,13 +38,20 @@ version of \code{x} whose chains had been passed through \code{FUN}. } \description{ -\code{colMeans.mcmc.list} is a "method" for (non-generic) \code{\link{colMeans}} applicable to \code{\link{mcmc.list}} objects. +\code{colMeans.mcmc.list} is a "method" for (non-generic) \code{\link[=colMeans]{colMeans()}} applicable to \code{\link{mcmc.list}} objects. + +\code{var.mcmc.list} is a "method" for (non-generic) +\code{\link[=var]{var()}} applicable to \code{\link{mcmc.list}} objects. Since MCMC chains +are assumed to all be sampling from the same underlying +distribution, pooled mean is used. This implementation should be +equivalent (within numerical error) to \code{var(as.matrix(x))} while +avoiding constructing the large matrix. \code{sweep.mcmc.list} is a "method" for (non-generic) -\code{\link{sweep}} applicable to \code{\link{mcmc.list}} objects. +\code{\link[=sweep]{sweep()}} applicable to \code{\link{mcmc.list}} objects. \code{lapply.mcmc.list} is a "method" for (non-generic) -\code{\link{lapply}} applicable to \code{\link{mcmc.list}} objects. +\code{\link[=lapply]{lapply()}} applicable to \code{\link{mcmc.list}} objects. } \examples{ data(line, package="coda") @@ -51,6 +61,12 @@ colMeans.mcmc.list(line) # "Method" stopifnot(isTRUE(all.equal(summary(line)$statistics[,"Mean"],colMeans.mcmc.list(line)))) } data(line, package="coda") +var(as.matrix(line)) # coda +var.mcmc.list(line) # "Method" +\dontshow{ +stopifnot(isTRUE(all.equal(var.mcmc.list(line), var(as.matrix(line))))) +} +data(line, package="coda") colMeans.mcmc.list(line)-1:3 colMeans.mcmc.list(sweep.mcmc.list(line, 1:3)) \dontshow{ @@ -64,9 +80,11 @@ stopifnot(isTRUE(all.equal(colMeans.mcmc.list(line)[c(2,3,1)],colMeans.mcmc.list } } \seealso{ -\code{\link{colMeans}}, \code{\link{mcmc.list}} +\code{\link[=colMeans]{colMeans()}}, \code{\link{mcmc.list}} + +\code{\link[=var]{var()}} -\code{\link{sweep}} +\code{\link[=sweep]{sweep()}} -\code{\link{lapply}} +\code{\link[=lapply]{lapply()}} }