[Rd] question about ... passed to two different functions

[Charles Geyer]

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
charlie at stat.umn.edu