[r-cran-coda] 03/60: Imported Upstream version 0.9-4
Andreas Tille
tille at debian.org
Fri Dec 16 12:11:21 UTC 2016
This is an automated email from the git hooks/post-receive script.
tille pushed a commit to branch master
in repository r-cran-coda.
commit 56c1069bf40f9e4f91012862bafba3f29208b28e
Author: Andreas Tille <tille at debian.org>
Date: Fri Dec 16 12:09:53 2016 +0100
Imported Upstream version 0.9-4
---
AUTHORS | 2 ++
CHANGELOG | 27 ++++++++++++++++++
DESCRIPTION | 6 ++--
INDEX | 7 +++--
NAMESPACE | 15 ++++++++++
R.Rcheck/00check.log | 3 ++
R/autocorrdiag.R | 25 +++++++++++++++++
R/batchSE.R | 50 +++++++++++++++++++++++++++++++++
R/gelman.R | 7 +++--
R/heidel.R | 6 ++--
R/mcmc.R | 2 +-
R/mcmclist.R | 21 ++++++++------
R/output.R | 78 +++++++++++++++++++++++++++++++++++++---------------
R/rejectionRate.R | 15 ++++++++++
man/as.ts.mcmc.Rd | 7 ++---
man/autocorr.Rd | 10 ++++++-
man/autocorr.diag.Rd | 36 ++++++++++++++++++++++++
man/autocorr.plot.Rd | 2 +-
man/batchSE.Rd | 40 +++++++++++++++++++++++++++
man/effectiveSize.Rd | 5 ++++
man/gelman.diag.Rd | 6 +++-
man/mcmc.convert.Rd | 33 +++++++++++-----------
man/plot.mcmc.Rd | 2 +-
man/rejectionRate.Rd | 31 +++++++++++++++++++++
24 files changed, 370 insertions(+), 66 deletions(-)
diff --git a/AUTHORS b/AUTHORS
index 8403eee..50ce339 100644
--- a/AUTHORS
+++ b/AUTHORS
@@ -35,3 +35,5 @@ contributed to the Statlib archive by Steven Lewis
Andrew Martin <admartin at wustl.edu> added namespace facilities and
fixed a large number of bugs that this change uncovered (see CHANGELOG)
+Russell Almond <almond at acm.org>, <ralmond at ets.org> provided a number
+of cleanups and additions (see CHANGELOG).
diff --git a/CHANGELOG b/CHANGELOG
index 0e429e9..2ebd83a 100644
--- a/CHANGELOG
+++ b/CHANGELOG
@@ -1,3 +1,30 @@
+0.9-4
+
+- Fixed documentation errors.
+
+0.9-3
+
+- Added date stamp check in read.openbugs to ensure files were created
+ at the same time.
+
+0.9-2
+- [RGA] Added an autoburnin flag (default) true to gelman.diag to suppress
+ automatic windowing for burn in (so I can do it manually).
+- [RGA] Fixed a problem where summary.mcmc.list would not give correct
+ pooled standard errors.
+- [RGA] Fixed propagation of standard par and titles to plot.mcmc and
+ plot.mcmc.list. Also fixed so that ask will default to value in
+ par().
+- [RGA] Fixed a problem where autocorr would apply thinning twice to
+ mcmc.list objects.
+- [RGA] Changed effectiveSize for mcmc.list to sum across all chains.
+ Original behavior can be recovered by using lapply(x,effectiveSize)
+- [RGA] Added autocorr.diag function.
+- [RGA] Patched summary.mcmc and summary.mcmc.list so it would give
+ standard error of NA when spectrum0 blows ups.
+- [RGA] Added a rejectionRate method.
+- [RGA] Added a batch Standard Error Function
+
0.9-1
- spectrum0 function now has default max.length argument of 200.
This means that the output will be batched to a length between
diff --git a/DESCRIPTION b/DESCRIPTION
index b6711ef..d8c6270 100644
--- a/DESCRIPTION
+++ b/DESCRIPTION
@@ -1,6 +1,6 @@
Package: coda
-Version: 0.9-1
-Date: 2004/12/20
+Version: 0.9-4
+Date: 2005/4/28
Title: Output analysis and diagnostics for MCMC
Author: Martyn Plummer, Nicky Best, Kate Cowles, Karen Vines
Maintainer: Martyn Plummer <plummer at iarc.fr>
@@ -8,4 +8,4 @@ Depends: R (>= 1.9)
Description: Output analysis and diagnostics for Markov Chain Monte Carlo simulations.
License: GPL Version 2 or later.
URL: http://www-fis.iarc.fr/coda/
-Packaged: Mon Dec 20 12:22:23 2004; martyn
+Packaged: Thu Apr 28 10:45:16 2005; martyn
diff --git a/INDEX b/INDEX
index abf2b96..377df45 100644
--- a/INDEX
+++ b/INDEX
@@ -1,11 +1,12 @@
-[.mcmc Extract or replace parts of MCMC objects
as.matrix.mcmc Conversions of MCMC objects
as.ts.mcmc Coerce mcmc object to time series
autocorr Autocorrelation function for Markov chains
+autocorr.diag Autocorrelation function for Markov chains
autocorr.plot Plot autocorrelations for Markov Chains
+batchSE Batch Standard Error
bugs2jags Convert WinBUGS data file to JAGS data file
-coda.options Options settings for the codamenu driver
codamenu Main menu driver for the coda package
+coda.options Options settings for the codamenu driver
crosscorr Cross correlations for MCMC output
crosscorr.plot Plot image of correlation matrix
cumuplot Cumulative quantile plot
@@ -20,6 +21,7 @@ heidel.diag Heidelberger and Welch's convergence
diagnostic
line Simple linear regression example
mcmc Markov Chain Monte Carlo Objects
+[.mcmc Extract or replace parts of MCMC objects
mcmc.list Replicated Markov Chain Monte Carlo Objects
mcmcUpgrade Upgrade mcmc objects in obsolete format
mcpar Mcpar attribute of MCMC objects
@@ -33,6 +35,7 @@ read.and.check Read data interactively and check that it
read.coda Read output files in CODA format
read.coda.interactive Read CODA output files interactively
read.openbugs Read CODA output files produced by OpenBUGS
+rejectionRate Rejection Rate for Metropolis-Hastings chains
spectrum0 Estimate spectral density at zero
spectrum0.ar Estimate spectral density at zero
summary.mcmc Summary statistics for Markov Chain Monte
diff --git a/NAMESPACE b/NAMESPACE
index 2216469..7dc871d 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -4,7 +4,11 @@ export(
as.mcmc,
as.ts.mcmc,
autocorr,
+ ##RGA
+ autocorr.diag,
autocorr.plot,
+ ##RGA
+ batchSE,
bugs2jags,
chanames,
coda.options,
@@ -37,6 +41,8 @@ export(
read.coda,
read.jags,
read.openbugs,
+ ##RGA
+ rejectionRate,
spectrum0,
spectrum0.ar,
thin,
@@ -71,6 +77,11 @@ S3method(summary, mcmc)
S3method(thin, mcmc)
S3method(time, mcmc)
S3method(window, mcmc)
+##RGA
+S3method(autocorr.diag, mcmc)
+S3method(rejectionRate, mcmc)
+S3method(batchSE,mcmc)
+S3method(as.data.frame,mcmc)
## mcmc.list
S3method("[", mcmc.list)
@@ -83,6 +94,10 @@ S3method(summary, mcmc.list)
S3method(thin, mcmc.list)
S3method(time, mcmc.list)
S3method(window, mcmc.list)
+##RGA
+S3method(autocorr.diag, mcmc.list)
+S3method(rejectionRate, mcmc.list)
+S3method(batchSE,mcmc.list)
## misc
S3method(print, gelman.diag)
diff --git a/R.Rcheck/00check.log b/R.Rcheck/00check.log
new file mode 100644
index 0000000..4bc706b
--- /dev/null
+++ b/R.Rcheck/00check.log
@@ -0,0 +1,3 @@
+* using log directory '/home/martyn/bayes/diagnostics/coda/0.9/coda/R.Rcheck'
+* using R version 2.1.0, 2005-04-18
+* checking for file 'R/DESCRIPTION' ... NO
diff --git a/R/autocorrdiag.R b/R/autocorrdiag.R
new file mode 100644
index 0000000..ac97b6b
--- /dev/null
+++ b/R/autocorrdiag.R
@@ -0,0 +1,25 @@
+### This is a replacement for the autocorrelation function
+### which only gives the autocorrelations and not the cross correlations.
+
+"autocorr.diag" <- function (mcmc.obj, ...) {
+ UseMethod("autocorr.diag")
+}
+
+### Only looks at the diagonal elements of the autocorrlation matrix.
+"autocorr.diag.mcmc" <- function (mcmc.obj,...) {
+ ac <- autocorr(mcmc.obj,...)
+ dd <- dim(ac)
+ result <- matrix(NA,nrow=dd[1],ncol=dd[2],dimnames =dimnames(ac)[1:2])
+ for (i in 1:dd[2])
+ result[,i] <- ac[,i,i]
+ return (result)
+}
+
+### Looks at the average autocorrelation for all chains series.
+"autocorr.diag.mcmc.list" <- function (mcmc.obj,...) {
+ ac <- lapply(mcmc.obj,autocorr.diag.mcmc,...)
+ result <- ac[[1]]
+ for (chain in 2:length(ac))
+ result <- result + ac[[chain]]
+ result/length(ac)
+}
diff --git a/R/batchSE.R b/R/batchSE.R
new file mode 100644
index 0000000..f307b4c
--- /dev/null
+++ b/R/batchSE.R
@@ -0,0 +1,50 @@
+
+"batchSE" <-
+ function(x,batchSize=100) {
+ UseMethod("batchSE")
+}
+
+"batchSE.mcmc" <-
+ function(x,batchSize=100) {
+ niter <- niter(x)
+ nbatch <- niter%/%batchSize
+ ## Truncate the odd lot observations
+ ## Do this off the front instead of the back??
+ niter <- nbatch*batchSize
+ ibatch <- rep(1:nbatch,each=batchSize)[1:niter]
+ batchMeans <- t(sapply(split(data.frame(x[1:niter,]),ibatch),
+ function(batch) apply(batch,2,mean)))
+ grandMean <- apply(batchMeans,2,mean)
+ mi2 <- sweep(batchMeans,2,grandMean,"-")^2
+ stds<-sqrt(apply(mi2,2,sum)*batchSize/(nbatch-1))
+ names(stds) <- dimnames(x)[[2]]
+ stds/sqrt(niter(x))
+}
+
+"batchSE.mcmc.list" <-
+ function(x,batchSize=100) {
+ nchain <- nchain(x)
+ niter <- niter(x)
+ nbatch <- niter%/%batchSize
+ ## truncate odd lot observations
+ niter <- nbatch*batchSize
+ ibatch <- rep(1:nbatch,each=batchSize)[1:niter]
+ batchMeans <- NULL
+ for (i in 1:nchain) {
+ batchMeans <- rbind(batchMeans,
+ t(sapply(split(data.frame(x[[i]][1:niter,]),ibatch),
+ function(batch) apply(batch,2,mean))))
+ }
+ #print(batchMeans)
+ grandMean <- apply(batchMeans,2,mean)
+ #cat("Grand Mean = ",grandMean,"\n")
+ mi2 <- sweep(batchMeans,2,grandMean,"-")^2
+ stds<-sqrt(apply(mi2,2,sum)*batchSize/(nchain*nbatch-1))
+ names(stds) <- dimnames(x[[1]])[[2]]
+ stds/sqrt(niter(x)*nchain(x))
+}
+
+## Needed for this function, but generally useful anyway.
+as.data.frame.mcmc <- function(x, row.names = NULL, optional=FALSE) {
+ as.data.frame.matrix(x,row.names,optional)
+}
diff --git a/R/gelman.R b/R/gelman.R
index fead51e..d3d91d0 100644
--- a/R/gelman.R
+++ b/R/gelman.R
@@ -1,4 +1,5 @@
-"gelman.diag" <- function (x, confidence = 0.95, transform = FALSE)
+"gelman.diag" <- function (x, confidence = 0.95, transform = FALSE,
+ autoburnin=TRUE)
## Gelman and Rubin's diagnostic
## Gelman, A. and Rubin, D (1992). Inference from iterative simulation
## using multiple sequences. Statistical Science, 7, 457-551.
@@ -12,7 +13,9 @@
x <- as.mcmc.list(x)
if (nchain(x) < 2)
stop("You need at least two chains")
- if (start(x) < end(x)/2)
+ ## RGA added an autoburnin parameter here, because if I have already
+ ## trimmed burn in, I don't want to do it again.
+ if (autoburnin && start(x) < end(x)/2 )
x <- window(x, start = end(x)/2 + 1)
Niter <- niter(x)
Nchain <- nchain(x)
diff --git a/R/heidel.R b/R/heidel.R
index 4282f2d..4dcec69 100644
--- a/R/heidel.R
+++ b/R/heidel.R
@@ -94,9 +94,9 @@ effectiveSize <- function(x)
{
if (is.mcmc.list(x))
{
- ans <- vector("list", nchain(x))
- for (i in 1:nchain(x))
- ans[[i]] <- effectiveSize(x[[i]])
+ ##RGA changed to sum across all chains
+ ess <- do.call("rbind",lapply(x,effectiveSize))
+ ans <- apply(ess,2,sum)
}
else
{
diff --git a/R/mcmc.R b/R/mcmc.R
index 026c1de..180be85 100644
--- a/R/mcmc.R
+++ b/R/mcmc.R
@@ -17,7 +17,7 @@
"as.mcmc.default" <- function (x)
if (is.mcmc(x)) x else mcmc(x)
-"as.ts.mcmc" <- function (x)
+"as.ts.mcmc" <- function (x, ...)
{
x <- as.mcmc(x)
y <- ts(x, start = start(x), end = end(x), deltat = thin(x))
diff --git a/R/mcmclist.R b/R/mcmclist.R
index b4d1fcd..c49d52d 100644
--- a/R/mcmclist.R
+++ b/R/mcmclist.R
@@ -72,8 +72,9 @@
"plot.mcmc.list" <-
function (x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
- auto.layout = TRUE, ask = TRUE, ...)
+ auto.layout = TRUE, ask = par("ask"), ...)
{
+## RGA fixed to use default ask value.
oldpar <- NULL
on.exit(par(oldpar))
if (auto.layout) {
@@ -83,17 +84,20 @@
}
for (i in 1:nvar(x)) {
if (trace)
- traceplot(x[, i, drop = FALSE], smooth = smooth)
+ ## RGA fixed to propagate ... argument.
+ traceplot(x[, i, drop = FALSE], smooth = smooth, ...)
if (density) {
if (missing(bwf))
- densplot(x[, i, drop = FALSE])
- else densplot(x[, i, drop = FALSE], bwf = bwf)
+ ## RGA fixed to propagate ... argument.
+ densplot(x[, i, drop = FALSE], ...)
+ else densplot(x[, i, drop = FALSE], bwf = bwf, ...)
}
if (i==1)
oldpar <- c(oldpar, par(ask = ask))
}
}
+
"summary.mcmc.list" <-
function (object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
{
@@ -105,12 +109,12 @@
if (is.matrix(x[[1]])) {
for (i in 1:nchain(x))
for(j in 1:nvar(x))
- xtsvar[i, j] <- spectrum0(x[[i]][,j])$spec
+ xtsvar[i, j] <- safespec0(x[[i]][,j])
xlong <- do.call("rbind", x)
}
else {
for (i in 1:nchain(x))
- xtsvar[i, ] <- spectrum0(x[[i]])$spec
+ xtsvar[i, ] <- safespec0(x[[i]])
xlong <- as.matrix(x)
}
@@ -120,8 +124,9 @@
varquant <- t(apply(xlong, 2, quantile, quantiles))
varstats[, 1] <- xmean
varstats[, 2] <- sqrt(xvar)
- varstats[, 3] <- sqrt(xvar/niter(x))
- varstats[, 4] <- sqrt(xtsvar/niter(x))
+ ##RGA fixed so now give correct std error for pooled (across chains).
+ varstats[, 3] <- sqrt(xvar/(niter(x)*nchain(x)))
+ varstats[, 4] <- sqrt(xtsvar/(niter(x)*nchain(x)))
varquant <- drop(varquant)
varstats <- drop(varstats)
out <- list(statistics = varstats, quantiles = varquant,
diff --git a/R/output.R b/R/output.R
index dea78e6..8099f68 100644
--- a/R/output.R
+++ b/R/output.R
@@ -1,13 +1,15 @@
"autocorr" <-
function (x, lags = c(0, 1, 5, 10, 50), relative = TRUE)
{
+ ## RGA moved MCMC list processing first, else thinning gets
+ ## applied twice. Thanks to Denise Chang for finding this.
+ if (is.mcmc.list(x))
+ return(lapply(x, autocorr, lags, relative))
if (relative)
lags <- lags * thin(x)
else if (any(lags%%thin(x) != 0))
stop("Lags do not conform to thinning interval")
lags <- lags[lags < niter(x) * thin(x)]
- if (is.mcmc.list(x))
- return(lapply(x, autocorr, lags, relative))
x <- as.mcmc(x)
y <- array(dim = c(length(lags), nvar(x), nvar(x)))
dimnames(y) <- list(paste("Lag", lags), varnames(x), varnames(x))
@@ -19,9 +21,10 @@ function (x, lags = c(0, 1, 5, 10, 50), relative = TRUE)
}
"autocorr.plot" <-
-function (x, lag.max, auto.layout = TRUE, ask = TRUE, ...)
+function (x, lag.max, auto.layout = TRUE, ask = par("ask"), ...)
{
- oldpar <- NULL
+## RGA fixed to use default ask value.
+ oldpar <- NULL
on.exit(par(oldpar))
if (auto.layout)
oldpar <- par(mfrow = set.mfrow(Nchains = nchain(x),
@@ -154,25 +157,38 @@ function (x, show.obs = TRUE, bwf, main = "", ylim, ...)
read.coda(file, index.file, start, end, thin, quiet)
}
-"read.openbugs" <- function(stem="", start, end, thin, quiet=FALSE)
+"read.openbugs" <-
+function (stem = "", start, end, thin, quiet = FALSE)
{
+ index.file <- paste(stem, "CODAindex.txt", sep = "")
+ if (!file.exists(index.file))
+ stop("No index file found")
+ index.date <- file.info(index.file)$ctime
+
nchain <- 0
while (TRUE) {
- output.file <- paste(stem,"CODAchain",nchain+1,".txt",sep="")
- if (file.exists(output.file))
+ output.file <- paste(stem, "CODAchain", nchain + 1, ".txt",
+ sep = "")
+ if (file.exists(output.file)) {
nchain <- nchain + 1
- else
- break
+ output.date <- file.info(output.file)$ctime
+ dt <- difftime(index.date, output.date, units="mins")
+ if(abs(as.numeric(dt)) > 1 ) {
+ warning(paste("Files \"",index.file,"\" and \"",output.file,
+ "\" were created at different times\n",sep=""))
+ }
+ }
+ else break
+
}
-
- if (nchain==0)
+ if (nchain == 0)
stop("No output files found")
- index.file <- paste(stem,"CODAindex.txt",sep="")
- ans <- vector("list",nchain)
+ ans <- vector("list", nchain)
for (i in 1:nchain) {
- output.file <- paste(stem,"CODAchain",i,".txt",sep="")
- ans[[i]] <- read.coda(output.file, index.file, start, end, thin, quiet)
+ output.file <- paste(stem, "CODAchain", i, ".txt", sep = "")
+ ans[[i]] <- read.coda(output.file, index.file, start,
+ end, thin, quiet)
}
return(mcmc.list(ans))
}
@@ -304,8 +320,9 @@ function (x, smooth = TRUE, col = 1:6, type = "l", ylab = "",
}
"plot.mcmc" <- function (x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
- auto.layout = TRUE, ask = TRUE, ...)
+ auto.layout = TRUE, ask = par("ask"), ...)
{
+## RGA fixed to use default ask value.
oldpar <- NULL
on.exit(par(oldpar))
if (auto.layout) {
@@ -316,18 +333,34 @@ function (x, smooth = TRUE, col = 1:6, type = "l", ylab = "",
for (i in 1:nvar(x)) {
y <- mcmc(as.matrix(x)[, i, drop=FALSE], start(x), end(x), thin(x))
if (trace)
- traceplot(y, smooth = smooth)
+ ## RGA fixed to propagate ... argument.
+ traceplot(y, smooth = smooth, ...)
if (density) {
if (missing(bwf))
- densplot(y)
+ ## RGA fixed to propagate ... argument.
+ densplot(y, ...)
else
- densplot(y, bwf = bwf)
+ densplot(y, bwf = bwf, ...)
}
if (i==1)
oldpar <- c(oldpar, par(ask=ask))
}
}
+
+### RGA This is a wrapper for spectrum0 which returns NA if
+### spectrum0 crashes. This has happened to me several times when
+### there was bug in my MCMC algorithm.
+"safespec0" <-
+ function (x) {
+ result <- try(spectrum0(x)$spec)
+ ## R
+ if (class(result) == "try-error") result <- NA
+ ## S-Plus
+ if (class(result) == "try") result <- NA
+ result
+}
+
"summary.mcmc" <-
function (object, quantiles = c(0.025, 0.25, 0.5, 0.75, 0.975), ...)
{
@@ -335,17 +368,18 @@ function (x, smooth = TRUE, col = 1:6, type = "l", ylab = "",
statnames <- c("Mean", "SD", "Naive SE", "Time-series SE")
varstats <- matrix(nrow = nvar(x), ncol = length(statnames),
dimnames = list(varnames(x), statnames))
- sp0 <- function(x) spectrum0(x)$spec
+ ## RGA replaced with safespec0
+ #sp0 <- function(x) spectrum0(x)$spec
if (is.matrix(x)) {
xmean <- apply(x, 2, mean)
xvar <- apply(x, 2, var)
- xtsvar <- apply(x, 2, sp0)
+ xtsvar <- apply(x, 2, safespec0)
varquant <- t(apply(x, 2, quantile, quantiles))
}
else {
xmean <- mean(x, na.rm = TRUE)
xvar <- var(x, na.rm = TRUE)
- xtsvar <- sp0(x)
+ xtsvar <- safespec(x)
varquant <- quantile(x, quantiles)
}
varstats[, 1] <- xmean
diff --git a/R/rejectionRate.R b/R/rejectionRate.R
new file mode 100644
index 0000000..8a3a29f
--- /dev/null
+++ b/R/rejectionRate.R
@@ -0,0 +1,15 @@
+
+rejectionRate.mcmc <- function (x) {
+ n <- nrow(x)
+ apply(x[1:(n-1),] == x[2:n,],2,mean)
+}
+
+rejectionRate.mcmc.list <- function (x) {
+ apply(sapply(x,rejectionRate.mcmc),1,mean)
+}
+
+rejectionRate <- function(x) {
+ UseMethod("rejectionRate")
+}
+
+
diff --git a/man/as.ts.mcmc.Rd b/man/as.ts.mcmc.Rd
index d5c60b2..e204f00 100644
--- a/man/as.ts.mcmc.Rd
+++ b/man/as.ts.mcmc.Rd
@@ -2,19 +2,16 @@
\alias{as.ts.mcmc}
\title{Coerce mcmc object to time series}
-\usage{as.ts.mcmc(x)}
+\usage{as.ts.mcmc(x,...)}
\arguments{
\item{x}{an mcmc object}
+ \item{...}{unused arguments for compatibility with generic \code{as.ts}}
}
\description{
\code{as.ts.mcmc} will coerce an mcmc object to a time series.}
-\note{
-This is pretending to be an \code{as.ts} method for mcmc objects. But
-\code{as.ts} is not generic.}
-
\author{Martyn Plummer}
\seealso{
diff --git a/man/autocorr.Rd b/man/autocorr.Rd
index 3f111ac..307b6d6 100644
--- a/man/autocorr.Rd
+++ b/man/autocorr.Rd
@@ -2,7 +2,15 @@
\alias{autocorr}
\title{Autocorrelation function for Markov chains}
-\usage{autocorr(mcmc.obj, lags = c(0, 1, 5, 10, 50), relative=TRUE }
+\usage{autocorr(x, lags = c(0, 1, 5, 10, 50), relative=TRUE)}
+
+\arguments{
+\item{x}{an mcmc object}
+\item{lags}{a vector of lags at which to calculate the autocorrelation}
+\item{relative}{a logical flag. TRUE if lags are relative to the thinning
+interval of the chain, or FALSE if they are absolute difference in iteration
+numbers}
+}
\description{
\code{autocorr} calculates the autocorrelation function for the
diff --git a/man/autocorr.diag.Rd b/man/autocorr.diag.Rd
new file mode 100644
index 0000000..3aec818
--- /dev/null
+++ b/man/autocorr.diag.Rd
@@ -0,0 +1,36 @@
+\name{autocorr.diag}
+\alias{autocorr.diag}
+\alias{autocorr.diag.mcmc}
+\alias{autocorr.diag.mcmc.list}
+\title{Autocorrelation function for Markov chains}
+
+\usage{autocorr.diag(mcmc.obj, ...)}
+\arguments{
+ \item{mcmc.obj}{an object of class \code{mcmc} or \code{mcmc.list}}
+ \item{...}{optional arguments to be passed to \code{autocorr}}
+}
+\description{
+\code{autocorr.diag} calculates the autocorrelation function for the
+Markov chain \code{mcmc.obj} at the lags given by \code{lags}.
+The lag values are taken to be relative to the thinning interval
+if \code{relative=TRUE}. Unlike \code{autocorr}, if \code{mcmc.obj}
+has many parmeters it only computes the autocorrelations with itself and
+not the cross correlations. In cases where \code{autocorr} would
+return a matrix, this function returns the diagonal of the matrix.
+Hence it is more useful for chains with many parameters, but may not
+be as helpful at spotting parameters.
+
+If \code{mcmc.obj} is of class \code{mcmc.list} then the returned
+vector is the average autocorrelation across all chains.
+}
+
+\value{
+A vector containing the autocorrelations.
+}
+
+\author{Russell Almond}
+
+\seealso{
+\code{\link{autocorr}, \link{acf}}, \code{\link{autocorr.plot}}.}
+}
+\keyword{ts}
diff --git a/man/autocorr.plot.Rd b/man/autocorr.plot.Rd
index a4cd5b0..0912f79 100644
--- a/man/autocorr.plot.Rd
+++ b/man/autocorr.plot.Rd
@@ -2,7 +2,7 @@
\alias{autocorr.plot}
\title{Plot autocorrelations for Markov Chains}
-\usage{autocorr.plot(x, lag.max, auto.layout = TRUE, ask=TRUE, ...)}
+\usage{autocorr.plot(x, lag.max, auto.layout = TRUE, ask=par("ask"), ...)}
\arguments{
\item{x}{A Markov Chain}
diff --git a/man/batchSE.Rd b/man/batchSE.Rd
new file mode 100644
index 0000000..a41dbf0
--- /dev/null
+++ b/man/batchSE.Rd
@@ -0,0 +1,40 @@
+\name{batchSE}
+\alias{batchSE}
+\title{Batch Standard Error}
+\description{
+Effective standard deviation of population to produce the correct
+standard errors.
+}
+\usage{
+batchSE(x, batchSize=100)
+}
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object.}
+ \item{batchSize}{Number of observations to include in each batch.}
+}
+\details{
+ Because of the autocorrelation, the usual method of taking
+ \code{var(x)/n} overstates the precision of the estimate. This method
+ works around the problem by looking at the means of batches of the
+ parameter. If the batch size is large enough, the batch means should
+ be approximately uncorrelated and the normal formula for computing the
+ standard error should work.
+
+ The batch standard error procedure is usually thought to be not as
+ accurate as the time series methods used in \code{summary} and
+ \code{effectiveSize}. It is included here for completeness.
+}
+\value{
+A vector giving the standard error for each column of \code{x}.
+}
+\references{
+Roberts, GO (1996) Markov chain concepts related to sampling algorithms,
+in Gilks, WR, Richardson, S and Spiegelhalter, DJ, \emph{Markov Chain
+ Monte Carlo in Practice}, Chapman and Hall, 45-58.
+}
+\seealso{
+ \code{\link{spectrum0.ar}}, \code{\link{effectiveSize}},
+ \code{\link{summary.mcmc}}
+}
+\author{Russell Almond}
+\keyword{ts}
diff --git a/man/effectiveSize.Rd b/man/effectiveSize.Rd
index 56359c2..30ef28d 100644
--- a/man/effectiveSize.Rd
+++ b/man/effectiveSize.Rd
@@ -18,6 +18,11 @@ size. \code{n = N} only when there is no autocorrelation.
Estimation of the effective sample size requires estimating the
spectral density at frequency zero. This is done by the function
\code{spectrum0.ar}
+
+For a \code{mcmc.list} object, the effective sizes are summed across
+chains. To get the size for each chain individually use
+\code{lapply(x,effectiveSize)}.
+
}
\value{
A vector giving the effective sample size for each column of \code{x}.
diff --git a/man/gelman.diag.Rd b/man/gelman.diag.Rd
index e58492f..97c8338 100644
--- a/man/gelman.diag.Rd
+++ b/man/gelman.diag.Rd
@@ -5,7 +5,7 @@
%\alias{print.gelman.diag}
\title{Gelman and Rubin's convergence diagnostic}
-\usage{gelman.diag(x, confidence = 0.95, transform=FALSE)}
+\usage{gelman.diag(x, confidence = 0.95, transform=FALSE, autoburnin=TRUE)}
\arguments{
\item{x}{An \code{mcmc.list} object with more than one chain,
@@ -17,6 +17,10 @@ potential scale reduction factor}
\code{x} should be transformed to improve the normality of the
distribution. If set to TRUE, a log transform or logit transform, as
appropriate, will be applied.}
+\item{autoburnin}{a logical flag indicating whether only the second half
+ of the series should be used in the computation. If set to TRUE
+ (default) and \code{start(x)} is less than \code{end(x)/2} then start
+ of series will be adjusted so that only second half of series is used.}
}
\description{
diff --git a/man/mcmc.convert.Rd b/man/mcmc.convert.Rd
index b671318..4c2cdc0 100644
--- a/man/mcmc.convert.Rd
+++ b/man/mcmc.convert.Rd
@@ -6,7 +6,8 @@
\title{Conversions of MCMC objects}
\usage{
-as.matrix(x, iters = FALSE, chains = FALSE))
+as.matrix.mcmc(x, iters = FALSE)
+as.matrix.mcmc.list(x, iters = FALSE, chains = FALSE)
as.array.mcmc.list(x, drop, ...)
}
@@ -20,21 +21,21 @@ as.array.mcmc.list(x, drop, ...)
}
\description{
- These are methods for the generic functions \code{as.matrix()},
- \code{as.array()} and \code{as.mcmc()}.
-
- \code{as.matrix()} strips the MCMC attributes from an \code{mcmc}
- object and returns a matrix. If \code{iters = TRUE} then a
- column is added with the iteration number. For \code{mcmc.list}
- objects, the rows of multiple chains are concatenated and, if
- \code{chains = TRUE} a column is added with the chain number.
-
- \code{mcmc.list} objects can be coerced to 3-dimensional arrays
- with the \code{as.array()} function.
-
- An \code{mcmc.list} object with a single chain can be coerced
- to an \code{mcmc} object with \code{as.mcmc()}. If the argument
- has multiple chains, this causes an error.
+ These are methods for the generic functions \code{as.matrix()},
+ \code{as.array()} and \code{as.mcmc()}.
+
+ \code{as.matrix()} strips the MCMC attributes from an \code{mcmc}
+ object and returns a matrix. If \code{iters = TRUE} then a
+ column is added with the iteration number. For \code{mcmc.list}
+ objects, the rows of multiple chains are concatenated and, if
+ \code{chains = TRUE} a column is added with the chain number.
+
+ \code{mcmc.list} objects can be coerced to 3-dimensional arrays
+ with the \code{as.array()} function.
+
+ An \code{mcmc.list} object with a single chain can be coerced
+ to an \code{mcmc} object with \code{as.mcmc()}. If the argument
+ has multiple chains, this causes an error.
}
\seealso{
diff --git a/man/plot.mcmc.Rd b/man/plot.mcmc.Rd
index 4a92491..a835b60 100644
--- a/man/plot.mcmc.Rd
+++ b/man/plot.mcmc.Rd
@@ -4,7 +4,7 @@
\usage{
\method{plot}{mcmc}(x, trace = TRUE, density = TRUE, smooth = TRUE, bwf,
- auto.layout = TRUE, ask = TRUE, ...)
+ auto.layout = TRUE, ask = par("ask"), ...)
}
\arguments{
diff --git a/man/rejectionRate.Rd b/man/rejectionRate.Rd
new file mode 100644
index 0000000..3df86a5
--- /dev/null
+++ b/man/rejectionRate.Rd
@@ -0,0 +1,31 @@
+\name{rejectionRate}
+\alias{rejectionRate}
+\alias{rejectionRate.mcmc}
+\alias{rejectionRate.mcmc.list}
+\title{Rejection Rate for Metropolis--Hastings chains}
+
+\usage{rejectionRate(x) }
+\arguments{
+ \item{x}{An \code{mcmc} or \code{mcmc.list} object.}
+}
+\description{
+\code{rejectionRate} calculates the fraction of time that a
+Metropolis--Hastings type chain rejected a proposed move. The rejection
+rate is calculates separately for each variable in the \code{mcmc.obj}
+argument, irregardless of whether the variables were drawn separately or
+in a block. In the latter case, the values returned should be the
+same.
+}
+\details{
+For the purposes of this function, a "rejection" has occurred if the
+value of the time series is the same at two successive time points.
+This test is done naively using \code{==} and may produce problems due
+to rounding error.
+}
+\value{
+A vector containing the rejection rates, one for each variable.
+}
+
+\author{Russell Almond}
+
+\keyword{ts}
--
Alioth's /usr/local/bin/git-commit-notice on /srv/git.debian.org/git/debian-science/packages/r-cran-coda.git
More information about the debian-science-commits
mailing list