The problem with 64-bit R-optimism under Windows 7

I am currently in the final phase of developing my first R package, which should conform to the multidimensional processing (MPT) models (see its homepage for the current version). Model fitting is achieved by the R optim function.
Today I first spoke with him on a Windows 7 machine and noticed something really strange: optim does not succeed when using the 64-bit version of R. This looks like an error for me (especially how nlminb converges for both versions of R). Since optim is at the core of my package, any help on this is greatly appreciated.

Here is a minimally reproducible example (usually the model is specified using expressions and not specified in the objective function, but for simplicity I put everything in the objective function here):

 # The objective function: llk.tree <- function (Q, data) { p <- Q[1] q <- Q[2] r <- Q[3] e <- rep(NA,4) e[1] <- p * q * r e[2] <- p * q * (1-r) e[3] <- p * (1-q) * r e[4] <- p * (1-q) * (1-r) + (1-p) llk <- sum(data * log(e)) if (is.na(llk)) llk <- -1e+19 if (llk == -Inf) llk <- -1e+19 return(-llk) } # The call to optim: optim(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3)) 

This example reproduces the MPT seed example from Riefer and Batchelder , namely row 1 in table 1 s. 327 (the expected parameter values ​​will be p = 1, q = .49 and r = .30).

Running it on a 32-bit R always gives the correct result (verified with versions 2.12.2 and 2.13.0):

 $par [1] 1.0000000 0.4944449 0.3000001 $value [1] 234.7110 $counts function gradient 11 11 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" 

(Please note that the quantity may vary due to random starting values.)

Running it on 64 bit R, on the other hand, can lead to this (wrong) result:

 $par [1] 0.8668081 0.6326655 0.1433857 $value [1] 257.7328 $counts function gradient 3 3 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" 

The return value of the objective function and the return values ​​of the parameters are different for each run, but the score is always 3!

Note that running nlminb gives the correct results for 32-bit and 64-bit R:

 nlminb(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), lower = 0, upper = 1) $par [1] 1.0000000 0.4944445 0.3000000 $objective [1] 234.711 $convergence [1] 0 $iterations [1] 14 $evaluations function gradient 19 55 $message [1] "relative convergence (4)" 

Finally, we have examples (this is our simple sample model) that worked on 64-bit R and optim , but more examples (for example, shown here) did not work.

And the score is always 3 ...

EDIT:

When fixing the initial values ​​(thanks to Joshua Ulrich) optim does not optim from these fixed values ​​to 64 bits of R:

 optim(c(0.5, 0.5, 0.5), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3)) $par [1] 0.5 0.5 0.5 $value [1] 276.1238 $counts function gradient 3 3 $convergence [1] 0 $message [1] "CONVERGENCE: REL_REDUCTION_OF_F <= FACTR*EPSMCH" 
+4
source share
1 answer

We did some more tests today and found the same problem as in the Linux question using 64-bit R.

However, thanks to Joachim Vandekerckhove for this inventive idea, we tried a simple replacement that solved the problem (although the problem remains suspicious). At the end of the objective function, if llk Inf , we set it to an extremely high value (it was 1e19 ).
Using a smaller value (e.g. 1e10 ) fixes the problem on 64-bit machines (still tested on Linux):

 llk.tree <- function (Q, data) { p <- Q[1] q <- Q[2] r <- Q[3] e <- rep(NA,4) e[1] <- p * q * r e[2] <- p * q * (1-r) e[3] <- p * (1-q) * r e[4] <- p * (1-q) * (1-r) + (1-p) llk <- sum(data * log(e)) if (is.na(llk)) llk <- -1e+10 if (llk == -Inf) llk <- -1e+10 return(-llk) } # The call to optim: optim(runif(3, 0.1, 0.9), llk.tree, data = c(24, 65, 30, 61), method = "L-BFGS-B", lower = rep(0, 3), upper = rep(1, 3)) 

This returns the correct result!

+3
source

All Articles