R devel - Sep 2009 - question about ... passed to two different functions

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

On 06/09/2009 1:50 PM, Charles Geyer wrote:
> 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.
The difference between foo and ludfun is that foo has a ... argument, 
and ludfun doesn't.  So when ludfun gets called with argument degree=2, 
it doesn't know what to do with it.  foo just puts it into the ... and 
ignores it.

Generally ... is very handy when there's only one function that may have 
optional arguments, but less so when you have two or more sets of 
optional arguments to handle.  You can use list(...) to get the list of 
args and spend a lot of work splitting them up, or you could tell users 
that if they want to pass optional arguments to outfun then ludfun 
should also be prepared to receive them, or you could add a parameter 
"ludcontrol" to be a list of things to pass to ludfun, and
"outcontrol"
to be a list of things to pass to outfun.  Or do away with optional 
parameters completely.

> 
> 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.
That's one solution.  A somewhat better one is to keep the values local, 
but still keep the signature of the function the same.  For example,

outfun <- local({  degree <- 2
                    function(beta) {
                       ...
                    }
                 })

or a function that creates outfun, e.g.

makeOutfun <- function(degree) {
   force(degree)  # make sure it gets evaluated
   function(beta) {  # Now create a function that can see it.
      ...
   }
}

These two solutions are better than a global in that you can have 
multiple outfuns with different degree values.
> 
> 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.
> 
I don't think so, but if my advice is way off base, maybe so.

Duncan Murdoch

Many thanks to those (Martin Morgan, Duncan Murdoch) who tried to straighten
me out on ... arguments.  It didn't work until I accidentally made two
examples
I thought were the same but one worked and the other didn't.  Finally I
achieved enlightenment.  The following section has been added to the
help page for the metrop function and will appear on CRAN when I get finished
with the simulated tempering function.

\section{Warning}{
If \code{outfun} is missing or not a function, then the log unnormalized
density can be defined without a \ldots argument and that works fine.
One can define it starting \code{ludfun <- function(state)} and that works
or \code{ludfun <- function(state, foo, bar)}, where \code{foo} and
\code{bar}
are supplied as additional arguments to \code{metrop}.

If \code{outfun} is a function, then both it and the log unnormalized
density function can be defined without \ldots arguments \emph{if they
have exactly the same arguments list} and that works fine.  Otherwise it
doesn't work.  Start the definitions \code{ludfun <- function(state,
foo)}
and \code{outfun <- function(state, bar)} and you get an error about
unused arguments.  Instead start the definitions
\code{ludfun <- function(state, foo, \ldots)}
and \code{outfun <- function(state, bar, \ldots)}, supply
\code{foo} and \code{bar} as additional arguments to \code{metrop},
and that works fine.

In short, the log unnormalized density function and \code{outfun} need
to have \ldots in their arguments list to be safe.  Sometimes it works
when \ldots is left out and sometimes it doesn't.

Of course, one can avoid this whole issue by always defining the log
unnormalized density function and \code{outfun} to have only one argument
\code{state} and use global variables (objects in the R global environment) to
specify any other information these functions need to use.  That too
follows the R way.  But some people consider that bad programming practice.
}

I hope that sums it up.  Apologies for submitting a rather stupid question
to the list.
-- 
Charles Geyer
Professor, School of Statistics
University of Minnesota
charlie at stat.umn.edu