我目前正处于开发我的第一个R包的最后阶段，它应该适合多项处理树（MPT）模型（参见当前版本的主页 ）。 R的optim函数实现了模型拟合。
今天我第一次在Windows 7机器上玩了一下，发现了一些非常奇怪的东西：在使用64位版本的R时， optim不会成功地收敛。这看起来像是一个bug（尤其是当nlminb收敛于R版本）。 由于optim是我的软件包的核心，在这个问题上的任何帮助，不胜感激。
这里有一个最小可重现的例子（通常模型是通过expression式指定的，而不是在目标函数中指定的，但为了简单起见，我把所有的东西放在目标函数中）：
# 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))
这个例子再现了Riefer＆Batchelder在MPT上的开创性论文的一个例子，即表1中的第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"
目标函数的返回值和返回的参数值在每次运行中都不相同，但count总是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"
我们今天做了一些更多的测试，并发现了与使用64位R的Linux下的问题中所描述的相同的问题。
不过，感谢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))
这将返回正确的结果！