[Rd] question about --with-valgrind-instrumentation=level
Does --with-valgrind-instrumentation=2 slow down R when valgrind or gctorture are not in use? I am thinking of compiling the R that the whole department uses for research and teachin with --with-valgrind-instrumentation=2. Is that a good idea or a bad idea? -- Charles Geyer Professor, School of Statistics University of Minnesota char...@stat.umn.edu __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] question about --with-valgrind-instrumentation=level
On Sun, 6 Sep 2009, Charles Geyer wrote: Does --with-valgrind-instrumentation=2 slow down R when valgrind or gctorture are not in use? I am thinking of compiling the R that the whole department uses for research and teachin with --with-valgrind-instrumentation=2. Is that a good idea or a bad idea? It is a slightly bad idea: it does add extra code in memory.c that would always be called. This only affects R.bin (or libR.so if you use the already slower shared-library version), so you can so as I do and simply use a different version of R.bin. There is a similar issue with memory profiling. -- Brian D. Ripley, rip...@stats.ox.ac.uk Professor of Applied Statistics, http://www.stats.ox.ac.uk/~ripley/ University of Oxford, Tel: +44 1865 272861 (self) 1 South Parks Road, +44 1865 272866 (PA) Oxford OX1 3TG, UKFax: +44 1865 272595 __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
[Rd] 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
char...@stat.umn.edu
__
R-devel@r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] question about ... passed to two different functions
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 d
[Rd] R v2.9.2 patched still referred to as RC
R v2.9.2 patched is still referred to as RC: http://cran.r-project.org/bin/windows/base/rpatched.html Since v2.9.2 is already out, this should be 'pat', correct? /Henrik __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
Re: [Rd] R v2.9.2 patched still referred to as RC
On 06/09/2009 5:57 PM, Henrik Bengtsson wrote: R v2.9.2 patched is still referred to as RC: http://cran.r-project.org/bin/windows/base/rpatched.html Since v2.9.2 is already out, this should be 'pat', correct? The label was correct: that was still the RC version. The problem was that I forgot to switch the build script back to the regular patch builds. I've done it now, so in a few hours you should get something more useful... Duncan Murdoch __ R-devel@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel
