Running 'accept-batch.R' [0s]
  Comparing 'accept-batch.Rout' to 'accept-batch.Rout.save' ... OK
  Running 'circle.R' [154s]
  Running 'initseq.R' [1s]
  Comparing 'initseq.Rout' to 'initseq.Rout.save' ... OK
  Running 'isotropic.R' [0s]
  Comparing 'isotropic.Rout' to 'isotropic.Rout.save' ... OK
  Running 'logit.R' [1s]
  Comparing 'logit.Rout' to 'logit.Rout.save' ... OK
  Running 'logitbat.R' [156s]
  Running 'logitfun.R' [1s]
  Comparing 'logitfun.Rout' to 'logitfun.Rout.save' ... OK
  Running 'logitfunarg.R' [1s]
  Comparing 'logitfunarg.Rout' to 'logitfunarg.Rout.save' ... OK
  Running 'logitidx.R' [1s]
  Comparing 'logitidx.Rout' to 'logitidx.Rout.save' ... OK
  Running 'logitlogidx.R' [1s]
  Comparing 'logitlogidx.Rout' to 'logitlogidx.Rout.save' ... OK
  Running 'logitmat.R' [1s]
  Comparing 'logitmat.Rout' to 'logitmat.Rout.save' ... OK
  Running 'logitnegidx.R' [1s]
  Comparing 'logitnegidx.Rout' to 'logitnegidx.Rout.save' ... OK
  Running 'logitsub.R' [141s]
  Running 'logitsubbat.R' [156s]
  Running 'logitvec.R' [1s]
  Comparing 'logitvec.Rout' to 'logitvec.Rout.save' ... OK
  Running 'morph.R' [154s]
  Running 'morph.metrop.R' [169s]
  Running 'morphtoo.R' [2s]
  Comparing 'morphtoo.Rout' to 'morphtoo.Rout.save' ... OK
  Running 'saveseed.R' [1s]
  Comparing 'saveseed.Rout' to 'saveseed.Rout.save' ... OK
  Running 'saveseedmorph.R' [172s]
  Running 'temp-par-witch.R' [2s]
  Comparing 'temp-par-witch.Rout' to 'temp-par-witch.Rout.save' ... OK
  Running 'temp-par.R' [174s]
  Running 'temp-ser-witch.R' [156s]
  Running 'temp-ser.R' [3s]
  Comparing 'temp-ser.Rout' to 'temp-ser.Rout.save' ... OK
  Running 'zero-error.R' [0s]
  Comparing 'zero-error.Rout' to 'zero-error.Rout.save' ... OK
Running the tests in 'tests/circle.R' failed.
Complete output:
  > 
  >  epsilon <- 1e-15
  > 
  >  library(mcmc)
  > 
  >  suppressWarnings(RNGkind("Marsaglia-Multicarry"))
  >  set.seed(42)
  > 
  >  d <- 5
  > 
  >  logh <- function(x) {
  +      if (! is.numeric(x)) stop("x not numeric")
  +      if (length(x) != d) stop("length(x) != d")
  +      fred <- 1 - sum(x^2)
  +      if (fred > 0) return(log(fred)) else return(-Inf)
  +  }
  > 
  >  out.metro <- metrop(logh, rep(0, d), 1e3, scale = 0.01)
  >  out.metro$accept
  [1] 0.979
  > 
  >  out.metro <- metrop(out.metro, scale = 0.1)
  >  out.metro$accept
  [1] 0.72
  > 
  >  out.metro <- metrop(out.metro, scale = 0.5)
  >  out.metro$accept
  [1] 0.16
  > 
  >  out.metro <- metrop(out.metro, scale = 0.4)
  >  out.metro$accept
  [1] 0.228
  > 
  >  out.metro <- metrop(out.metro, nbatch = 1e2, debug = TRUE)
  > 
  >  all(out.metro$batch[- out.metro$nbatch, ] == out.metro$current[- 1, ])
  [1] TRUE
  >  all(out.metro$current[1, ] == out.metro$initial)
  [1] TRUE
  >  all(out.metro$batch[out.metro$nbatch, ] == out.metro$final)
  [1] TRUE
  > 
  >  .Random.seed <- out.metro$initial.seed
  >  d <- ncol(out.metro$proposal)
  >  n <- nrow(out.metro$proposal)
  >  my.proposal <- matrix(NA, n, d)
  >  my.u <- double(n)
  >  ska <- out.metro$scale
  >  for (i in 1:n) {
  +      my.proposal[i, ] <- out.metro$current[i, ] + ska * rnorm(d)
  +      if (is.na(out.metro$u[i])) {
  +          my.u[i] <- NA
  +      } else {
  +          my.u[i] <- runif(1)
  +      }
  +  }
  >  max(abs(out.metro$proposal - my.proposal)) < epsilon
  [1] TRUE
  >  all(is.na(out.metro$u) == is.na(my.u))
  [1] TRUE
  >  all(out.metro$u[!is.na(out.metro$u)] == my.u[!is.na(my.u)])
  [1] TRUE
  > 
  >  my.curr.log.green <- apply(out.metro$current, 1, logh)
  >  my.prop.log.green <- apply(out.metro$proposal, 1, logh)
  >  all(is.na(out.metro$u) == ((my.prop.log.green == -Inf) |
  +      (my.prop.log.green > my.curr.log.green)))
  [1] TRUE
  >  foo <- my.prop.log.green - my.curr.log.green
  >  blurfle <- foo - out.metro$log.green
  >  blurfle[foo == -Inf & out.metro$log.green == -Inf] <- 0
  >  max(blurfle) < epsilon
  [1] TRUE
  > 
  >  my.accept <- (my.prop.log.green > -Inf) & (is.na(my.u) | my.u < exp(foo))
  >  sum(my.accept) == round(n * out.metro$accept)
  [1] TRUE
  > 
  >  my.path <- matrix(NA, n, d)
  >  my.path[my.accept, ] <- out.metro$proposal[my.accept, ]
  >  my.path[! my.accept, ] <- out.metro$current[! my.accept, ]
  > 
  >  all(my.path == out.metro$batch)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.28    0.10    0.32 
