针对某些初始值产生 NaN 的函数的优化答案

【问题标题】：Optimization for functions that produce NaN for some initial values针对某些初始值产生 NaN 的函数的优化
【发布时间】：2013-01-11 07:32:58
【问题描述】：

我想找到以下目标函数的所有局部最小值

func <- function(b){Mat=matrix(c(+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2),2,2);d=(det(Mat));return(d)}

'func' 是 Logistic 回归模型的 Fisher 信息矩阵的行列式，是参数 b1 和 b2 的函数，其中 b1 属于 [-.3, .3]，b2 属于 [6, 8]

假设这两个初始值为 b = c(b1, b2)

> in1 <- c(-0.04785405, 6.42711047)
> in2 <- c(0.2246729, 7.5211575)

初始值为in1的局部最小值为：

> optim(in1, fn = func, lower = c(-.3, 6), upper = c(.3, 8), method = "L-BFGS-B")

$par
[1] -0.04785405  6.42711047

$value
[1] 3.07185e-27

$counts
function gradient 
   1        1 

$convergence
[1] 52

$message
[1] "ERROR: ABNORMAL_TERMINATION_IN_LNSRCH"

从$massage 中可以看出，优化过程中发生了终止，无法计算最小值，optim 返回 in1 作为局部最优值。

对于 'in2' 也会出现错误：

> optim(in2, fn = func, lower = c(-.3, 6), upper = c(.3, 8), method = "L-BFGS-B")

Error in optim(in2, fn = func, lower = c(-0.3, 6), upper = c(0.3, 8),  : 
L-BFGS-B needs finite values of 'fn'

发生此错误是因为func 的值为in2' isNaN`：

> func(in2)
[1] NaN

但是对于in1，目标函数在in1 的值被计算，但优化终止，因为optim 无法继续计算另一个初始值：

> func(in1)
[1] 3.07185e-27

让我定义没有 det 的 func ，就像矩阵一样，看看发生了什么：

Mat.func <- function(b){Mat=matrix(c(+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5)/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5)/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2,+0.5*1/((1/(exp(-b[1]-b[2]*-5)+1))*(1-(1/(exp(-b[1]-b[2]*-5)+1))))*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2*exp(-b[1] - b[2] * -5) * -5/(exp(-b[1] - b[2] * -5) + 1)^2+0.5*1/((1/(exp(-b[1]-b[2]*5)+1))*(1-(1/(exp(-b[1]-b[2]*5)+1))))*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2*exp(-b[1] - b[2] * 5) * 5/(exp(-b[1] - b[2] * 5) + 1)^2),2,2);d=Mat;return(d)}

我们得到

         > Mat.func(in1)
              [,1]         [,2]
         [1,] 1.109883e-14 2.784007e-15
         [2,] 2.784007e-15 2.774708e-13

        > Mat.func(in2)
              [,1] [,2]
          [1,]  Inf  Inf
          [2,]  Inf  Inf

因此，按双精度计算，Mat.func(in2) 元素的值是 Inf。我还用mpfr函数重写了Mat.func：

Mat.func.mpfr <-function(b, prec){ d=c(+0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2+0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2,
                                   +0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) * -5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2+0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) * 5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2,
                                   +0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) * -5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2+0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5)/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) * 5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2,
                                   +0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*-5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) * -5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) * -5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * -5) + 1)^2+0.5*1/((1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))*(1-(1/(exp(-mpfr(b[1], precBits = prec)-mpfr(b[2], precBits = prec)*5)+1))))*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) * 5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2*exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) * 5/(exp(-mpfr(b[1], precBits = prec) - mpfr(b[2], precBits = prec) * 5) + 1)^2)
                               Mat = new("mpfrMatrix", d, Dim = c(2L, 2L))
                               return(Mat)}

因此：

