I have hit a problem with the design of the mcmc package I can't

figure out, possibly because I don't really understand the R function

call mechanism. The function metrop in the mcmc package has a ... argument

that it passes to one or two user-supplied functions, which are other

arguments to metrop. When the two functions don't have the same arguments,

this doesn't work. Here's an example.

library(mcmc)

library(MASS)

set.seed(42)

n <- 100

rho <- 0.5

beta0 <- 0.25

beta1 <- 0.5

beta2 <- 1

beta3 <- 1.5

Sigma <- matrix(rho, 3, 3)

diag(Sigma) <- 1

Sigma <- 0.75 * Sigma

Mu <- rep(0, 3)

foo <- mvrnorm(n, Mu, Sigma)

x1 <- foo[ , 1]

x2 <- foo[ , 2]

x3 <- foo[ , 3]

modmat <- cbind(1, foo)

eta <- beta0 + beta1 * x1 + beta2 * x2 + beta3 * x3

p <- 1 / (1 + exp(- eta))

y <- as.numeric(runif(n) < p)

out <- glm(y ~ x1 + x2 + x3, family = binomial())

summary(out)

### now we want to do a Bayesian analysis of the model, so we write

### a function that evaluates the log unnormalized density of the

### Markov chain we want to run (log likelihood + log prior)

ludfun <- function(beta) {

stopifnot(is.numeric(beta))

stopifnot(length(beta) == ncol(modmat))

eta <- as.numeric(modmat %*% beta)

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)

}

beta.initial <- as.vector(out$coefficients)

out <- metrop(ludfun, initial = beta.initial, nbatch = 20,

blen = 10, nspac = 5, scale = 0.56789)

### Works fine. Here are the Monte Carlo estimates of the posterior

### means for each parameter with Monte Carlo standard errors.

apply(out$batch, 2, mean)

sqrt(apply(out$batch, 2, function(x) initseq(x)$var.con) / out$nbatch)

### Now suppose I want Monte Carlo estimates of some function of

### the parameters other than the identity function. I write a function

### outfun that does that. Also suppose I want some extra arguments

### to outfun. This example is a bit forced, but I hit on natural

### examples with a new function (not yet released) that works like

### metrop but does simulated tempering.

outfun <- function(beta, degree) {

stopifnot(is.numeric(beta))

stopifnot(length(beta) == ncol(modmat))

stopifnot(is.numeric(degree))

stopifnot(length(degree) == 1)

stopifnot(degree == as.integer(degree))

stopifnot(length(degree) > 0)

result <- NULL

for (i in 1:degree)

result <- c(result, beta^i)

return(result)

}

out <- metrop(out, outfun = outfun, degree = 2)

### Oops! Try it and you get

###

### Error in obj(state, ...) : unused argument(s) (degree = 2)

I don't understand what the problem is (mostly because of ignorance). Because

foo <- function(x, ...) x

foo(x = 2, y = 3)

does work. The error is happening when ludfun is called, and I assume

the complaint is that it doesn't have an argument "degree", but then

why doesn't the simple example just above fail? So clearly I don't

understand what's going on.

An obvious solution is to ignore ... and just use global variables, i. e.,

define degree <- 2 in the global environment and make the signature of

outfun function(beta). That does work. But I don't want to have to

explain this issue on the help pages if I can actually fix the problem.

I have no idea whether one needs to look at the source code for the

mcmc package to diagnose the issue. If one does, it's on CRAN.

--

Charles Geyer

Professor, School of Statistics

University of Minnesota

[hidden email]
______________________________________________

[hidden email] mailing list

https://stat.ethz.ch/mailman/listinfo/r-devel