使用optim()求二项检验临界区时N≥116不收敛问题

我来帮你拆解下这个问题：你想用optim()的Brent方法找二项检验的临界值，结果发现样本量N≥116时函数就卡壳了，没法收敛到正确结果，但N更小的时候一切正常。你的代码是这样的：

N <- 116 
optim(1, function(x) abs(1 - pbinom(x, N, 0.1) - 0.05), method = "Brent", lower = 1, upper = N)

为什么会出现这个问题？ #

核心原因是Brent方法是为连续函数设计的极小值查找工具，但你的目标函数本质是离散的——毕竟x是二项分布的计数，只能取整数。

• 当N比较小时，pbinom(x, N, 0.1)的变化步长（也就是每增加1个x，累积概率的变化量）比较大，目标函数的“起伏”很明显，Brent方法能轻松定位到最优的整数点；
• 但当N增大到116时，这个步长变得非常小，目标函数看起来就像“平缓的连续函数”，Brent方法没法识别出离散的跳跃点，自然就找不到精确的整数临界值了。

我们可以手动验证下N=116时的正确临界值：我们要找最小的x，使得上侧概率1 - pbinom(x, 116, 0.1)≤0.05。计算后会发现：

• 1 - pbinom(17, 116, 0.1) ≈ 0.048（符合≤0.05的要求）
• 1 - pbinom(16, 116, 0.1) ≈ 0.072（超过0.05）
  所以正确的临界值是17，但optim因为函数的离散性，没法精准定位到这个点。

怎么解决？ #

方案1：用R内置的qbinom()直接计算（最推荐）

R专门为离散分布提供了分位数函数qbinom()，完全可以替代手动用optim查找，效率高还准确：

N <- 116
# 找上侧0.05的临界值，等价于找最小的x使得P(X ≤ x) ≥ 0.95
critical_x <- qbinom(0.95, N, 0.1)

# 验证一下结果
1 - pbinom(critical_x - 1, N, 0.1) # 输出≈0.072，大于0.05
1 - pbinom(critical_x, N, 0.1)     # 输出≈0.048，小于等于0.05

这是最稳妥的方法，因为qbinom()本身就是为二项分布的分位数（临界值）设计的，完美适配你的需求。

方案2：修改optim的目标函数（仅当你一定要用optim时）

如果出于某些原因必须用optim()，可以修改目标函数，强制它只在整数点计算：

N <- 116
optim(1, function(x) {
  # 把连续的x转成整数，同时限制范围
  x_int <- round(x)
  if (x_int < 1 || x_int > N) return(Inf)
  # 只在整数点计算目标值
  abs(1 - pbinom(x_int, N, 0.1) - 0.05)
}, method = "Brent", lower = 1, upper = N)

不过这种方法不如qbinom()可靠，因为Brent方法还是针对连续函数设计的，可能会在整数点之间震荡，所以优先选方案1。

内容的提问来源于stack exchange，提问作者Ross