MATLAB hist函数直方图算法及C++转换边缘偏移逻辑问询

问题背景 #

将MATLAB代码转换为C时，自定义的等效于[counts,centers]=hist(___)的直方图函数结果与原MATLAB代码不符，定位bug困难。使用MATLAB Coder生成C代码对比后，其中一段分箱边缘调整的代码逻辑不明，需解析其调整方式、设计意图，同时明确MATLAB hist函数的核心计算算法。

测试用MATLAB函数及调用脚本如下：

测试MATLAB函数 #

function [counts, centers] = my_hist(values, bins)
    [counts, centers] = hist(values, bins);
    disp(centers);
    disp(counts);
end

调用脚本 #

values = rand(1,1000);
bins = linspace(0.05, 0.95, 10);

[counts, centers] = my_hist(values, bins);

待解析C++代码逻辑分析 #

生成的C++核心代码片段用于调整分箱边缘，具体逻辑如下：

for (k = 0; k < 9; k++) {
    double absx_tmp;
    absx_tmp = edges[k + 1];
    absx = std::abs(absx_tmp);
    if ((!std::isinf(absx)) && (!std::isnan(absx))) {
      if (absx <= 2.2250738585072014E-308) {
        absx = 4.94065645841247E-324;
      } else {
        frexp(absx, &low_i);
        absx = std::ldexp(1.0, low_i - 53);
      }
    } else {
      absx = rtNaN;
    }
    edges[k + 1] = absx_tmp + absx;
}

边缘调整方式 #

这段代码的核心是将每个分箱的右边缘向右偏移一个该数值对应的最小可表示单位（ULP，Unit in the Last Place），具体步骤：

1. 遍历第2到第10个分箱边缘（对应输入10个bin中心生成的11个边缘中的中间9个右边缘）
2. 对当前边缘值，先判断是否为有效浮点数（非无穷、非NaN）：
   • 若边缘值绝对值小于等于2.2250738585072014E-308（double类型的最小正归一化数DBL_MIN），则直接使用4.94065645841247E-324（double类型的最小正非归一化数DBL_TRUE_MIN）作为偏移量
   • 否则，通过frexp拆分边缘值的尾数与指数，计算出该数值量级下的最小ULP（ldexp(1.0, low_i -53)，其中53是double类型的有效位数）
3. 将原边缘值加上计算出的偏移量，得到调整后的右边缘

设计意图 #

浮点数存在精度限制，当数据值刚好等于分箱边缘时，可能因二进制表示的微小误差导致归属错误（比如某值理论上等于边缘，但实际存储略大，被错误分到下一个箱子）。通过将右边缘向右偏移一个ULP，确保刚好等于原右边缘的值会被归入当前分箱，彻底避免浮点数精度误差引发的分箱错误，和MATLAB hist函数默认的「左闭右开（最后一个区间为闭区间）」规则完全匹配。

MATLAB hist函数核心算法 #

当输入为hist(values, bins)（bins为中心向量）时，算法流程如下：

1. 生成初始分箱边缘：
   • 若bins是等距向量：首边缘为第一个中心减去半间距，尾边缘为最后一个中心加上半间距；中间边缘为相邻两个中心的中点
   • 若bins是非等距向量：首边缘为第一个中心减去与第二个中心间距的一半，尾边缘为最后一个中心加上与倒数第二个中心间距的一半；中间边缘依次为相邻中心的中点
2. 分箱计数：
   • 对每个values中的元素，判断其落在哪个区间：前N-1个区间为左闭右开[edges(i), edges(i+1))，最后一个区间为闭区间[edges(N), edges(N+1)]（N为bin数量）
   • 统计每个区间内的元素数量，得到counts
3. 返回结果：counts为各分箱的元素数量，centers为输入的bins向量

若输入的bins是整数（表示分箱数量），则MATLAB会先根据values的范围自动计算等距的分箱中心，再执行上述流程。

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