2012-05-04 66 views
6

我在PHP中實現了Pythagorean means,算術和幾何平均值是小菜一碟,但我很難找到一個可靠的harmonic mean實現。諧波平均值計算和浮點精度

這是WolframAlpha definition

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相同ñ-1pow(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); 
} 

兩種實現很好地工作在大多數情況下(例如v1v2WolframAlpha),但他們失敗了壯觀如果滿足的1/n 總和系列是0,我應該警告0拿到另一個部門,但我不...

考慮以下設置:-2, 3, 6WolframAlpha說,這是一個複雜的無限):

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,它仍然是不一樣的0INF我所期待的。我甚至嘗試在返回值上調用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擴展的情況下解決此問題?

回答

4

你是對的。您發現的數字是浮點算術特性的人爲因素。

添加更多的精度不會幫助你。你所做的只是移動目標帖子。

底線是計算完成有限精度。這意味着在某個時候,中間結果將被舍入。那中間結果就不再確切。錯誤通過計算傳播,並最終成爲您的最終結果。當精確結果爲零時,通常可以得到大約1e-16的雙精度數字結果。

出現這種情況每次你的計算涉及與分母的分數是不是2

圍繞它是表達的計算在整數或有理數的條款(如果可以的話),只有這樣,一個電源時間,並使用任意精度的整數包進行計算。這就是Wolfram | Alpha所做的。

請注意,幾何平均值的計算也不是微不足道的。嘗試20次1e20的序列。由於數字都是一樣的,結果應該是1e20。但是你會發現結果是無限的。原因在於這20個數字(10e400)的乘積超出了雙精度浮點數的範圍,因此它被設置爲無窮大。無限的第20根仍然無限。

最後,一個元觀察:Pythogarian意味着真正只對正數有意義。什麼是3和-3的幾何平均值?它是虛構的嗎?鏈接到的維基百科頁面上的不平等鏈只有在所有值都是正值時纔有效。

+0

非常好的答案和觀察傑弗裏,使用任意精度庫做的伎倆,也四捨五入到最大精度('圓(array_sum($參數),ini_get('精度'))')返回'-0'這也可能是避免'gmp'或'bcmath'的依賴的好方法。關於你的元觀察,你是對的。我應該過濾負值還是使用它們的絕對值? –

+0

@AlixAxel舍入是移動的目標職位。它可能適用於完全爲零的值,但在某些時候,對於非常接近0的值會給出錯誤的結果。以'H(999999,-999998,-999997,999996)'爲例。結果在'1e + 18'左右,但四捨五入到最大。雙精度將給出0. –

+0

@AlixAxel您如何處理負面投入取決於您的要求。如果純粹是爲了提供信息,那麼我只會提出警告。 –

3

是的,這是浮點精度的問題。 -1/2可以精確表示,但1/3和1/6不能。因此,當你添加它們時,你不會得到零。

你可以使用你提到的使用通用分母方法(你發佈的H2和H3公式),但是這只是在一段時間內將罐頭踢出去,一旦總和就會得到不準確的結果產品術語開始四捨五入。

爲什麼要採取可能爲負數的諧波平均值呢?這是一個固有的不穩定計算(H(-2,3,6 + epsilon)對於很小的epsilon變化很大)。

+0

謝謝Keith,關於負數,我只是爲了完整性而已,但我認爲它沒有多大意義。我應該過濾負數還是隻使用它們的絕對值? –

+1

@AlixAxel:我會拋出一個異常,如果你可以在PHP中做到這一點。如果沒有,則返回錯誤代碼。默默地忽略不好的輸入是一個壞主意。 –