## Function to implement consistent Batch Means procedure
## (Jones, Haran, Caffo and Neath, 2006, JASA)
## Galin L. Jones, Murali Haran, Brian S. Caffo, and Ronald Neath,
##"Fixed-Width Output Analysis for Markov Chain Monte Carlo"
## input: vals, a vector of N values (from a Markov chain),
## bs=batch size and g, a function 
## output: estimate of E(g(x)) and an estimate of the Monte Carlo
## standard error of estimate of E(g(x))
id <- function(x) return(x)  # default: identity function

bm <- function(vals,bs="sqroot",g=id,warn=FALSE)
  {
    N <- length(vals)
    if (N<1000)
      {
        if (warn) # if warning
          cat("WARNING: too few samples (less than 1000)\n")
        if (N<10)
          return(NA)
      }

    if (bs=="sqroot") 
      {
        b <- floor(sqrt(N)) # batch size
        a <- floor(N/b) # number of batches
      }
    else
      if (bs=="cuberoot") 
        {
          b <- floor(N^(1/3)) # batch size
          a <- floor(N/b) # number of batches
        }
    else # batch size provided
      {
        stopifnot(is.numeric(bs))  
        b <- floor(bs) # batch size
        if (b > 1) # batch size valid
          a <- floor(N/b) # number of batches
        else
          stop("batch size invalid (bs=",bs,")")
      }
    
    Ys <- sapply(1:a,function(k) return(mean(g(vals[((k-1)*b+1):(k*b)]))) )

    muhat <- mean(g(Ys))
    sigmahatsq <- b*sum((Ys-muhat)^2)/(a-1)

    bmse <- sqrt(sigmahatsq/N)

    return(list(est=muhat,se=bmse,bs=bs))
  }