R/computeacf.R

In hadron: Analysis Framework for Monte Carlo Simulation Data in Physics

#' Computes The ACF and Integrated AC Time
#' 
#' Computes the ACF and integrated autocorrelation time of a time series. It
#' also estimates the corresponding standard errors.
#' 
#' The standard error of the ACF is computed using equation (E.11) of M.
#' Luescher, hep-lat/0409106. The error of the integrated autocorrelation time
#' using the Madras Sokal formula, see also hep-lat/0409106.
#' 
#' @param tseries the time series.
#' @param W.max maximal time lag to be used.
#' @param Lambda cut-off needed to estimate the standard error of the ACF.
#' @return It returns a list of class \code{hadronacf} with members
#' \item{lags}{time lags of the integrated autocorrelation function}
#' \item{Gamma}{normalised autocorrelation function} \item{dGamma}{error of
#' normalised autocorrelation function} \item{W.max}{max time lag used for
#' the call of \code{acf}} \item{W}{the cut-off up to which the ACF is
#' integrated for the integrated autocorrelation time} \item{tdata}{the
#' original time series} \item{tau}{the estimated integrated autocorrelation
#' time} \item{dtau}{the estimated error of the integrated autocorrelation
#' time}
#' @author Carsten Urbach, \email{curbach@@gmx.de}
#' @seealso \code{\link{uwerr}}, \code{\link{acf}}
#' \code{\link{bootstrap.analysis}}
#' @references 'Monte Carlo errors with less errors', Ulli Wolff,
#' http://arxiv.org/abs/hep-lat/0306017hep-lat/0306017
#' 
#' 'Schwarz-preconditioned HMC algorithm for two-flavour lattice QCD', Martin
#' Luescher, http://arxiv.org/abs/hep-lat/0409106hep-lat/0409106
#' 
#' N. Madras, A. D. Sokal, J. Stat. Phys. 50 (1988) 109
#' @keywords autocorrelationfunction autocorrelationtime timeseries
#' @examples
#' 
#' data(plaq.sample)
#' myacf <- computeacf(plaq.sample, 300)
#' plot(myacf)
#' summary(myacf)
#' 
#' @export computeacf
computeacf <- function(tseries, W.max, Lambda=100) {
  N <- length(tseries)
  if(missing(W.max)) W.max <- floor(sqrt(N))
  Gamma.tmp <- rep(0, times=2*W.max+W.max+1)
  dGamma <- rep(0., times=W.max+1)
  ## for t > W.max we set Gamma to 0 for simplicity
  Gamma.tmp[1:(W.max+1)] <- acf(tseries, lag.max = W.max, plot=FALSE)$acf

  ## now we determine the error using (E.11) from Luescher, hep-lat/0409106
  for(t in 0:W.max) {
    k <- c(max(1,(t-Lambda)):(t+Lambda))
    dGamma[t+1] <- sum((Gamma.tmp[(k+t+1)]+Gamma.tmp[(abs(k-t)+1)]-2*Gamma.tmp[t+1]*Gamma.tmp[(k+1)])^2);
    dGamma[t+1] <- sqrt(dGamma[t+1]/N)
  }
  Gamma <- Gamma.tmp[1:(W.max+1)]

  ## and we determine where to stop summing Gamma -> W
  ## for a more sophisticated approach see Wolff, hep-lat/0306017
  W <- 0
  for(t in 0:W.max) {
    W <- t
    if(Gamma[t+1] < dGamma[t+1]) break
  }

  ## compute integrated autocorrelation time
  tau <- 0.5 + sum(Gamma[2:(W+1)])
  ## Madras, Sokal approximation for the error, (E.14) in Luescher, hep-lat/0409106
  dtau <- sqrt((4*W + 2) * tau^2 / N)
  
  res <- list(lags = c(0:(W.max)), Gamma=Gamma, dGamma=dGamma,
              W.max=W.max, W=W, tseries=tseries, tau=tau, dtau=dtau)
  attr(res, "class") <- c("hadronacf", "list")
  return(invisible(res))
}
#' plot.hadronacf
#'
#' @description
#' generic function to plot an object of class "myGamma"
#'
#' @param x Object of type `hadronacf` generated by \link{computeacf}
#' @param ... Generic graphical parameters to be passed on
#' @param col String. Color to be used for the data points.
#'
#' @return
#' No return value.
#' 
#' @export
plot.hadronacf <- function (x, ..., col = "black") {
  Gamma <- x
  ## this is to avoid a warning
  Gamma$dGamma[1] <- 0.001
  ## determine ylimits from data
  ylim <- c(min(Gamma$Gamma - 2 * Gamma$dGamma, na.rm = TRUE), 
            max(Gamma$Gamma + 2 * Gamma$dGamma, na.rm = TRUE))
  ## data points
  plot(Gamma$lags, Gamma$Gamma, ylim = ylim,
       col = col, xlab="t", ylab="Gamma[t]", ...)
  ## errors
  arrows(Gamma$lags, Gamma$Gamma - Gamma$dGamma, Gamma$lags, Gamma$Gamma + Gamma$dGamma,
         length = 0.01, angle = 90, code = 3, col = col)
  # cuttoff and baseline
  abline(h=0, col="red")
  abline(v=Gamma$W, col="red")
}

#' summary.hadronacf
#'
#' @description
#' generic function to summarise an object of class "myGamma"
#'
#' @param object Object of type `hadronacf` generated by \link{computeacf}
#' @param ... Generic parameters to be passed on
#'
#' @return
#' No return value.
#' 
#' @export
summary.hadronacf <- function (object, ...) {
  Gamma <- object
  cat("Analysis based on Autocorrelation function\n")
  cat("cut-off parameter W:\t", Gamma$W, "\n")
  cat("tauint:\t\t\t", Gamma$tau, "\n")
  cat("dtauint:\t\t", Gamma$dtau, "\n")
  cat("data mean:\t\t", mean(Gamma$tseries), "\n")
  cat("data error (naive):\t", sqrt(var(Gamma$tseries)/length(Gamma$tseries)), "\n")
  cat("data error (corrected):\t", sqrt(2*Gamma$tau*var(Gamma$tseries)/length(Gamma$tseries)), "\n")
}