Running the tests in 'tests/logitbat.R' failed.
Complete output:
  > 
  >  # test batching (blen)
  > 
  >  epsilon <- 1e-15
  > 
  >  library(mcmc)
  > 
  >  suppressWarnings(RNGkind("Marsaglia-Multicarry"))
  >  set.seed(42)
  > 
  >  n <- 100
  >  rho <- 0.5
  >  beta0 <- 0.25
  >  beta1 <- 1
  >  beta2 <- 0.5
  > 
  >  x1 <- rnorm(n)
  >  x2 <- rho * x1 + sqrt(1 - rho^2) * rnorm(n)
  >  eta <- beta0 + beta1 * x1 + beta2 * x2
  >  p <- 1 / (1 + exp(- eta))
  >  y <- as.numeric(runif(n) < p)
  > 
  >  out <- glm(y ~ x1 + x2, family = binomial())
  > 
  >  logl <- function(beta) {
  +      if (length(beta) != 3) stop("length(beta) != 3")
  +      beta0 <- beta[1]
  +      beta1 <- beta[2]
  +      beta2 <- beta[3]
  +      eta <- beta0 + beta1 * x1 + beta2 * x2
  +      p <- exp(eta) / (1 + exp(eta))
  +      return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0])))
  +  }
  > 
  >  out.metro <- metrop(logl, coefficients(out), 1e3, scale = 0.01)
  >  out.metro$accept
  [1] 0.982
  > 
  >  out.metro <- metrop(out.metro, scale = 0.1)
  >  out.metro$accept
  [1] 0.795
  > 
  >  out.metro <- metrop(out.metro, scale = 0.5)
  >  out.metro$accept
  [1] 0.264
  > 
  >  apply(out.metro$batch, 2, mean)
  [1] 0.06080257 1.42304941 0.52634149
  > 
  >  out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
  +      scale = 0.5, debug = TRUE, blen = 5)
  > 
  >  niter <- out.metro$nbatch * out.metro$blen * out.metro$nspac
  >  niter == nrow(out.metro$current)
  [1] TRUE
  >  niter == nrow(out.metro$proposal)
  [1] TRUE
  >  all(out.metro$current[1, ] == out.metro$initial)
  [1] TRUE
  >  all(out.metro$current[niter, ] == out.metro$final) |
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  [1] TRUE
  > 
  >  .Random.seed <- out.metro$initial.seed
  >  d <- ncol(out.metro$proposal)
  >  n <- nrow(out.metro$proposal)
  >  my.proposal <- matrix(NA, n, d)
  >  my.u <- double(n)
  >  ska <- out.metro$scale
  >  for (i in 1:n) {
  +      my.proposal[i, ] <- out.metro$current[i, ] + ska * rnorm(d)
  +      if (is.na(out.metro$u[i])) {
  +          my.u[i] <- NA
  +      } else {
  +          my.u[i] <- runif(1)
  +      }
  +  }
  >  max(abs(out.metro$proposal - my.proposal)) < epsilon
  [1] TRUE
  >  all(is.na(out.metro$u) == is.na(my.u))
  [1] TRUE
  >  all(out.metro$u[!is.na(out.metro$u)] == my.u[!is.na(my.u)])
  [1] TRUE
  > 
  >  my.curr.log.green <- apply(out.metro$current, 1, logl)
  >  my.prop.log.green <- apply(out.metro$proposal, 1, logl)
  >  all(is.na(out.metro$u) == (my.prop.log.green > my.curr.log.green))
  [1] TRUE
  >  foo <- my.prop.log.green - my.curr.log.green
  >  max(abs(foo - out.metro$log.green)) < epsilon
  [1] TRUE
  > 
  >  my.accept <- is.na(my.u) | my.u < exp(foo)
  >  sum(my.accept) == round(n * out.metro$accept)
  [1] TRUE
  >  if (my.accept[niter]) {
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  +  } else {
  +      all(out.metro$current[niter, ] == out.metro$final)
  +  }
  [1] TRUE
  > 
  >  my.current <- out.metro$current
  >  my.current[my.accept, ] <- my.proposal[my.accept, ]
  >  my.current <- rbind(out.metro$initial, my.current[- niter, ])
  >  max(abs(out.metro$current - my.current)) < epsilon
  [1] TRUE
  > 
  >  my.path <- matrix(NA, n, d)
  >  my.path[my.accept, ] <- out.metro$proposal[my.accept, ]
  >  my.path[! my.accept, ] <- out.metro$current[! my.accept, ]
  >  nspac <- out.metro$nspac
  > 
  >  my.path <- my.path[seq(nspac, niter, by = nspac), ]
  > 
  >  foom <- array(as.vector(t(my.path)), c(d, out.metro$blen, out.metro$nbatch))
  >  boom <- t(apply(foom, c(1, 3), mean))
  > 
  >  all(dim(boom) == dim(out.metro$batch))
  [1] TRUE
  >  max(abs(boom - out.metro$batch)) < epsilon
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.54    0.10    0.62 
Running the tests in 'tests/logitsub.R' failed.
Complete output:
  > 
  >  # test spacing (nspac)
  > 
  >  epsilon <- 1e-15
  > 
  >  library(mcmc)
  > 
  >  suppressWarnings(RNGkind("Marsaglia-Multicarry"))
  >  set.seed(42)
  > 
  >  n <- 100
  >  rho <- 0.5
  >  beta0 <- 0.25
  >  beta1 <- 1
  >  beta2 <- 0.5
  > 
  >  x1 <- rnorm(n)
  >  x2 <- rho * x1 + sqrt(1 - rho^2) * rnorm(n)
  >  eta <- beta0 + beta1 * x1 + beta2 * x2
  >  p <- 1 / (1 + exp(- eta))
  >  y <- as.numeric(runif(n) < p)
  > 
  >  out <- glm(y ~ x1 + x2, family = binomial())
  > 
  >  logl <- function(beta) {
  +      if (length(beta) != 3) stop("length(beta) != 3")
  +      beta0 <- beta[1]
  +      beta1 <- beta[2]
  +      beta2 <- beta[3]
  +      eta <- beta0 + beta1 * x1 + beta2 * x2
  +      p <- exp(eta) / (1 + exp(eta))
  +      return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0])))
  +  }
  > 
  >  out.metro <- metrop(logl, coefficients(out), 1e3, scale = 0.01)
  >  out.metro$accept
  [1] 0.982
  > 
  >  out.metro <- metrop(out.metro, scale = 0.1)
  >  out.metro$accept
  [1] 0.795
  > 
  >  out.metro <- metrop(out.metro, scale = 0.5)
  >  out.metro$accept
  [1] 0.264
  > 
  >  apply(out.metro$batch, 2, mean)
  [1] 0.06080257 1.42304941 0.52634149
  > 
  >  out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
  +      scale = 0.5, debug = TRUE, nspac = 3)
  > 
  >  niter <- out.metro$nbatch * out.metro$blen * out.metro$nspac
  >  niter == nrow(out.metro$current)
  [1] TRUE
  >  niter == nrow(out.metro$proposal)
  [1] TRUE
  >  all(out.metro$current[1, ] == out.metro$initial)
  [1] TRUE
  >  all(out.metro$current[niter, ] == out.metro$final) |
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  [1] TRUE
  > 
  >  .Random.seed <- out.metro$initial.seed
  >  d <- ncol(out.metro$proposal)
  >  n <- nrow(out.metro$proposal)
  >  my.proposal <- matrix(NA, n, d)
  >  my.u <- double(n)
  >  ska <- out.metro$scale
  >  for (i in 1:n) {
  +      my.proposal[i, ] <- out.metro$current[i, ] + ska * rnorm(d)
  +      if (is.na(out.metro$u[i])) {
  +          my.u[i] <- NA
  +      } else {
  +          my.u[i] <- runif(1)
  +      }
  +  }
  >  max(abs(out.metro$proposal - my.proposal)) < epsilon
  [1] TRUE
  >  all(is.na(out.metro$u) == is.na(my.u))
  [1] TRUE
  >  all(out.metro$u[!is.na(out.metro$u)] == my.u[!is.na(my.u)])
  [1] TRUE
  > 
  >  my.curr.log.green <- apply(out.metro$current, 1, logl)
  >  my.prop.log.green <- apply(out.metro$proposal, 1, logl)
  >  all(is.na(out.metro$u) == (my.prop.log.green > my.curr.log.green))
  [1] TRUE
  >  foo <- my.prop.log.green - my.curr.log.green
  >  max(abs(foo - out.metro$log.green)) < epsilon
  [1] TRUE
  > 
  >  my.accept <- is.na(my.u) | my.u < exp(foo)
  >  sum(my.accept) == round(n * out.metro$accept)
  [1] TRUE
  >  if (my.accept[niter]) {
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  +  } else {
  +      all(out.metro$current[niter, ] == out.metro$final)
  +  }
  [1] TRUE
  > 
  >  my.current <- out.metro$current
  >  my.current[my.accept, ] <- my.proposal[my.accept, ]
  >  my.current <- rbind(out.metro$initial, my.current[- niter, ])
  >  max(abs(out.metro$current - my.current)) < epsilon
  [1] TRUE
  > 
  >  my.path <- matrix(NA, n, d)
  >  my.path[my.accept, ] <- out.metro$proposal[my.accept, ]
  >  my.path[! my.accept, ] <- out.metro$current[! my.accept, ]
  >  nspac <- out.metro$nspac
  >  my.path <- my.path[seq(nspac, niter, by = nspac), ]
  >  all(dim(my.path) == dim(out.metro$batch))
  [1] TRUE
  > 
  >  all(my.path == out.metro$batch)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.42    0.15    0.53 