require(Rmpfr)
> Mat.func.mpfr(c(in1), prec = 54)
'mpfrMatrix' of dim(.) =  (2, 2) of precision  54   bits 
     [,1]                   
 [1,] 1.10988301365972506e-14
 [2,] 2.78400749725484580e-15
      [,2]                   
 [1,] 2.78400749725484580e-15
 [2,] 2.77470753414931256e-13

 > Mat.func.mpfr(c(in2), prec = 54)
 'mpfrMatrix' of dim(.) =  (2, 2) of precision  54   bits 
      [,1] [,2]
 [1,]  Inf  Inf
 [2,]  Inf  Inf

 > Mat.func.mpfr(c(in2), prec = 55)
 'mpfrMatrix' of dim(.) =  (2, 2) of precision  55   bits 
      [,1]                    
 [1,]  4.16032108702067276e-17
 [2,] -8.34300174643550123e-17
      [,2]                    
 [1,] -8.34300174643550154e-17
 [2,]  1.04008027175516816e-15

因此，精确到 55，矩阵元素的值不再是 Inf。很遗憾， mpfr 函数改变了目标的类，det 也没有 r 优化函数不能应用，为了澄清，我提供了两个例子：

> class(mpfr (1/3, 54))
[1] "mpfr"
attr(,"package")
[1] "Rmpfr"

## determinant
example1 <- function(x){
  d <- c(mpfr(x, prec = 54), 3 * mpfr(x, prec = 54), 5 * mpfr(x, prec = 54), 7 * mpfr(x, prec = 54))
  Mat = new("mpfrMatrix", d, Dim = c(2L, 2L))
  return(det(Mat))
}

> example1(2)
Error in UseMethod("determinant") : 
no applicable method for 'determinant' applied to an object of class "c('mpfrMatrix',    'mpfrArray', 'Mnumber', 'mNumber', 'mpfr', 'list', 'vector')"

##optimization 
example2 <- function(x)  ## Rosenbrock Banana function
   100 * (mpfr(x[2], prec = 54) - mpfr(x[1], prec = 54) * mpfr(x[1], prec = 54 ))^2 + (1 - mpfr(x[1], prec = 54))^2

> example2(c(-1.2, 1))
1 'mpfr' number of precision  54   bits 
[1] 24.1999999999999957
> optim(c(-1.2,1), example2)
Error in optim(c(-1.2, 1), example2) : 
(list) object cannot be coerced to type 'double'

因此，使用 mpfr 无法解决问题。

要找到所有局部最小值，应该编写一个应用不同随机初始值的算法。但是可以看出，对于某些初始值，该函数会产生NaN （忽略这些值并不是一个好主意，因为它通常会导致丢失一些局部最小值，特别是对于具有大量局部最优值的函数）。

我想知道是否有任何 R 包可以任意精度进行优化过程以避免目标函数的NaN？

谢谢

【问题讨论】：

您或许可以重新制定您的目标函数以获得更少的NaN 值（例如，通过重新排列以尽量减少下溢/上溢的可能性）
所以你试图最小化log(det(some_matrix))。 1) 每当 det(some_matrix) < 0 时，您都会得到 NaN，因为 log(x) 没有为 x < 0 定义；你对此有何意义？ 2）优化器会尝试找到det(some_matrix) == 0的位置；也许将您的目标更改为abs(det(some_matrix)) 将修复渐近行为。很难说你想做什么。
@flodel 实际上我在没有日志的情况下尝试过，但答案是一样的。抱歉，我在问题中写的 func(in1) 和 func(in2) 的值是针对 det(some_matrix) 而不是 log(det(some_matrix)) 我删除了 log

标签： r optimization nan arbitrary-precision

【解决方案1】：

回答您的问题，使用Rmpfr - 生成的矩阵：（虽然效率不高......！......）：

是的，determinant() 不适用于 mpfr 矩阵，然而你可以简单地使用类似的东西

M <- Mat.func.mpfr(in2, prec = 55)
m <- as(M, "matrix")
ldm <- determinant(m) # is already  log() !

然后使用

 { r <- determinant(., log=TRUE) ; if(r$sign <= 0) -Inf else r$modulus }

我在上面提到过......比使用 log(det(.)) 的“设计错误”要好得多

【讨论】：

【解决方案2】：

使用mpfr 有助于避免在函数中计算NaN（以及在优化算法中停止）。但 mpfr 输出是一个“mpfr”类，某些 R 函数（例如 optim 和 det）可能不适用于此类。像往常一样，as.numeric 可用于将“mpfr”类转换为“数字”类。

exp(9000)
[1] Inf

require(Rmpfr)
number <- as.numeric(exp(mpfr(9000, prec = 54)))

class(number)
[1] "numeric"

round(number)
[1] 1.797693e+308

number * 1.797692e-308
[1] 3.231699

number * 1.797693e-307
[1] 32.317

number * (1/number)
[1] 1

number * .2
[1] 3.595386e+307

number * .9
[1] 1.617924e+308

number * 1.1
[1] Inf

number * 2
[1] Inf

number / 2
[1] 8.988466e+307

number + 2
[1] 1.797693e+308

number + 2 * 10 ^ 291
[1] 1.797693e+308

number + 2 * 10 ^ 292
[1] Inf

number - 2
[1] 1.797693e+308

number - 2 * 10 ^ 307
[1] 1.597693e+308

number - 2 * 10 ^ 308
[1] -Inf

现在考虑以下矩阵函数：

mat <- function(x){
x1 <- x[1]
x2 <- x[2]
d = matrix(c(exp(5 * x1+ 4 * x2), exp(9 * x1), exp(2 * x2 + 4 * x1),
           exp(3 * x1)), 2, 2)
         return(d)
}

这个矩阵的元素很有可能产生Inf：

mat(c(300, 1))
    [,1] [,2]
[1,]  Inf  Inf
[2,]  Inf  Inf

所以如果在函数环境中返回det，我们得到的不是数字结果NaN，optim 函数肯定会被终止。为了解决这个问题，这个函数的行列式由mpfr编写：

func <- function (x){
  x1 <- mpfr(x[1], prec = precision)
  x2 <- mpfr(x[2], prec = precision)
  mat <- new("mpfrMatrix",c(exp(5 * x1+ 4 * x2), exp(9 * x1), exp(2 * x2 + 4 * x1),   exp(3 * x1)), Dim = c(2L,2L))
  d <- mat[1, 1] * mat[2, 2] - mat[2, 1] * mat[1, 2]
  return(as.numeric(-d))
}

那么对于 x1 = 3 和 x2 = 1，我们有：

func(c(3,1))
[1] 6.39842e+17

optim(c(3, 1),func)

$par
[1] 0.4500 1.4125

$value
[1] -4549.866

$counts
function gradient 
  13       NA 

$convergence
[1] 0

$message
NULL

对于 x1 = 300 和 x2 = 1：

func(c(300,1))
[1] 1.797693e+308

optim(c(300, 1),func)
$par
[1] 300   1

$value
[1] 1.797693e+308

$counts
function gradient 
   3       NA 

$convergence
[1] 0

$message
NULL

可以看出，没有停止，甚至optim 声称在优化过程中收敛。但是，似乎没有迭代，optim 只是将初始值作为局部最小值返回（当然，1.797693e+308 不是这个函数的局部最小值！！）。在这种情况下，应用mpfr 可以防止优化过程终止，但是如果我们真的期望优化算法从它们的值是Inf的这些点开始以R双精度并继续迭代以达到局部最小值，除了定义一个具有'mpfr'类的函数外，优化函数还应该有这种能力与“mpfr”类一起工作。

【讨论】：

随意编写一个附加包来执行此操作。一个困难是大多数 R 的优化器在内部使用较低级别的（C 或 FORTRAN）代码，这不容易适应使用mfpr。一种想法是您可以尝试optimx 包中的一些优化器，它们具有“纯R”版本的共轭梯度和准牛顿方法，可能适用于mfpr 对象。
@Ben Bolker 实际上我想编写这个包，我还选择了一个名为 Rsolnp 的合适包，它是纯 R 并且基于 Matlab 代码。

【解决方案3】：

我认为答案（我认为“agstudy”也给出了）是：确保您最小化的函数确实 NOT 返回 NaN（或 NA），而是返回 +Inf（如果您最小化，或 -Inf 如果你最大化）。

