Running 'onechain.R' [165s]
  Running 'testchoosepsrf.R' [1s]
  Comparing 'testchoosepsrf.Rout' to 'testchoosepsrf.Rout.save' ... OK
  Running 'testgelman.bm.R' [14s]
  Comparing 'testgelman.bm.Rout' to 'testgelman.bm.Rout.save' ... OK
  Running 'testmaxeigen.R' [182s]
  Running 'testmpsrf.R' [185s]
  Running 'testneff.R' [14s]
  Comparing 'testneff.Rout' to 'testneff.Rout.save' ... OK
Running the tests in 'tests/onechain.R' failed.
Complete output:
  > set.seed(10)
  > library(mvtnorm)
  > library(mcmcse)
  > library(stableGR)
  > 
  > ################ 
  > # Start by making a few chains 
  > #
  > # Details on the chain construction
  > p <- 5
  > N <- 50^2
  > tail.ind <- floor(N*.80):N
  > foo <- matrix(.50, nrow=p, ncol=p)
  > sigma <- foo^(abs(col(foo)-row(foo)))
  > mu <- sample(10:20, p)
  > mu2 <- mu[p]
  > 
  > # Create the chain
  > mvn.gibbs <- stableGR:::mvn.gibbs
  > out.gibbs1 <- mvn.gibbs(N = N, mu = mu, sigma = sigma, p = p)
  > obj <- list(out.gibbs1)
  > 
  > ##############################################
  > # size = NULL, for a single chain (ideal batch size)
  > 
  > # Calculate psrfs with stable.GR
  > size <- NULL
  > method <- "lug"
  > results1 <- stable.GR(obj, blather = TRUE, size = size, method = method)
  > 
  > # Isolate the first component of the  chain
  > chain1 <- matrix(out.gibbs1[ ,1], ncol = 1)
  > GR1chain <- stable.GR(list(chain1), method = method, multivariate = FALSE, size = size, autoburnin=FALSE, blather = TRUE)
  > 
  > # Calculate psrfs by for this chain by hand
  > # Calculate tau^2 for each chain
  > tausq <- (mcse(chain1, method = method)$se)^2 * N 
  > tau2 <- asym.var(list(chain1), multivariate = FALSE, method = method, size = size, autoburnin = FALSE)
  > all.equal(as.numeric(tausq), as.numeric(tau2))
  [1] TRUE
  > all.equal(as.numeric(tausq), as.numeric(GR1chain$blather$asymVars))
  [1] TRUE
  > 
  > # Calculate s^2 
  > ssquared <- var(chain1)
  > 
  > 
  > # Calculate sigma^2 estimate
  > Nneeded <- GR1chain$blather$Nneeded
  > sigsq <- ((Nneeded-1) * ssquared + tausq)/Nneeded
  > othersigsq <- as.numeric(GR1chain$blather$sigmasq)
  > all.equal(as.numeric(sigsq), othersigsq)
  [1] TRUE
  > 
  > # Calculate the diagnostic
  > Rhat <- sigsq / ssquared
  > that <- sqrt(Rhat)
  > all.equal(as.numeric(that[1]) , as.numeric(GR1chain$psrf))
  [1] TRUE
  > 
  > 
  > ######################################## 
  > # Calculate psrfs with stable.GR
  > size <- "sqroot"
  > method <- "lug"
  > results2 <- stable.GR(obj, blather = TRUE, size = size, method = method)
  > results3 <- stable.GR(chain1, blather = TRUE, size = size, method = method)
  > all.equal(as.numeric(results2$psrf[1]), as.numeric(results3$psrf[1]))
  [1] TRUE
  > 
  > 
  > # Calculate psrfs by for this chain by hand
  > # Calculate tau^2 for each chain
  > tausq <- (mcse(chain1, method = method, size = size)$se)^2 * N 
  > tau2 <- asym.var(list(chain1), multivariate = FALSE, method = method, size = size, autoburnin = FALSE)
  > all.equal(as.numeric(tausq), as.numeric(tau2))
  [1] TRUE
  > all.equal(as.numeric(tausq), as.numeric(results2$blather$asymVars[1]))
  [1] TRUE
  > all.equal(as.numeric(tausq), as.numeric(results3$blather$asymVars))
  [1] TRUE
  > 
  > # Calculate s^2 
  > ssquared <- var(chain1)
  > 
  > 
  > # Calculate sigma^2 estimate
  > Nneeded <- results2$blather$Nneeded
  > sigsq <- ((Nneeded-1) * ssquared + tausq)/Nneeded
  > othersigsq <- as.numeric(results3$blather$sigmasq)
  > all.equal(as.numeric(sigsq), othersigsq)
  [1] TRUE
  > othersigsq <- as.numeric(results2$blather$sigmasq[1])
  > all.equal(as.numeric(sigsq), othersigsq)
  [1] TRUE
  > 
  > 
  > # Calculate the diagnostic
  > Rhat <- sigsq / ssquared
  > that <- sqrt(Rhat)
  > all.equal(as.numeric(that[1]) , as.numeric(results2$psrf[1]))
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.17    0.12    2.28 
Running the tests in 'tests/testmaxeigen.R' failed.
Complete output:
  > set.seed(10)
  > library(mvtnorm)
  > library(mcmcse)
  > library(stableGR)
  > 
  > ################ 
  > # Start by making a few chains to work with
  > 
  > # Details on the chain construction
  > p <- 5
  > N <- 10000
  > tail.ind <- floor(N*.80):N
  > foo <- matrix(.50, nrow=p, ncol=p)
  > sigma <- foo^(abs(col(foo)-row(foo)))
  > mu <- sample(10:20, p)
  > mu2 <- mu[p]
  > 
  > # Create the chains
  > mvn.gibbs <- stableGR:::mvn.gibbs
  > out.gibbs1 <- mvn.gibbs(N = N, mu = mu, sigma = sigma, p = p)
  > out.gibbs2 <- mvn.gibbs(N = N, mu = mu, sigma = sigma, p = p)
  > 
  > obj <- list(out.gibbs1, out.gibbs2)
  > 
  > ################ 
  > # Perform unit test using the two chains in obj
  > 
  > withfun <- stable.GR(obj, mapping = "maxeigen", size = "sqroot")$mpsrf
  > 
  > 
  > # Calculate Tmat for each chain
  > stacked <- rbind(out.gibbs1, out.gibbs2)
  > That <- mcse.multi(stacked, method = "lug", size = sqrt(N))$cov
  > #Tmat1 <- mcse.multi(out.gibbs1, method = "lug")$cov
  > #Tmat2 <- mcse.multi(out.gibbs2, method = "lug")$cov
  > #That <- .5*(Tmat1 + Tmat2) #good
  > 
  > #Calc Smat
  > cov1 <- var(out.gibbs1)
  > cov2 <- var(out.gibbs2)
  > Smat <- .5*(cov1 + cov2) #good
  > 
  > Sinv <- qr.solve(Smat)
  > thingy <- Sinv %*% That
  > maxeigen <- max(eigen(thingy)$values)
  > 
  > Nchain <- length(obj)
  > all.equal(2, Nchain)
  [1] TRUE
  > 
  > rhat <- (N-1)/N +  maxeigen/N
  > byhand <- sqrt(rhat)
  > all.equal(byhand, withfun)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
    13.23    0.43   13.67 