Running the tests in 'tests/logitsubbat.R' failed.
Complete output:
  > 
  >  # test batching (blen) and spacing (nspac) together
  > 
  >  epsilon <- 1e-15
  > 
  >  library(mcmc)
  > 
  >  suppressWarnings(RNGkind("Marsaglia-Multicarry"))
  >  set.seed(42)
  > 
  >  n <- 100
  >  rho <- 0.5
  >  beta0 <- 0.25
  >  beta1 <- 1
  >  beta2 <- 0.5
  > 
  >  x1 <- rnorm(n)
  >  x2 <- rho * x1 + sqrt(1 - rho^2) * rnorm(n)
  >  eta <- beta0 + beta1 * x1 + beta2 * x2
  >  p <- 1 / (1 + exp(- eta))
  >  y <- as.numeric(runif(n) < p)
  > 
  >  out <- glm(y ~ x1 + x2, family = binomial())
  > 
  >  logl <- function(beta) {
  +      if (length(beta) != 3) stop("length(beta) != 3")
  +      beta0 <- beta[1]
  +      beta1 <- beta[2]
  +      beta2 <- beta[3]
  +      eta <- beta0 + beta1 * x1 + beta2 * x2
  +      p <- exp(eta) / (1 + exp(eta))
  +      return(sum(log(p[y == 1])) + sum(log(1 - p[y == 0])))
  +  }
  > 
  >  out.metro <- metrop(logl, coefficients(out), 1e3, scale = 0.01)
  >  out.metro$accept
  [1] 0.982
  > 
  >  out.metro <- metrop(out.metro, scale = 0.1)
  >  out.metro$accept
  [1] 0.795
  > 
  >  out.metro <- metrop(out.metro, scale = 0.5)
  >  out.metro$accept
  [1] 0.264
  > 
  >  apply(out.metro$batch, 2, mean)
  [1] 0.06080257 1.42304941 0.52634149
  > 
  >  out.metro <- metrop(logl, as.numeric(coefficients(out)), 1e2,
  +      scale = 0.5, debug = TRUE, blen = 5, nspac = 3)
  > 
  >  niter <- out.metro$nbatch * out.metro$blen * out.metro$nspac
  >  niter == nrow(out.metro$current)
  [1] TRUE
  >  niter == nrow(out.metro$proposal)
  [1] TRUE
  >  all(out.metro$current[1, ] == out.metro$initial)
  [1] TRUE
  >  all(out.metro$current[niter, ] == out.metro$final) |
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  [1] TRUE
  > 
  >  .Random.seed <- out.metro$initial.seed
  >  d <- ncol(out.metro$proposal)
  >  n <- nrow(out.metro$proposal)
  >  my.proposal <- matrix(NA, n, d)
  >  my.u <- double(n)
  >  ska <- out.metro$scale
  >  for (i in 1:n) {
  +      my.proposal[i, ] <- out.metro$current[i, ] + ska * rnorm(d)
  +      if (is.na(out.metro$u[i])) {
  +          my.u[i] <- NA
  +      } else {
  +          my.u[i] <- runif(1)
  +      }
  +  }
  >  max(abs(out.metro$proposal - my.proposal)) < epsilon
  [1] TRUE
  >  all(is.na(out.metro$u) == is.na(my.u))
  [1] TRUE
  >  all(out.metro$u[!is.na(out.metro$u)] == my.u[!is.na(my.u)])
  [1] TRUE
  > 
  >  my.curr.log.green <- apply(out.metro$current, 1, logl)
  >  my.prop.log.green <- apply(out.metro$proposal, 1, logl)
  >  all(is.na(out.metro$u) == (my.prop.log.green > my.curr.log.green))
  [1] TRUE
  >  foo <- my.prop.log.green - my.curr.log.green
  >  max(abs(foo - out.metro$log.green)) < epsilon
  [1] TRUE
  > 
  >  my.accept <- is.na(my.u) | my.u < exp(foo)
  >  sum(my.accept) == round(n * out.metro$accept)
  [1] TRUE
  >  if (my.accept[niter]) {
  +      all(out.metro$proposal[niter, ] == out.metro$final)
  +  } else {
  +      all(out.metro$current[niter, ] == out.metro$final)
  +  }
  [1] TRUE
  > 
  >  my.current <- out.metro$current
  >  my.current[my.accept, ] <- my.proposal[my.accept, ]
  >  my.current <- rbind(out.metro$initial, my.current[- niter, ])
  >  max(abs(out.metro$current - my.current)) < epsilon
  [1] TRUE
  > 
  >  my.path <- matrix(NA, n, d)
  >  my.path[my.accept, ] <- out.metro$proposal[my.accept, ]
  >  my.path[! my.accept, ] <- out.metro$current[! my.accept, ]
  >  nspac <- out.metro$nspac
  > 
  >  my.path <- my.path[seq(nspac, niter, by = nspac), ]
  > 
  >  foom <- array(as.vector(t(my.path)), c(d, out.metro$blen, out.metro$nbatch))
  >  boom <- t(apply(foom, c(1, 3), mean))
  > 
  >  all(dim(boom) == dim(out.metro$batch))
  [1] TRUE
  >  max(abs(boom - out.metro$batch)) < epsilon
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.62    0.14    0.75 
Running the tests in 'tests/morph.R' failed.
Complete output:
  > library(mcmc)
  > isotropic <- mcmc:::isotropic
  > isotropic.logjacobian <- mcmc:::isotropic.logjacobian
  > 
  > # make sure morph identity works properly
  > TestMorphIdentity <- function(m.id) {
  +   ident.func <- function(x) x
  +   if (!all.equal(m.id$transform(1:10), 1:10))
  +     return(FALSE)
  +   if (!all.equal(m.id$inverse(1:10), 1:10))
  +     return(FALSE)
  +   x <- seq(-1,1, length.out=15)
  +   if (!all.equal(sapply(x, m.id$lud(function(x) dnorm(x, log=TRUE))),
  +                  dnorm(x, log=TRUE)))
  +     return(FALSE)
  +   if (!all.equal(m.id$outfun(ident.func)(x), x))
  +     return(FALSE)
  +   return(TRUE)
  + }
  > 
  > TestMorphIdentity(morph())
  [1] TRUE
  > TestMorphIdentity(morph.identity())
  [1] TRUE
  > 
  > TestMorphIdentityOutfun <- function(m) {
  +   f <- m$outfun(NULL)
  +   x <- 1:20
  +   if (!identical(x, f(x)))
  +     return(FALSE)
  +   f <- m$outfun(c(6, 8))
  +   if (!identical(x[c(6, 8)], f(x)))
  +     return(FALSE)
  +   i <- rep(FALSE, 20)
  +   i[c(1, 3, 5)] <- TRUE
  +   f <- m$outfun(i)
  +   if (!identical(x[i], f(x)))
  +     return(FALSE)
  +   return(TRUE)
  + }
  > 
  > TestMorphIdentityOutfun(morph())
  [1] TRUE
  > TestMorphIdentityOutfun(morph.identity())
  [1] TRUE
  > 
  > # make sure that morph and morph.identity give back the same things
  > all.equal(sort(names(morph.identity())), sort(names(morph(b=1))))
  [1] TRUE
  > 
  > # test center parameter, univariate version
  > zero.func <- function(x) 0
  > center <- 2
  > x <- seq(-1,1, length.out=15)
  > morph.center <- morph(center=center)
  > all.equal(sapply(x, morph.center$transform), x-center)
  [1] TRUE
  > all.equal(sapply(x, morph.center$inverse), x+center)
  [1] TRUE
  > all.equal(sapply(x, morph.center$lud(function(y) dnorm(y, log=TRUE))),
  +           dnorm(x, log=TRUE, mean=-2))
  [1] TRUE
  > 
  > # test center parameter, multivariate version
  > center <- 1:4
  > x <- rep(0, 4)
  > morph.center <- morph(center=center)
  > lud.mult.dnorm <- function(x) prod(dnorm(x, log=TRUE))
  > all.equal(morph.center$transform(x), x-center)
  [1] TRUE
  > all.equal(morph.center$inverse(x), x+center)
  [1] TRUE
  > all.equal(morph.center$lud(lud.mult.dnorm)(x),
  +           lud.mult.dnorm(x - center))
  [1] TRUE
  > # test 'r'.
  > r <- 1
  > morph.r <- morph(r=r)
  > x <- seq(-1, 1, length.out=20)
  > all.equal(sapply(x, morph.r$lud(function(x) dnorm(x, log=TRUE))),
  +           dnorm(x, log=TRUE))
  [1] TRUE
  > x <- seq(1.1, 2, length.out=10)
  > all(sapply(x, morph.r$lud(function(x) dnorm(x, log=TRUE)))
  +     !=
  +     dnorm(x, log=TRUE))
  [1] TRUE
  > 
  > TestExponentialEvenPWithRInverse <- function() {
  +   r <- 0.3
  +   p <- 2.2
  +   morph.r <- morph(r=r, p=p)
  +   x <- seq(0, r, length.out=20)
  +   all.equal(x, sapply(x, morph.r$inverse))
  + }
  > 
  > TestExponentialEvenPWithRInverse()
  [1] TRUE
  > 
  > # make sure morph$lud passes '...' arguments.
  > mean <- 2
  > ident.morph <- morph()
  > dnorm.morph <- ident.morph$lud(function(x, mean=0)
  +                                  dnorm(x, mean=mean, log=TRUE))
  > all.equal(dnorm.morph(2, mean), dnorm(2, mean=mean, log=TRUE))
  [1] TRUE
  > x <- seq(-3, 3, length.out=20)
  > m2 <- morph(r=10)
  > dnorm.morph <- m2$lud(function(x, mean)
  +                         dnorm(x, mean=mean, log=TRUE))
  > all.equal(sapply(x, function(y) dnorm.morph(y, 2)),
  +           dnorm(x, mean=2, log=TRUE))
  [1] TRUE
  > 
  > # make sure morph$outfun passes '...' arguments.
  > outfun.orig <- function(x, mean) x + mean
  > ident.morph <- morph()
  > mean <- 1
  > outfun.morph <- ident.morph$outfun(outfun.orig)
  > all.equal(outfun.morph(1:10, mean), 1:10+mean)
  [1] TRUE
  > 
  > m2 <- morph(r=10)
  > outfun.morph <- m2$outfun(outfun.orig)
  > all.equal(sapply(1:10, function(x) outfun.morph(x, mean)), 1:10+mean)
  [1] TRUE
  > 
  > ###########################################################################
  > # test built-in exponential and polynomial transformations.
  > f <- morph(b=3)
  > x <- seq(0, 10, length.out=100)
  > all.equal(x, sapply(sapply(x, f$transform), f$inverse))
  [1] TRUE
  > 
  > f <- morph(p=3)
  > all.equal(x, sapply(sapply(x, f$transform), f$inverse))
  [1] TRUE
  > 
  > f <- morph(p=3, r=10)
  > all.equal(-10:10, Vectorize(f$transform)(-10:10))
  [1] TRUE
  > 
  > f <- morph(p=3, b=1)
  > all.equal(x, sapply(sapply(x, f$transform), f$inverse))
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.23    0.04    0.26 
Running the tests in 'tests/morph.metrop.R' failed.
Complete output:
  > library(mcmc)
  > 
  > .morph.unmorph <- mcmc:::.morph.unmorph
  > 
  > ###########################################################################
  > # basic functionality check, can morph.metro run?  Can we change the
  > # transformation?
  > set.seed(42)
  > obj <- morph.metrop(function(x) dt(x, df=3, log=TRUE),
  +                     100, 100, morph=morph(b=3))
  > obj <- morph.metrop(obj, morph=morph(b=1))
  > 
  > obj <- morph.metrop(function(x) prod(dt(x, df=3, log=TRUE)),
  +                     rep(100, 3), 100, morph=morph(p=3, b=1))
  > obj <- morph.metrop(obj, morph=morph(r=1, p=3, b=1))
  > 
  > all.equal(class(obj), c("mcmc", "morph.metropolis"))
  [1] TRUE
  > 
  > ###########################################################################
  > # check .morph.unmorph
  > obj <- list(final=10)
  > outfun <- function(x) x
  > m <- morph(p=3)
  > obj <- .morph.unmorph(obj, m, outfun)
  > all.equal(class(obj), c("mcmc", "morph.metropolis"))
  [1] TRUE
  > all.equal(sort(names(obj)),
  +           sort(c("final", "morph", "morph.final", "outfun")))
  [1] TRUE
  > all.equal(c(obj$final, obj$morph.final), c(m$inverse(10), 10))
  [1] TRUE
  > all.equal(obj$outfun, outfun)
  [1] TRUE
  > all.equal(obj$morph, m)
  [1] TRUE
  > 
  > proc.time()
     user  system elapsed 
     0.32    0.09    0.40 