第二个：你应该真的使用 log(det(.)) 而不是
{ r

这也更准确。 {提示：看看 det 在 R 中是如何定义的！}

现在到Rmpfr，我会单独回复。应该像标准 R 一样使用“mpfr”-numbers， .... Rmpfr 的作者说 .... 但你可能需要一点小心。但是，不应需要 tryCatch()。

【讨论】：

【解决方案4】：

我试图重新制定你可怕的（对不起这个词）目标函数。我很确定 w 我们可以找到更简单的形式。希望其他人可以使用它来找到您的优化问题的解决方案...

func1 <- function(b){
  A <- exp(-b[1]+5*b[2])
  C <- exp(-b[1]-5*b[2])
  A1 <- A + 1
  C1 <- C + 1
  D <- 1/A1
  H <- 1/C1
  K <- D*(1-D)
  J <- H*(1-H)
  M <- (A/A1^2)^2/K
  N <- (C/C1^2)^2/J


Mat <- matrix(c( 1 *M    + 1  *N,
                -5 *M    + 5  *N,
                -5 *M    + 5  *N,
                25 *M    + 25 *N),2,2)

  Mat <- 0.5*Mat
  d <- log(det(Mat))
  return(d)
}

编辑

正如我所说，您可以再次简化您的功能。看起来好多了

func1 <- function(b){
  A <- exp(-b[1]+5*b[2])
  C <- exp(-b[1]-5*b[2])
  A1 <- A + 1
  C1 <- C + 1
  M <- A/A1^2
  N <- C/C1^2
  det.Mat <-25*M*N
  log(det.Mat)
}

这里有两个函数之间的一些测试。

func1(c(1,2))
[1] -16.7814
> func1(c(8,2))
[1] -17.03498
> func1(c(10,2))
[1] -18.16742
> func(c(10,2))
[1] -18.16742
> func(c(10,5))
[1] -46.83608

重新制定最小化了下溢/上溢的可能性（不能将中间结果存储在寄存器中）..这就是我们得到 Inf 而不是 NA（见下文）的原因，它是无限的，但仍然是 numeric，适用于与NaN相反的更远的计算，就像NA值..

函数(c(10,100))
[1] 南 func1(c(10,100)) [1] -Inf

现在我以更简单的形式测试您的优化指令，并且如您所见，它会收敛：

in1 <- c(-0.04785405, 6.42711047)
in2 <- c(0.2246729, 7.5211575)
ll <- optim(in1, fn = func1, lower = c(-.3, 6), upper = c(.3, 8), method = "L-BFGS-B")
 do.call(rbind,ll)


            function                                           gradient                                          
par         "-0.04785405"                                      "8"                                               
value       "-76.7811241751318"                                "-76.7811241751318"                               
counts      "2"                                                "2"                                               
convergence "0"                                                "0"                                               
message     "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL" "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

in2 也一样

optim(in2, fn = func1, lower = c(-.3, 6), upper = c(.3, 8), method = "L-BFGS-B")
$par
[1] 0.2246729 8.0000000

$value
[1] -76.78112

$counts
function gradient 
       2        2 

$convergence
[1] 0

$message
[1] "CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL"

【讨论】：

我同意我的函数公式很糟糕，感谢您的修改。其实不是我写的！！事实上，我应用了函数gsub、D和paste 来构建所有非线性模型的信息矩阵，而无需用户干预，并且您看到的公式是自动构建的。
@EhsanMasoudi 好的，但是当你在这里问答案时，你必须付出更多的努力......或者至少告诉我们你必须如何拥有这样的功能......据说我更新了我的答案更简单的功能。
如果你能解释为什么 func(c(10, 100)) 是 "Nan" 但 func(c(10, 100)) 是 `-Inf'，我将不胜感激
我添加一些解释...这是我对事物的解释..我不是数字专家..

【解决方案5】：

对于任意精度：gmp 和/或Rmpfr。不过，在代码中添加一些 tryCatch 可能会更好（以避免在给定尝试导致 NaN 错误时崩溃）

【讨论】：

我用 Rmpfr 编辑了这个问题，可以看出它没有用。此外，正如我所说，通过 tryCathch 忽略 NaN 可能会导致丢失一些局部最小值。
arb 精度 = 任意精度或 arb 库中的精度？
@Adam "arbitrary" ... 在我发布此内容时，我什至不知道 arb 包。