R语言base::findInterval()函数异常行为的技术问询

一、异常行为的根本原因 #

这并非findInterval()的bug，而是浮点数精度误差导致的。十进制的0.06无法被二进制浮点数精确表示，seq(1/100, 1, 1/100)生成的第6个断点（理论值0.06）与直接输入的0.06在内存中的存储值存在微小差异：

nbins <- 100
probs <- seq(1/nbins, 1, 1/nbins)
# 验证浮点数差异
probs[6] == 0.06  # 返回 FALSE
sprintf("%.30f", probs[6])   # 输出类似 0.060000000000000005329
sprintf("%.30f", 0.06)       # 输出类似 0.059999999999999997780

当调用findInterval(0.06, probs)时，由于probs[6]略大于输入的0.06，函数会返回最大的索引i使得probs[i] <= 0.06，也就是5，而非预期的6。

二、解决浮点数精度问题的方案 #

方案1：改用整数运算（最优） #

由于你的p是整数累积和除以总数，可直接将比较转换为整数运算，完全规避浮点数误差：

ones <- rep(1, 2000)
p_int <- cumsum(ones)  # 整数累积和，避免浮点数
N <- sum(ones)
nbins <- 100
probs_int <- seq(N/nbins, N, N/nbins)  # 整数断点：20,40,...,2000

# 测试验证
x_int <- c(0.0095*N, 0.01*N, 0.0105*N, 0.011*N)
z_int <- c(0.0595*N, 0.06*N, 0.0605*N, 0.061*N)
findInterval(x_int, probs_int)  # 结果：[1] 0 1 1 1
findInterval(z_int, probs_int)  # 结果：[1] 5 6 6 6（符合预期）

方案2：给断点添加微小偏移量 #

如果必须使用浮点数，可给probs添加极小的偏移量，抵消精度误差：

probs_adj <- probs + 1e-12
findInterval(z, probs_adj)  # 结果：[1] 5 6 6 6（符合预期）

三、高效获取区间变化索引的方法 #

你当前的代码通过diff()查找区间变化，对于长向量可优化为直接利用findInterval()的反向查找，效率更高：

# 原方法
ind_old <- findInterval(p, probs) |> diff() |> { \(.) which(. == 1) }() + 1

# 优化方法：直接用probs反向查找p中的位置
ind_new <- findInterval(probs, p) + 1

优化原理 #

由于p和probs都是严格递增的，findInterval(probs[i], p)返回的是最大的索引k使得p[k] <= probs[i]，k+1就是p中第一个属于区间i+1的元素索引，与原方法结果完全一致，但时间复杂度从O(n)（n为p的长度）降至O(m log n)（m为probs的长度），在p非常长时优势明显。

内容的提问来源于stack exchange，提问作者R.Andres Castaneda