Running the tests in 'tests/saveseedmorph.R' failed.
Complete output:
  > 
  >  library(mcmc)
  > 
  >  set.seed(42)
  > 
  >  h <- function(x) if (all(x >= 0) && sum(x) <= 1) return(1) else return(-Inf)
  >  out <- morph.metrop(obj = h, initial = rep(0, 5), nbatch = 100, blen = 17,
  +      nspac = 3, scale = 0.1)
  > 
  >  out1 <- morph.metrop(out)
  >  out2 <- morph.metrop(out1)
  >  out3 <- morph.metrop(out, nbatch = 2 * out$nbatch)
  > 
  >  fred <- rbind(out1$batch, out2$batch)
  >  identical(fred, out3$batch)
  [1] TRUE
  > 
  >  out <- morph.metrop(out, morph = morph(p = 2.2, r = 0.3))
  > 
  >  out1 <- morph.metrop(out)
  >  out2 <- morph.metrop(out1)
  >  out3 <- morph.metrop(out, nbatch = 2 * out$nbatch)
  > 
  >  fred <- rbind(out1$batch, out2$batch)
  >  identical(fred, out3$batch)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     0.64    0.04    0.68 
Running the tests in 'tests/temp-par.R' failed.
Complete output:
  > 
  >  library(mcmc)
  > 
  >  set.seed(42)
  > 
  >  data(foo)
  >  attach(foo)
  > 
  >  out <- glm(y ~ x1 + x2 + x3, family = binomial, x = TRUE)
  > ## IGNORE_RDIFF_BEGIN
  >  summary(out)
  
  Call:
  glm(formula = y ~ x1 + x2 + x3, family = binomial, x = TRUE)
  
  Coefficients:
              Estimate Std. Error z value Pr(>|z|)    
  (Intercept)   0.5772     0.2766   2.087 0.036930 *  
  x1            0.3362     0.4256   0.790 0.429672    
  x2            0.8475     0.4701   1.803 0.071394 .  
  x3            1.5143     0.4426   3.422 0.000622 ***
  ---
  Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
  
  (Dispersion parameter for binomial family taken to be 1)
  
      Null deviance: 134.602  on 99  degrees of freedom
  Residual deviance:  86.439  on 96  degrees of freedom
  AIC: 94.439
  
  Number of Fisher Scoring iterations: 5
  
  > ## IGNORE_RDIFF_END
  > 
  >  modmat <- out$x
  > 
  >  models <- cbind(rep(0:1, each = 4), rep(rep(0:1, times = 2), each = 2),
  +                rep(0:1, times = 4))
  > 
  >  exes <- paste("x", 1:3, sep = "")
  >  betas <- NULL
  >  for (i in 1:nrow(models)) {
  +      inies <- as.logical(models[i, ])
  +      foo <- exes[inies]
  +      bar <- paste("y ~", paste(foo, collapse = " + "))
  +      if (! any(inies)) bar <- paste(bar, "1")
  +      baz <- glm(as.formula(bar), family = binomial)
  +      beta <- rep(0, 4)
  +      beta[c(TRUE, inies)] <- baz$coef
  +      betas <- rbind(betas, beta)
  +  }
  > 
  >  neighbors <- matrix(FALSE, nrow(models), nrow(models))
  >  for (i in 1:nrow(neighbors)) {
  +      for (j in 1:ncol(neighbors)) {
  +          foo <- models[i, ]
  +          bar <- models[j, ]
  +          if (sum(foo != bar) == 1) neighbors[i, j] <- TRUE
  +      }
  +  }
  > 
  >  ludfun <- function(state, ...) {
  +      stopifnot(is.numeric(state))
  +      stopifnot(length(state) == ncol(models) + 2)
  +      stopifnot(length(state) == ncol(models) + 2)
  +      icomp <- state[1]
  +      stopifnot(icomp == as.integer(icomp))
  +      stopifnot(1 <= icomp && icomp <= nrow(models))
  +      beta <- state[-1]
  +      inies <- c(TRUE, as.logical(models[icomp, ]))
  +      beta.logl <- beta
  +      beta.logl[! inies] <- 0
  +      eta <- as.numeric(modmat %*% beta.logl)
  +      logp <- ifelse(eta < 0, eta - log1p(exp(eta)), - log1p(exp(- eta)))
  +      logq <- ifelse(eta < 0, - log1p(exp(eta)), - eta - log1p(exp(- eta)))
  +      logl <- sum(logp[y == 1]) + sum(logq[y == 0])
  +      val <- logl - sum(beta^2) / 2
  +      return(val)
  +  }
  > 
  >  ludval <- NULL
  >  for (i in 1:nrow(models)) ludval <- c(ludval, ludfun(c(i, betas[i, ])))
  >  all(is.finite(ludval))
  [1] TRUE
  > 
  > 
  >  out <- temper(ludfun, initial = betas, neighbors = neighbors, nbatch = 20,
  +      blen = 10, nspac = 5, scale = 0.56789, parallel = TRUE, debug = TRUE)
  > 
  >  names(out)
   [1] "lud"           "neighbors"     "nbatch"        "blen"         
   [5] "nspac"         "scale"         "outfun"        "debug"        
   [9] "parallel"      "initial.seed"  "final.seed"    "time"         
  [13] "batch"         "acceptx"       "accepti"       "initial"      
  [17] "final"         "which"         "unif.which"    "state"        
  [21] "log.hastings"  "unif.hastings" "proposal"      "acceptd"      
  [25] "norm"          "unif.choose"   "coproposal"   
  > 
  >  ### check decision about within-component or jump/swap
  > 
  >  identical(out$unif.which < 0.5, out$which)
  [1] TRUE
  > 
  >  identical(out$which, out$proposal[ , 1] == out$coproposal[ , 1])
  [1] TRUE
  > 
  >  ### check proposal and coproposal are actually current state or part thereof
  > 
  >  prop <- out$proposal
  >  coprop <- out$coproposal
  >  prop.i <- prop[ , 1]
  >  coprop.i <- coprop[ , 1]
  >  alt.prop <- prop
  >  alt.coprop <- coprop
  >  for (i in 1:nrow(prop)) {
  +      alt.prop[i, ] <- c(prop.i[i], out$state[i, prop.i[i], ])
  +      alt.coprop[i, ] <- c(coprop.i[i], out$state[i, coprop.i[i], ])
  +  }
  >  identical(coprop, alt.coprop)
  [1] TRUE
  >  identical(prop[! out$which, ], alt.prop[! out$which, ])
  [1] TRUE
  >  identical(prop[out$which, 1], alt.prop[out$which, 1])
  [1] TRUE
  > 
  >  ### check hastings ratio calculated correctly
  > 
  >  foo <- apply(prop, 1, ludfun)
  >  fooco <- apply(coprop, 1, ludfun)
  >  prop[ , 1] <- out$coproposal[ , 1]
  >  coprop[ , 1] <- out$proposal[ , 1]
  >  foo.swap <- apply(prop, 1, ludfun)
  >  fooco.swap <- apply(coprop, 1, ludfun)
  >  log.haste <- ifelse(out$which, foo - fooco,
  +      foo.swap + fooco.swap - foo - fooco)
  >  all.equal(log.haste, out$log.hastings)
  [1] TRUE
  > 
  >  ### check hastings rejection decided correctly
  > 
  >  identical(out$log.hastings >= 0, is.na(out$unif.hastings))
  [1] TRUE
  >  all(out$log.hastings < 0 | out$acceptd)
  [1] TRUE
  >  identical(out$acceptd,
  +      out$log.hastings >= 0 | out$unif.hastings < exp(out$log.hastings))
  [1] TRUE
  > 
  >  ### check acceptance carried out or not (according to decision) correctly
  > 
  >  before <- out$state
  >  after <- before
  >  after[- dim(after)[1], , ] <- before[-1, , ]
  >  after[dim(after)[1], , ] <- out$final
  >  my.after <- before
  >  for (i in 1:length(out$acceptd)) {
  +      if (out$acceptd[i]) {
  +          if (out$which[i]) {
  +              j <- out$proposal[i, 1]
  +              my.after[i, j, ] <- out$proposal[i, -1]
  +          } else {
  +              j <- out$proposal[i, 1]
  +              k <- out$coproposal[i, 1]
  +              my.after[i, j, ] <- out$coproposal[i, -1]
  +              my.after[i, k, ] <- out$proposal[i, -1]
  +          }
  +      }
  +  }
  >  identical(after, my.after)
  [1] TRUE
  > 
  >  ### check within-component proposal
  > 
  >  my.coproposal.within <- out$coproposal[out$which, ]
  >  proposal.within <- out$proposal[out$which, ]
  >  my.z <- out$norm[out$which, ]
  >  my.proposal.within <- my.coproposal.within
  >  my.proposal.within[ , -1] <- my.coproposal.within[ , -1] + out$scale * my.z
  >  all.equal(proposal.within, my.proposal.within)
  [1] TRUE
  > 
  >  my.unif.choose <- out$unif.choose[out$which, 1]
  >  my.i <- floor(nrow(models) * my.unif.choose) + 1
  >  all(1 <= my.i & my.i <= nrow(models))
  [1] TRUE
  >  identical(my.i, my.coproposal.within[ , 1])
  [1] TRUE
  > 
  >  ### check swap proposal
  > 
  >  coproposal.swap <- out$coproposal[! out$which, ]
  >  proposal.swap <- out$proposal[! out$which, ]
  >  unif.choose.swap <- out$unif.choose[! out$which, ]
  >  my.i <- floor(nrow(models) * unif.choose.swap[ , 1]) + 1
  >  nneighbors <- apply(out$neighbors, 1, sum)
  >  my.nneighbors <- nneighbors[my.i]
  >  my.k <- floor(my.nneighbors * unif.choose.swap[ , 2]) + 1
  >  my.j <- my.k
  >  foo <- seq(1, ncol(out$neighbors))
  >  for (i in seq(along = my.j)) {
  +      my.j[i] <- (foo[out$neighbors[my.i[i], ]])[my.k[i]]
  +  }
  >  identical(coproposal.swap[ , 1], my.i)
  [1] TRUE
  >  identical(proposal.swap[ , 1], my.j)
  [1] TRUE
  > 
  >  ### check standard normal and uniform random numbers are as purported
  > 
  >  save.Random.seed <- .Random.seed
  >  .Random.seed <- out$initial.seed
  > 
  >  nx <- ncol(out$initial)
  >  niter <- out$nbatch * out$blen * out$nspac
  >  my.norm <- matrix(NA, nrow = nrow(out$norm), ncol = ncol(out$norm))
  >  my.unif.which <- rep(NA, niter)
  >  my.unif.hastings <- rep(NA, niter)
  >  my.unif.choose <- matrix(NA, niter, 2)
  >  for (iiter in 1:niter) {
  +      my.unif.which[iiter] <- runif(1)
  +      if (out$which[iiter]) {
  +          my.unif.choose[iiter, 1] <- runif(1)
  +          my.norm[iiter, ] <- rnorm(nx)
  +          if (out$log.hastings[iiter] < 0) my.unif.hastings[iiter] <- runif(1) 
  +      } else {
  +          my.unif.choose[iiter, ] <- runif(2)
  +          if (out$log.hastings[iiter] < 0) my.unif.hastings[iiter] <- runif(1) 
  +      }
  +  }
  >  identical(my.norm, out$norm)
  [1] TRUE
  >  identical(my.unif.which, out$unif.which)
  [1] TRUE
  >  identical(my.unif.hastings, out$unif.hastings)
  [1] TRUE
  >  identical(my.unif.choose, out$unif.choose)
  [1] TRUE
  > 
  >  .Random.seed <- save.Random.seed
  > 
  >  ### check batch means
  > 
  >  foo <- after[seq(1, niter) %% out$nspac == 0, , ]
  >  foo <- array(as.vector(foo), dim = c(out$blen, out$nbatch, dim(foo)[2:3]))
  >  foo <- apply(foo, c(2, 3, 4), mean)
  >  all.equal(foo, out$batch)
  [1] TRUE
  > 
  >  ### check acceptance rates
  > 
  >  accept.within <- out$acceptd[out$which]
  >  my.i.within <- out$coproposal[out$which, 1]
  >  my.acceptx <- as.vector(sapply(split(accept.within, my.i.within), mean))
  >  identical(my.acceptx, out$acceptx)
  [1] TRUE
  > 
  >  accept.swap <- out$acceptd[! out$which]
  >  my.i.swap <- out$coproposal[! out$which, 1]
  >  my.j.swap <- out$proposal[! out$which, 1]
  >  nmodel <- nrow(out$neighbors)
  >  my.accepti <- matrix(NA, nmodel, nmodel)
  >  for (i in 1:nmodel) {
  +      for (j in 1:nmodel) {
  +          if (out$neighbors[i, j]) {
  +              my.accepti[i, j] <-
  +                  mean(accept.swap[my.i.swap == i & my.j.swap == j])
  +          }
  +      }
  +  }
  >  identical(my.accepti, out$accepti)
  [1] TRUE
  > 
  >  ### check scale vector
  > 
  >  nx <- ncol(models) + 1
  >  newscale <- rnorm(nx, 0.5, 0.1)
  > 
  >  out <- temper(out, scale = newscale)
  > 
  >  my.coproposal.within <- out$coproposal[out$which, ]
  >  proposal.within <- out$proposal[out$which, ]
  >  my.z <- out$norm[out$which, ]
  >  my.proposal.within <- my.coproposal.within
  >  my.proposal.within[ , -1] <- my.coproposal.within[ , -1] +
  +      sweep(my.z, 2, out$scale, "*")
  >  all.equal(proposal.within, my.proposal.within)
  [1] TRUE
  > 
  >  ### check scale matrix
  > 
  >  matscale <- matrix(rnorm(nx * nx, 0.0, 0.1), nx, nx)
  >  diag(matscale) <- 0.56789
  > 
  >  out <- temper(out, scale = matscale)
  > 
  >  my.coproposal.within <- out$coproposal[out$which, ]
  >  proposal.within <- out$proposal[out$which, ]
  >  my.z <- out$norm[out$which, ]
  >  my.proposal.within <- my.coproposal.within
  >  my.proposal.within[ , -1] <- my.coproposal.within[ , -1] +
  +      my.z %*% t(out$scale)
  >  all.equal(proposal.within, my.proposal.within)
  [1] TRUE
  > 
  >  ### check scale list
  > 
  >  lisztscale <- list(0.56789, newscale, matscale, matscale, newscale,
  +      0.98765, 0.98765, newscale)
  > 
  >  out <- temper(out, scale = lisztscale)
  > 
  >  my.coproposal.within <- out$coproposal[out$which, ]
  >  proposal.within <- out$proposal[out$which, ]
  >  my.z <- out$norm[out$which, ]
  >  my.proposal.within <- my.coproposal.within
  >  for (iiter in 1:nrow(my.z)) {
  +      my.i <- my.coproposal.within[iiter, 1]
  +      my.scale <- out$scale[[my.i]]
  +      if (is.matrix(my.scale)) {
  +          my.proposal.within[iiter, -1] <- my.coproposal.within[iiter, -1] +
  +              my.z[iiter, , drop = FALSE] %*% t(my.scale)
  +      } else {
  +          my.proposal.within[iiter, -1] <- my.coproposal.within[iiter, -1] +
  +              my.z[iiter, ] * my.scale
  +      }
  +  }
  >  all.equal(proposal.within, my.proposal.within)
  [1] TRUE
  > 
  >  ### check outfun
  > 
  >  outfun <- function(state, icomp, ...) {
  +      stopifnot(is.matrix(state))
  +      stopifnot(is.numeric(state))
  +      nx <- ncol(betas)
  +      ncomp <- nrow(betas)
  +      stopifnot(ncol(state) == nx)
  +      stopifnot(nrow(state) == ncomp)
  +      stopifnot(1 <= icomp && icomp <= ncomp)
  +      foo <- state[icomp, ]
  +      bar <- foo^2
  +      return(c(foo, bar))
  +  }
  > 
  >  out <- temper(out, outfun = outfun, icomp = 4)
  > 
  >  before <- out$state
  >  after <- before
  >  after[- dim(after)[1], , ] <- before[-1, , ]
  >  after[dim(after)[1], , ] <- out$final
  >  outies <- apply(after, 1, outfun, icomp = 4)
  >  outies <- t(outies)
  > 
  >  foo <- outies[seq(1, niter) %% out$nspac == 0, ]
  >  foo <- array(as.vector(foo), dim = c(out$blen, out$nbatch, dim(foo)[2]))
  >  foo <- apply(foo, c(2, 3), mean)
  >  all.equal(foo, out$batch)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.17    0.29    2.45 
