在 R 中使用 L-BFGS-B 方法最小化函数的问题？

Question

here 提出了类似的问题，但没有收到任何答案。

我有一个需要最小化的函数。我的函数有 3 个参数，x, y, and w。给定w，我需要找到最佳的x and y 对，以最小化我的功能。但是x and y被限制在(-0.5 ,0.5]区间内

我尝试使用optim 函数，当我没有限制时它可以工作。但是，当我应用限制并使用 L-BFGS-B 方法时，我得到一个错误。我不知道为什么会出现这个错误。

以下是我的功能示例和我尝试过的一些方法。这是我试图最小化的功能：

fn <- function(par, w) {
  x <- par[1]
  y <- par[2]

  p1 <- w^x * (w + 1) * beta(x + 1, w + 1 - x)
  p2 <- w^y * (w + 1) * beta(y + 1, w + 1 - y)
  p3 <- w^(2 * x) * (w + 1) * (beta(2 * x + 1, w + 1 - 2 * x) - (w + 1) * (beta(x + 1, w + 1 - x)^2))
  p4 <- w^(2 * y) * (w + 1) * (beta(2 * y + 1, w + 1 - 2 * y) - (w + 1) * (beta(y + 1, w + 1 - y)^2))
  p5 <- w^(y + x) * (w + 1) * (beta(y + x + 1, w + 1 - y - x) - (w + 1) * beta(y + 1, w + 1 - y) * beta(x + 1, w + 1 - x))
  p6 <- (y^2) * p3 * p1^(2 * y - 2) * p2^(-2 * x)
  p7 <- y * x * p1^(2 * y - 1) * p2^(-2 * x - 1) * p5
  p8 <- (x^2) * p4 * p1^(2 * y) * p2^(-2 * x - 2)
  mu <- (p1^y) / (p2^x)
  sigma <- p6 - 2 * p7 + p8
  v <- (gamma(x + 1)^y / gamma(y + 1)^x)

  quant <- ceiling((qnorm(0.95)*sqrt(sigma)/(v-mu))^2)
  
  return(quant)
}

这是我的问题：

# Here is an example of it working when there are no restrictions
oo = optim(par = c(0.1, 0.2), fn, w = 0.1)
oo$par
# In this example, I apply restrictions and it returns an error
oo1 = optim(par = c(0.1, 0.2), fn, w = 0.1, lower = c(-0.49, -0.49), upper = c(0.5, 0.5), method="L-BFGS-B")
oo1$par

在第二个示例中，我应用了限制并收到以下错误：

优化错误(par = c(0.1, 0.2), fn, w = 0.1, lower = c(-0.49, -0.49), : L-BFGS-B 需要 'fn' 的有限值

此外，我还尝试了optimx 包，它似乎不起作用。它只返回NAs。例如：

library(optimx)
opt <- optimx(par = c(0.1, 0.2), fn = fn, lower = c(-0.4, -0.4), upper = c(0.5, 0.5), w = 0.1)
opt
> opt
         p1 p2         value fevals gevals niter convcode kkt1 kkt2 xtime
L-BFGS-B NA NA 8.988466e+307     NA     NA    NA     9999   NA   NA 0.001

有没有人知道为什么我在应用限制时不能最小化我的功能？ （以及为什么optimx 包返回NAs？）

Answer 1

这是一个开始。我加了一行

cat(x,y, w, quant, "\n")

就在return(quant) 之前。您的第一次（无边界）运行：

0.1 0.2 0.1 40 
0.12 0.2 0.1 40 
0.1 0.22 0.1 41 
0.12 0.18 0.1 39 
0.13 0.16 0.1 39 
0.14 0.18 0.1 40 
0.11 0.195 0.1 40 
0.11 0.19 0.1 39 
0.12 0.19 0.1 40 
0.11 0.18 0.1 39 
0.1125 0.1825 0.1 39 

你的第二个（有界）给了：

0.1 0.2 0.1 40 
0.101 0.2 0.1 40 
0.099 0.2 0.1 39 
0.1 0.201 0.1 40 
0.1 0.199 0.1 39 
-0.49 -0.49 0.1 NaN 

换句话说，R 尝试在下界评估函数给出了NaN。

再四处寻找：

library(emdbook)
cc <- curve3d(fn(c(x,y), w = 0.1), xlim = c(-0.5,0.5), ylim = c(-0.5, 0.5),
              sys3d = "image")

如果 x=y 或 x=0 或 y=0 ，您的函数似乎给出了非有限值？您可能想捕捉这些情况并返回一个有限值...

事实上，这导致了一种解决方法，即不使下限或上限完全相等，以阻止 R 尝试在不好的地方评估函数：

oo1 = optim(par = c(0.1, 0.2), fn, w = 0.1, 
   lower = c(-0.49, -0.48), 
   upper = c(0.5, 0.49), method="L-BFGS-B")

也许可以通过查看目标函数并认真思考它在哪里会有非有限值来诊断，但是"thought is irksome and three minutes is a long time"