[Rd] question about ... passed to two different functions
Charles Geyer
charlie at stat.umn.edu
Sun Sep 6 19:50:46 CEST 2009
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
More information about the R-devel
mailing list