谐波平均值计算和浮点精度

我在PHP中实现了Pythagorean means，算术和几何平均值是小菜一碟，但我很难找到一个可靠的harmonic mean实现。谐波平均值计算和浮点精度

Harmonic Mean Definition from 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与相同ñ^-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 _我总和系列是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。

虽然键盘上的返回结果为-108086391056890000我的dev的机器（32位）说，这是-1.0808639105689E+17，它仍然是不一样的0或INF我所期待的。我甚至尝试在返回值上调用is_infinite()，但它按预期返回为false。

我还发现了stats_harmonic_mean()功能那是stats PECL扩展的一部分，但让我惊喜，我得到了完全相同的马车结果：-1.0808639105689E+17，如果任何参数为0，返回0但没有检查做的该系列的总和，as you can see on line 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扩展的情况下解决此问题？

来源

2012-05-04 Alix Axel

你是对的。您发现的数字是浮点算术特性的人为因素。

添加更多的精度不会帮助你。你所做的只是移动目标帖子。

底线是计算完成有限精度。这意味着在某个时候，中间结果将被舍入。那中间结果就不再确切。错误通过计算传播，并最终成为您的最终结果。当精确结果为零时，通常可以得到大约1e-16的双精度数字结果。

出现这种情况每次你的计算涉及与分母的分数是不是2

围绕它是表达的计算在整数或有理数的条款（如果可以的话），只有这样，一个电源时间，并使用任意精度的整数包进行计算。这就是Wolfram | Alpha所做的。

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

最后，一个元观察：Pythogarian意味着真正只对正数有意义。什么是3和-3的几何平均值？它是虚构的吗？链接到的维基百科页面上的不平等链只有在所有值都是正值时才有效。

来源

2012-05-04 05:30:37

非常好的答案和观察杰弗里，使用任意精度库做的伎俩，也四舍五入到最大精度（'圆（array_sum（$参数），ini_get（'精度'））'）返回'-0'这也可能是避免'gmp'或'bcmath'的依赖的好方法。关于你的元观察，你是对的。我应该过滤负值还是使用它们的绝对值？ –

@AlixAxel舍入是移动的目标职位。它可能适用于完全为零的值，但在某些时候，对于非常接近0的值会给出错误的结果。以'H（999999，-999998，-999997,999996）'为例。结果在'1e + 18'左右，但四舍五入到最大。双精度将给出0. –

@AlixAxel您如何处理负面投入取决于您的要求。如果纯粹是为了提供信息，那么我只会提出警告。 –

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

你可以使用你提到的使用通用分母方法（你发布的H2和H3公式），但是这只是在一段时间内将罐头踢出去，一旦总和就会得到不准确的结果产品术语开始四舍五入。

为什么要采取可能为负数的谐波平均值呢？这是一个固有的不稳定计算（H（-2,3,6 + epsilon）对于很小的epsilon变化很大）。

来源

2012-05-04 04:22:18

谢谢Keith，关于负数，我只是为了完整性而已，但我认为它没有多大意义。我应该过滤负数还是只使用它们的绝对值？ –

@AlixAxel：我会抛出一个异常，如果你可以在PHP中做到这一点。如果没有，则返回错误代码。默默地忽略不好的输入是一个坏主意。 –

谐波平均值计算和浮点精度

回答

相关问题