windows - Windows 7下64位R的优化问题

Question

我目前正处于开发我的第一个 R 包的最后阶段，该包应该适合多项式处理树 (MPT) 模型(请参阅其 homepage 了解当前版本)。模型拟合是通过R的optim函数实现的。
今天我第一次在 Windows 7 机器上玩它并注意到一些非常奇怪的事情:optim 在使用 64 位版本的 R 时没有成功收敛。这对我来说像是一个错误(特别是当 nlminb 收敛于两个 R 版本时)。由于 optim 是我的包的核心，因此非常感谢对这个问题的任何帮助。

这是一个最小可重现的例子(通常模型是通过表达式指定的，而不是在目标函数中指定的，但为了简单起见，我把所有东西都放在目标函数中):

# 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))

此示例重现了 seminal paper on MPTs by Riefer & Batchelder 中的示例，即表 1 p 中的第 1 行。 327(预期参数值为 p = 1、q = .49 和 r = .30)。

在 32 位 R 上运行它总是给出正确的结果(尝试使用版本 2.12.2 和 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"

(请注意，由于随机起始值，计数可能会有所不同。)

另一方面，在 64 位 R 上运行它可能会产生这样的(错误的)结果:

$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"

目标函数的返回值和返回的参数值每次运行都不同，但计数总是3!

请注意，运行 nlminb 会在 32 位和 64 位 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)"

最后一点:我们有适用于 64 位 R 和 optim 的示例(这是我们最简单的示例模型)，但更多示例(如此处显示的示例)不起作用。

计数总是 3...

编辑:

固定起始值时(感谢 Joshua Ulrich)optim 不会偏离 64 位 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"

Answer 1

我们今天做了一些更多的测试，发现了与 Linux 下使用 64 位 R 的问题中描述的相同的问题。

但是，感谢Joachim Vandekerckhove对于这个巧妙的想法，我们尝试了一个简单的更改来解决问题(尽管问题仍然很可疑)。在目标函数的末尾，如果 llk 是 Inf，我们将它设置为一个非常高的值(是 1e19)。
使用较小的值(例如，1e10)可消除 64 位机器上的问题(目前已在 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))

这将返回正确的结果!