Running the tests in 'tests/temp-ser-witch.R' failed.
Complete output:
  > 
  >  library(mcmc)
  > 
  >  set.seed(42)
  > 
  >  d <- 3
  >  witch.which <- 1 - (1 / 2)^(1 / d) * (1 / 4)^(seq(0, 5) / d)
  >  witch.which
  [1] 0.2062995 0.5000000 0.6850197 0.8015749 0.8750000 0.9212549
  > 
  >  ncomp <- length(witch.which)
  > 
  >  neighbors <- matrix(FALSE, ncomp, ncomp)
  >  neighbors[row(neighbors) == col(neighbors) + 1] <- TRUE
  >  neighbors[row(neighbors) == col(neighbors) - 1] <- TRUE
  >  neighbors[row(neighbors) == col(neighbors) + 2] <- TRUE
  >  neighbors[row(neighbors) == col(neighbors) - 2] <- TRUE
  > 
  >  ludfun <- function(state) {
  +      stopifnot(is.numeric(state))
  +      stopifnot(length(state) == d + 1)
  +      icomp <- state[1]
  +      stopifnot(icomp == as.integer(icomp))
  +      stopifnot(1 <= icomp && icomp <= ncomp)
  +      theta <- state[-1]
  +      if (any(abs(theta) > 1.0)) return(-Inf)
  +      bnd <- witch.which[icomp]
  +      if(bnd >= 1.0)
  +          stop(sprintf("witch.which[%d] >= 1.0", icomp))
  +      if(bnd <= 0.0)
  +          stop(sprintf("witch.which[%d] <= 0.0", icomp))
  +      if (all(abs(theta) > bnd))
  +          return(- (d + 1) * log(2) - d * log(1 - bnd))
  +      return(- (d + 1) * log(2) - log1p(- (1 - bnd)^d))
  +  }
  > 
  >  initial <- c(1, rep(0, d))
  > 
  >  out <- temper(ludfun, initial = initial, neighbors = neighbors,
  +      nbatch = 50, blen = 13, nspac = 7, scale = 0.3456789)
  > 
  >  names(out)
   [1] "lud"          "neighbors"    "nbatch"       "blen"         "nspac"       
   [6] "scale"        "outfun"       "debug"        "parallel"     "initial.seed"
  [11] "final.seed"   "time"         "batch"        "acceptx"      "accepti"     
  [16] "initial"      "final"        "ibatch"      
  > 
  >  out$acceptx
  [1] 0.6388889 0.4385246 0.3631714 0.4885246 0.4709677 0.4735516
  > 
  >  out$accepti
            [,1]      [,2]      [,3]      [,4]      [,5]      [,6]
  [1,]        NA 0.5071770 0.2727273        NA        NA        NA
  [2,] 0.7070064        NA 0.4355828 0.4186047        NA        NA
  [3,] 0.5816327 0.8039216        NA 0.5888889 0.5662651        NA
  [4,]        NA 0.7415730 0.8571429        NA 0.7857143 0.6626506
  [5,]        NA        NA 0.5204082 0.6516854        NA 0.8378378
  [6,]        NA        NA        NA 0.3515152 0.5056818        NA
  > 
  >  colMeans(out$ibatch)
  [1] 0.1830769 0.2153846 0.1630769 0.1369231 0.1353846 0.1661538
  > 
  >  ### check that have prob 1 / 2 for corners
  > 
  >  outfun <- function(state) {
  +      stopifnot(is.numeric(state))
  +      icomp <- state[1]
  +      stopifnot(icomp == as.integer(icomp))
  +      stopifnot(1 <= icomp && icomp <= length(witch.which))
  +      theta <- state[-1]
  +      foo <- all(abs(theta) > witch.which[icomp])
  +      bar <- rep(0, length(witch.which))
  +      baz <- rep(0, length(witch.which))
  +      bar[icomp] <- as.numeric(foo)
  +      baz[icomp] <- 1
  +      return(c(bar, baz))
  +  }
  > 
  >  out <- temper(out, blen = 103, outfun = outfun, debug = TRUE)
  > 
  >  eta.batch <- out$batch[ , seq(1, ncomp)]
  >  noo.batch <- out$batch[ , seq(ncomp + 1, ncomp + ncomp)]
  >  eta <- colMeans(eta.batch)
  >  noo <- colMeans(noo.batch)
  >  mu <- eta / noo
  >  eta
  [1] 0.06660194 0.06388350 0.05766990 0.06563107 0.10368932 0.22912621
  >  noo
  [1] 0.1365049 0.1258252 0.1293204 0.1370874 0.1716505 0.2996117
  >  mu
  [1] 0.4879090 0.5077160 0.4459459 0.4787535 0.6040724 0.7647440
  > 
  >  eta.batch.rel <- sweep(eta.batch, 2, eta, "/")
  >  noo.batch.rel <- sweep(noo.batch, 2, noo, "/")
  >  mu.batch.rel <- eta.batch.rel - noo.batch.rel
  > 
  >  mu.mcse.rel <- apply(mu.batch.rel, 2, sd) / sqrt(out$nbatch)
  >  mu.mcse.rel
  [1] 0.05010927 0.07897321 0.09678339 0.12636113 0.11261781 0.07082685
  > 
  >  foo <- cbind(mu, mu * mu.mcse.rel)
  >  colnames(foo) <- c("means", "MCSE")
  >  foo
           means       MCSE
  [1,] 0.4879090 0.02444876
  [2,] 0.5077160 0.04009596
  [3,] 0.4459459 0.04316016
  [4,] 0.4787535 0.06049584
  [5,] 0.6040724 0.06802931
  [6,] 0.7647440 0.05416441
  > 
  >  ### check decision about within-component or jump/swap
  > 
  >  identical(out$unif.which < 0.5, out$which)
  [1] TRUE
  > 
  >  identical(out$which, out$proposal[ , 1] == out$state[ , 1])
  [1] TRUE
  > 
  >  ### check hastings ratio calculated correctly
  > 
  >  n <- apply(neighbors, 1, sum)
  >  i <- out$state[ , 1]
  >  istar <- out$proposal[ , 1]
  >  foo <- apply(out$state, 1, ludfun)
  >  bar <- apply(out$proposal, 1, ludfun)
  >  my.log.hastings <- bar - foo - log(n[istar]) + log(n[i])
  >  all.equal(my.log.hastings, out$log.hastings)
  [1] TRUE
  > 
  > 
  > proc.time()
     user  system elapsed 
     2.73    0.10    2.82