Cálculo del promedio armónico y Precisión del flotador

Question

Implemento el Pythagorean means en PHP, los medios aritméticos y geométricos son sencillos, pero estoy teniendo un momento difícil para encontrar una implementación confiable de harmonic mean.

Este es el WolframAlpha definition:

---

Y esta es la implementación equivalente en PHP:

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

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

    return func_num_args()/$result; 
} 

Ahora, si alguno de los argumentos es 0 esto va a lanzar una división por 0 advertencia, pero desde 1/n es lo mismo que n^-1 y pow(0, -1) regresa con gracia la constante INF sin lanzar cualquier error que pudiera volver a escribir que a la siguiente (que todavía va a lanzar errores si no hay argumentos, pero le permite ignorar que por ahora):

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

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

Ambas implementaciones funcionan bien para la mayoría de los casos (ejemplo v1, v2 y WolframAlpha), pero fallan espectacularmente si la suma de la1/n _iserie es 0, debería tener otra división por 0 advertencia, pero yo no ...

Considere el siguiente conjunto: -2, 3, 6 (WolframAlpha dice que es un complejo infinita):

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

= 0 

Sin embargo, mis dos implementaciones devuelven -2.7755575615629E-17 como la suma (v1, v2) en lugar de 0.

Si bien el resultado retorno de la CodePad es -108086391056890000 mi máquina dev (32 bits) dice que es -1.0808639105689E+17, todavía no es nada como el 0 o INF que estaba esperando. Incluso traté de llamar al is_infinite() en el valor de retorno, pero volvió como false como se esperaba.

También encontré la función stats_harmonic_mean() que es parte de la extensión stats PECL, pero para mi sorpresa que tiene exactamente el mismo resultado con errores: -1.0808639105689E+17, si alguno de los argumentos es 0, 0 se devuelve pero no cheques se hacen a la suma de la 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 /* }}} */ 

Esto parece un error típico flotante de precisión, pero realmente no puedo entender la razón por la ya que los cálculos individuales son bastante precisos:

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

¿Es posible evitar este problema sin tener que volver a las extensiones gmp/bcmath?

Answer 1

Estás en la correcta. Los números que está encontrando son un artefacto de las peculiaridades de la aritmética de coma flotante.

Agregar más precisión no lo ayudará. Todo lo que estás haciendo es mover los mensajes de la meta.

La conclusión es que los cálculos se realizan con precisión finita. Eso significa que en algún punto, se redondeará un resultado intermedio. Ese resultado intermedio ya no es exacto. El error se propaga a través de los cálculos, y finalmente lo convierte en su resultado final. Cuando el resultado exacto es cero, por lo general se obtiene un resultado numérico de alrededor de 1e-16 con números de precisión doble.

Esto ocurre cada vez que el cálculo consiste en una fracción con un denominador que no es una potencia de 2.

La única manera de evitarlo es expresar los cálculos en términos de números enteros o números racionales (si puede) y usa un paquete entero de precisión arbitraria para hacer los cálculos. Esto es lo que hace Wolfram | Alpha.

Tenga en cuenta que el cálculo de la media geométrica tampoco es trivial. Pruebe una secuencia de 20 veces 1e20. Dado que los números son todos iguales, el resultado debe ser 1e20. Pero lo que encontrarás es que el resultado es infinito. La razón es que el producto de esos 20 números (10e400) está fuera del rango de números de coma flotante de precisión doble, por lo que está configurado en infinito. La vigésima raíz del infinito sigue siendo infinita.

Finalmente, una metaobservación: el Pythogarian significa que realmente solo tiene sentido para números positivos. ¿Cuál es la media geométrica de 3 y -3? ¿Es imaginario? La cadena de desigualdades en la página de Wikipedia a la que se vincula solo es válida si todos los valores son positivos.

Answer 2

Sí, este es un problema con la precisión de coma flotante. -1/2 se puede representar exactamente, pero 1/3 y 1/6 no se pueden representar. Por lo tanto, cuando los agrega, no acaba de obtener cero.

Puede hacer el enfoque de denominador común que mencionó (las fórmulas H2 y H3 que publicó), pero eso le da un vuelco a la lata un poco, aún obtendrá resultados inexactos una vez que la suma -el término de los productos comienza a redondearse.

¿Por qué toman la media armónica de los números que podrían ser negativos, de todos modos? Ese es un cálculo inherentemente inestable (H (-2, 3, 6 + epsilon) varía ampliamente para épsilon muy pequeño).