From 733aa116cfd480b6154a4d89f143a02e0c72ad2f Mon Sep 17 00:00:00 2001 From: "Pavel N. Krivitsky" Date: Mon, 23 Sep 2024 13:53:07 +1000 Subject: [PATCH] Implementation and documentation improvements to the colMeans.mcmc.list() family. --- DESCRIPTION | 6 +++--- NAMESPACE | 1 + R/mcmc-utils.R | 30 +++++++++++++++++++----------- man/mcmc-utilities.Rd | 20 ++++++++++++-------- 4 files changed, 35 insertions(+), 22 deletions(-) diff --git a/DESCRIPTION b/DESCRIPTION index 0f81886..d8dee4d 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: statnet.common -Version: 4.10.0-442 -Date: 2024-01-30 +Version: 4.10.0-447 +Date: 2024-09-23 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")), @@ -13,7 +13,7 @@ BugReports: https://github.com/statnet/statnet.common/issues License: GPL-3 + file LICENSE URL: https://statnet.org Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.1 +RoxygenNote: 7.3.2 Encoding: UTF-8 Suggests: covr, rlang (>= 1.1.1), diff --git a/NAMESPACE b/NAMESPACE index 2f76992..0a92d52 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -152,6 +152,7 @@ importFrom(coda,as.mcmc) importFrom(coda,as.mcmc.list) importFrom(methods,is) importFrom(stats,as.formula) +importFrom(stats,var) importFrom(utils,capture.output) importFrom(utils,getAnywhere) importFrom(utils,modifyList) diff --git a/R/mcmc-utils.R b/R/mcmc-utils.R index 60f419e..7e78cea 100644 --- a/R/mcmc-utils.R +++ b/R/mcmc-utils.R @@ -13,30 +13,37 @@ #' @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 -#' \code{\link{sweep}}. +#' @param \dots additional arguments to the functions evaluated on each chain. #' @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`] +#' +#' @details These implementations should be equivalent (within +#' numerical error) to the same function being called on +#' `as.matrix(x)`, while avoiding construction of the large matrix. +#' +#' @seealso [`mcmc.list`] +#' +#' [colMeans()] #' @examples #' data(line, package="coda") -#' summary(line) # coda +#' colMeans(as.matrix(line)) # also coda #' colMeans.mcmc.list(line) # "Method" #' \dontshow{ -#' stopifnot(isTRUE(all.equal(summary(line)$statistics[,"Mean"],colMeans.mcmc.list(line)))) +#' stopifnot(isTRUE(all.equal(colMeans(as.matrix(line)),colMeans.mcmc.list(line)))) #' } #' @export colMeans.mcmc.list -colMeans.mcmc.list<-function(x,...) colMeans(as.matrix(x),...) +colMeans.mcmc.list <- function(x,...){ + nchain <- length(x) + Reduce(`+`, lapply(x, colMeans, ...)) / nchain +} #' @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. +#' distribution, their pooled mean is used. #' #' @seealso [var()] #' @examples @@ -46,12 +53,13 @@ colMeans.mcmc.list<-function(x,...) colMeans(as.matrix(x),...) #' \dontshow{ #' stopifnot(isTRUE(all.equal(var.mcmc.list(line), var(as.matrix(line))))) #' } +#' @importFrom stats var #' @export var.mcmc.list var.mcmc.list <- function(x, ...){ nchain <- length(x) niter <- NROW(x[[1]]) - SSW <- Reduce(`+`, lapply(x, cov)) * (niter-1) - SSB <- if(nchain > 1) cov(t(sapply(x, colMeans))) * niter else 0 + SSW <- Reduce(`+`, lapply(x, var, ...)) * (niter-1) + SSB <- if(nchain > 1) var(t(sapply(x, colMeans, ...))) * niter else 0 (SSW+SSB) / (niter*nchain-1) } diff --git a/man/mcmc-utilities.Rd b/man/mcmc-utilities.Rd index 55deec3..29d9613 100644 --- a/man/mcmc-utilities.Rd +++ b/man/mcmc-utilities.Rd @@ -19,8 +19,7 @@ lapply.mcmc.list(X, FUN, ...) \arguments{ \item{x}{a \code{\link{mcmc.list}} object.} -\item{\dots}{additional arguments to \code{\link{colMeans}} or -\code{\link{sweep}}.} +\item{\dots}{additional arguments to the functions evaluated on each chain.} \item{STATS, FUN, check.margin}{See help for \code{\link[=sweep]{sweep()}}.} @@ -43,9 +42,7 @@ whose chains had been passed through \code{FUN}. \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. +distribution, their pooled mean is used. \code{sweep.mcmc.list} is a "method" for (non-generic) \code{\link[=sweep]{sweep()}} applicable to \code{\link{mcmc.list}} objects. @@ -53,12 +50,17 @@ avoiding constructing the large matrix. \code{lapply.mcmc.list} is a "method" for (non-generic) \code{\link[=lapply]{lapply()}} applicable to \code{\link{mcmc.list}} objects. } +\details{ +These implementations should be equivalent (within +numerical error) to the same function being called on +\code{as.matrix(x)}, while avoiding construction of the large matrix. +} \examples{ data(line, package="coda") -summary(line) # coda +colMeans(as.matrix(line)) # also coda colMeans.mcmc.list(line) # "Method" \dontshow{ -stopifnot(isTRUE(all.equal(summary(line)$statistics[,"Mean"],colMeans.mcmc.list(line)))) +stopifnot(isTRUE(all.equal(colMeans(as.matrix(line)),colMeans.mcmc.list(line)))) } data(line, package="coda") var(as.matrix(line)) # coda @@ -80,7 +82,9 @@ stopifnot(isTRUE(all.equal(colMeans.mcmc.list(line)[c(2,3,1)],colMeans.mcmc.list } } \seealso{ -\code{\link[=colMeans]{colMeans()}}, \code{\link{mcmc.list}} +\code{\link{mcmc.list}} + +\code{\link[=colMeans]{colMeans()}} \code{\link[=var]{var()}}