optim (method=Brent) 和 optimization 没有给出正确的二项式分布最小值 (N > 1000)

Question

我正在使用optim()(和optimize()) 尝试查找二项式分布的分位数，但是对于 N ~ 2000 (N = 2135)，函数不会给出正确的值。

optim(21, function(x) abs(1 - pbinom(x, 2135, 21/2135) - 0.1), 
      method = "Brent", lower = 1, upper = 2135) 

optimize(function(x) abs(1 - pbinom(x, 2135, 21/2135) - 0.1), c(1,2135))

PS：我也尝试将min参数设置为等于概率，但我仍然得到错误的答案。

Answer 1

问题在于optimize()，假设参数的微小变化将提供有关是否达到最小值（以及如果未达到最小值应朝哪个方向）的可靠信息。（我最初说该函数需要可微，这可能不是真的：请参阅有关布伦特方法的维基百科文章。）换句话说，大多数容易获得的优化算法可能会在分段常数的目标函数上失败，如下所示这个是...

IMO 对于这个几乎相同的问题所接受的答案是完全错误的。（它指出“起点处的梯度几乎为 0”，而实际上它恰好为零；正如您所发现的那样，使用optimize()并没有帮助，并且选择不同的起点或多或少取决于运气。 .)

我编了一个较小的例子来说明：找到 N=10、prob=0.2 的二项式分布的 0.6 分位数。R 可以直接、非常轻松地做到这一点：qbinom(0.6, size=10, prob=0.2)！但是假设您想解决类似形式的其他问题，这只是一个例子，或者约束是由家庭作业问题给出的，或者......

稍微简化的目标函数（使用平方差而不是绝对值）：

fx <- function(x) (pbinom(x, size=10, prob=0.2) - 0.6)^2

这看起来像什么？

curve(fx, from = 0, to =10, n=501)

所以正确的答案是 2 到 3 之间的任何值。在这种特殊情况下，optimize(fx, interval=c(1,10))恰好可以正常工作（返回 2.313，您可以将floor()其转换为 2），但如果我使用更宽的间隔（optimize(fx, interval=c(1,100))返回 99.99996），它将失败，或者如果我用更大的size. 让我尝试

fx2 <- function(x) pbinom(x, size=1000, prob=0.2) - 0.6
qbinom(0.6, size=1000, prob=0.2)  ## answer: 203
optimize(fx2, interval=c(1,1000)) ## 999.9999

问题是，如果优化方法的初始步骤跳跃少于一个单位，算法将得出结论：它已找到最小值。

一种可能的解决方案是寻找根而不是最小值：

fx3 <- function(x) pbinom(x, size=1000, prob=0.2) - 0.6
uniroot(fx3, interval=c(1,1000)) ## 203

我不知道解决这个优化问题的好方法。随机全局优化器可以工作，但通常效率非常低。请参阅此处，了解涉及 R 中非线性离散优化的一个特定问题。您还可以查看优化任务视图，尽管我发现它没有用...