Calcolo della media armonica e precisione del galleggiante

Question

Sto implementando lo Pythagorean means in PHP, i metodi aritmetici e geometrici sono un pezzo di torta ma sto avendo davvero difficoltà a trovare un'implementazione affidabile harmonic mean.

Questa è la WolframAlpha definition:

---

E questo è l'implementazione equivalente in PHP:

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

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

    return func_num_args()/$result; 
} 

Ora, se uno degli argomenti è 0 questo sarà lanciare un divisione per avviso 0, ma poiché 1/n corrisponde a n^-1 e pow(0, -1) restituisce con garbo la costante INF senza gettare eventuali errori ho potuto riscrivere che a quanto segue (che sarà ancora gettare errori se non ci sono argomenti, ma lascia ignorare che per ora):

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

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

Entrambe le implementazioni funzionano bene per la maggior parte dei casi (esempio v1, v2 e WolframAlpha), ma non riescono spettacolare se la somma del1/n _iserie è 0, dovrei ottenere un'altra divisione per 0 avvertimento, ma non lo faccio ...

Si consideri il seguente insieme: -2, 3, 6 (WolframAlpha dice che è un complesso infinito):

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

= 0 

Tuttavia, entrambe le implementazioni restituiscono -2.7755575615629E-17 come la somma (v1, v2) anziché 0.

Mentre il risultato ritorno sul CodePad è -108086391056890000 mia macchina dev (32-bit) dice che è -1.0808639105689E+17, ancora è niente come il 0 o INF mi aspettavo. Ho anche provato a chiamare is_infinite() sul valore restituito, ma è tornato come false come previsto.

Ho trovato anche la funzione stats_harmonic_mean() che fa parte dell'estensione stats PECL, ma con mia sorpresa ho ottenuto esattamente lo stesso risultato buggy: -1.0808639105689E+17, se uno degli argomenti è 0, 0 viene restituito, ma nessun controllo sono fatto per la somma della serie, 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 /* }}} */ 

Questo appare come un tipico errore di precisione galleggiante, ma non posso davvero capire il motivo per cui dal momento che i singoli calcoli sono piuttosto precisi:

0.123.

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

È possibile aggirare questo problema senza ripristinare le estensioni gmp/bcmath?

Answer 1

Sei corretto. I numeri che trovi sono un artefatto delle peculiarità dell'aritmetica in virgola mobile.

Aggiungere più precisione non ti aiuterà. Tutto quello che stai facendo è spostare i pali della porta.

La riga di fondo è che i calcoli vengono eseguiti con precisione finita. Ciò significa che ad un certo punto, un risultato intermedio sarà arrotondato. Questo risultato intermedio non è più esatto. L'errore si propaga attraverso i calcoli e alla fine lo fa diventare il risultato finale. Quando il risultato esatto è zero, di solito ottieni un risultato numerico di circa 1e-16 con numeri a precisione doppia.

Questo accade ogni volta che il calcolo comporta una frazione con un denominatore che non è una potenza di 2.

L'unico modo intorno ad esso è quello di esprimere i calcoli in termini di numeri interi o numeri razionali (se è possibile) e utilizzare un pacchetto intero di precisione arbitraria per eseguire i calcoli. Questo è ciò che Wolfram | Alpha fa.

Si noti che anche il calcolo della media geometrica non è banale. Prova una sequenza di 20 volte 1e20. Poiché i numeri sono tutti uguali, il risultato dovrebbe essere 1e20. Ma quello che troverai è che il risultato è infinito. Il motivo è che il prodotto di quei 20 numeri (10e400) non rientra nell'intervallo dei numeri a virgola mobile a precisione doppia e quindi è impostato su infinito. La ventesima radice dell'infinito è ancora l'infinito.

Infine, una meta-osservazione: il Pythogarian significa davvero solo un senso per i numeri positivi. Qual è la media geometrica di 3 e -3? È immaginario ?? La catena di disuguaglianze sulla pagina di Wikipedia a cui si fa riferimento è valida solo se tutti i valori sono positivi.

Answer 2

Sì, questo è un problema con la precisione in virgola mobile. -1/2 può essere rappresentato esattamente, ma 1/3 e 1/6 non possono. Quindi quando li aggiungi, non ottieni assolutamente zero.

Si può fare l'approccio "usa-un-comune denominatore" che hai menzionato (le formule H2 e H3 che hai postato), ma che semplicemente scalfisce il barattolo lungo la strada un po ', otterrai comunque risultati inaccurati una volta la somma termine dei prodotti inizia a arrotondare.

Perché stai prendendo la media armonica di numeri che potrebbero essere negativi, comunque? Questo è un calcolo intrinsecamente instabile (H (-2,3,6 + epsilon) varia notevolmente per epsilon molto piccolo).