谐均值计算和浮点精度

Question

我正在用PHP实现勾股平均数，算术平均数和几何平均数很容易实现，但是我很难想出一个可靠的调和平均数实现。

这是WolframAlpha的定义：

---

以下是在PHP中等效的实现：

function harmonicMeanV1()
{
    $result = 0;
    $arguments = func_get_args();

    foreach ($arguments as $argument)
    {
        $result += 1 / $argument;
    }

    return func_num_args() / $result;
}

现在，如果任何一个参数是0，这将引发除以0的警告。但由于1 / n与n^-1相同，并且pow(0, -1)优雅地返回INF常量而不会引发任何错误，因此我可以将其重写为以下形式（如果没有参数仍然会引发错误，但现在先忽略它）：

function harmonicMeanV2()
{
    $arguments = func_get_args();
    $arguments = array_map('pow', $arguments, array_fill(0, count($arguments), -1));

    return count($arguments) / array_sum($arguments);
}

对于大部分情况，两个实现都能很好地工作（例如v1、v2和WolframAlpha），但是当1 / n_i系列的总和为0时，它们会惨遭失败。我应该收到另一个除以0的警告，但事实并非如此......

考虑以下集合：-2, 3, 6（WolframAlpha说它是一个复杂的无限）：

  1 / -2    // -0.5
+ 1 / 3     // 0.33333333333333333333333333333333
+ 1 / 6     // 0.16666666666666666666666666666667

= 0

然而，我的两种实现都返回-2.7755575615629E-17作为总和（v1，v2），而不是0。

虽然CodePad上的返回结果是-108086391056890000，但我的开发机器（32位）显示为-1.0808639105689E+17，但它与我预期的0或INF完全不同。我甚至尝试在返回值上调用is_infinite()，但像预期的那样返回了false。

我还发现了stats_harmonic_mean()函数，它是stats PECL扩展的一部分，但令我惊讶的是，我得到了完全相同的错误结果：-1.0808639105689E+17，如果任何一个参数是0，则返回0，但没有对序列的总和进行检查，如您可以在第3585行看到的那样。

3557    /* {{{ proto float stats_harmonic_mean(array a)
3558       Returns the harmonic mean of an array of values */
3559    PHP_FUNCTION(stats_harmonic_mean)
3560    {
3561        zval *arr;
3562        double sum = 0.0;
3563        zval **entry;
3564        HashPosition pos;
3565        int elements_num;
3566    
3567        if (zend_parse_parameters(ZEND_NUM_ARGS() TSRMLS_CC, "a",  &arr) == FAILURE) {
3568            return;
3569        }
3570        if ((elements_num = zend_hash_num_elements(Z_ARRVAL_P(arr))) == 0) {
3571            php_error_docref(NULL TSRMLS_CC, E_WARNING, "The array has zero elements");
3572            RETURN_FALSE;
3573        }
3574    
3575        zend_hash_internal_pointer_reset_ex(Z_ARRVAL_P(arr), &pos);
3576        while (zend_hash_get_current_data_ex(Z_ARRVAL_P(arr), (void **)&entry, &pos) == SUCCESS) {
3577            convert_to_double_ex(entry);
3578            if (Z_DVAL_PP(entry) == 0) {
3579                RETURN_LONG(0);
3580            }
3581            sum += 1 / Z_DVAL_PP(entry);
3582            zend_hash_move_forward_ex(Z_ARRVAL_P(arr), &pos);   
3583        }
3584    
3585        RETURN_DOUBLE(elements_num / sum);
3586    }
3587    /* }}} */

这看起来像是一个典型的浮点精度问题，但我无法理解原因，因为每个计算都非常精确：

Array
(
    [0] => -0.5
    [1] => 0.33333333333333
    [2] => 0.16666666666667
)

是否有办法解决这个问题而不需要退回到使用gmp / bcmath扩展程序？

- Alix Axel

Answer 1

你是正确的。你发现的数字是浮点运算的特殊性质造成的产物。

增加更多精度并不会有所帮助。你只是在改变目标。

底线是计算是以有限的精度完成的。这意味着在某个时刻，中间结果将被四舍五入。那个中间结果就不再是精确的了。误差通过计算传播，并最终进入你的最终结果。当精确结果为零时，使用双精度数通常会得到约为1e-16的数值结果。

每当你的计算涉及分母不是2的幂的分数时，这种情况都会发生。

唯一的解决方法是用整数或有理数表达计算（如果可以），并使用任意精度整数包进行计算。这就是Wolfram|Alpha所做的。

请注意，计算几何平均值也不是简单的。尝试一个20次1e20的序列。由于所有数字都相同，结果应该是1e20。但你会发现结果是无穷大。原因是这20个数字的乘积（10e400）超出了双精度浮点数的范围，因此被设置为无穷大。无穷大的20次方根仍然是无穷大。

最后，一个元观察：毕达哥拉斯平均数只对正数有意义。3和-3的几何平均值是什么？它是虚数吗？你链接到的维基百科页面上的不等式链只在所有值都为正数时才有效。

- Jeffrey Sax

---

非常好的回答和观察，Jeffrey。使用任意精度库就可以解决问题，同时将结果四舍五入到最大精度(round(array_sum($arguments), ini_get('precision')))会返回-0，这也是避免依赖于gmp或bcmath的好方法。关于你的元观察，你是对的。我应该只过滤负值还是使用它们的绝对值？ - Alix Axel

@AlixAxel 四舍五入是移动目标杆。对于确切为零的值可能有效，但对于非常接近0的值，在某些时候会给出错误的结果。以H(999999,-999998,-999997,999996)为例。结果约为1e+18，但将其四舍五入到最大双精度会得到0。 - Jeffrey Sax

@AlixAxel 如何处理负面输入取决于你的需求。如果仅仅是为了提供信息，那么我只会给出警告。 - Jeffrey Sax

Answer 2

是的，这是浮点精度问题。-1/2 可以准确地表示，但 1/3 和 1/6 不能。因此当你把它们加起来时，你并不能得到零。

你可以采用你提到的通分的方法（你发布的 H2 和 H3 公式），但这只会把问题推迟一些时间，一旦积和之和开始四舍五入，你仍然会得到不准确的结果。

为什么你要计算可能是负数的数字的调和平均值呢？这是一个不稳定的计算（H(-2,3,6+epsilon)对于非常小的 epsilon 会变化很大）。

- Keith Randall

---

谢谢Keith，关于负数，我只是想要完整性，但我认为这没有太多意义。我应该过滤负数还是只使用它们的绝对值？ - Alix Axel

@AlixAxel：如果在PHP中可以的话，我会抛出一个异常。如果不行，就返回一个错误代码。默默地忽略错误输入是一个坏主意。 - Keith Randall