Running the tests in 'tests/testmpsrf.R' failed.
Complete output:
  > set.seed(10)
  > library(mvtnorm)
  > library(mcmcse)
  > library(stableGR)
  > 
  > ################ 
  > # Start by making a few chains to work with
  > 
  > # Details on the chain construction
  > p <- 5
  > N <- 10000
  > tail.ind <- floor(N*.80):N
  > foo <- matrix(.50, nrow=p, ncol=p)
  > sigma <- foo^(abs(col(foo)-row(foo)))
  > mu <- sample(10:20, p)
  > mu2 <- mu[p]
  > 
  > # Create the chains
  > mvn.gibbs <- stableGR:::mvn.gibbs
  > out.gibbs1 <- mvn.gibbs(N = N, mu = mu, sigma = sigma, p = p)
  > out.gibbs2 <- mvn.gibbs(N = N, mu = mu, sigma = sigma, p = p)
  > 
  > obj <- list(out.gibbs1, out.gibbs2)
  > 
  > ################ 
  > # Perform unit test using the two chains in obj
  > 
  > outwithfun <- stable.GR(obj, blather = TRUE, size = "sqroot")
  > withfun <- outwithfun$mpsrf
  > blather <- outwithfun$blather
  > 
  > # Calculate Tmat for each chain
  > stacked <- rbind(out.gibbs1, out.gibbs2)
  > That <- mcse.multi(stacked, method = "lug", size = sqrt(N))$cov
  > Tmat1 <- mcse.multi(out.gibbs1, method = "lug", size = "sqroot")$cov
  > 
  > all.equal(That, blather$AsymVarMatrix)
  [1] TRUE
  > 
  > #Calc Smat
  > sloan <- stableGR:::size.and.trim(obj, size = "sqroot")
  > trimmedchains <- sloan$trimmedchains
  > stackedchains <- do.call(rbind, trimmedchains)
  > Smat <- var(stackedchains)
  > all.equal(Smat, blather$S)
  [1] TRUE
  > 
  > 
  > detratio <- det(That)/det(Smat) 
  > all.equal(detratio, det(solve(blather$S, blather$AsymVarMatrix)))
  [1] TRUE
  > 
  > Nchain <- length(obj)
  > all.equal(2, Nchain)
  [1] TRUE
  > 
  > top <- (N-1) + ((detratio)^(1/p))
  > byhand <- sqrt(top/N)
  > 
  > all.equal(byhand, withfun)
  [1] TRUE
  > 
  > ################ ################ ################
  > # Perform unit test using a SINGLE chain (just in case)
  > onechain <- list(out.gibbs1)
  > withfun <- stable.GR(onechain, size = "sqroot")$mpsrf
  > 
  > #calculate determinants
  > Teigen <- eigen(Tmat1)$values
  > cov1 <- var(out.gibbs1)
  > Seigen <- eigen(cov1)$values
  > detT <- (prod(Teigen))
  > detS <- (prod(Seigen))
  > detratio <- detT/detS #good
  > 
  > Nchain <- length(onechain)
  > all.equal(1, Nchain)
  [1] TRUE
  > 
  > rhat <- (N-1)/N + ((detratio)^(1/p))/N
  > byhand <- sqrt(rhat)
  > all.equal(byhand, withfun)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
    13.59    0.28   13.92