Расчет среднего гармонического и точность с плавающей точкой

Я внедряю пифагорейские средние в PHP, арифметические и геометрические средние несложны, но мне действительно трудно придумать надежную реализацию гармонического среднего.

Это определение WolframAlpha:

Определение среднего гармонического от 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 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 расширения?

Источник

2 ответа

Решение

Ты прав. Числа, которые вы находите, являются артефактом особенностей арифметики с плавающей точкой.

Добавление большей точности вам не поможет. Все, что вы делаете, это перемещаете посты ворот.

Суть в том, что расчеты производятся с конечной точностью. Это означает, что в какой-то момент промежуточный результат будет округлен. Этот промежуточный результат тогда уже не точен. Ошибка распространяется через вычисления и в конечном итоге превращается в ваш конечный результат. Когда точный результат равен нулю, вы обычно получаете числовой результат около 1e-16 с числами двойной точности.

Это происходит каждый раз, когда ваш расчет включает дробь со знаменателем, который не является степенью 2.

Единственный способ обойти это - выразить вычисления в виде целых или рациональных чисел (если вы можете) и использовать целочисленный пакет произвольной точности для выполнения вычислений. Это то, что делает Wolfram|Alpha.

Обратите внимание, что вычисление среднего геометрического также не тривиально. Попробуйте последовательность из 20 раз 1e20. Поскольку числа все одинаковые, результат должен быть 1e20. Но вы обнаружите, что результат бесконечен. Причина в том, что произведение этих 20 чисел (10e400) находится за пределами диапазона чисел с плавающей запятой двойной точности, и поэтому оно установлено в бесконечность. 20-й корень бесконечности все еще бесконечен.

Наконец, мета-наблюдение: средства Пифогаря действительно имеют смысл только для положительных чисел. Что такое среднее геометрическое 3 и -3? Это мнимое?? Цепочка неравенств на странице Википедии, на которую вы ссылаетесь, действительна, только если все значения положительные.

Источник

Да, это проблема точности с плавающей запятой. -1/2 может быть представлено точно, но 1/3 и 1/6 не могут. Таким образом, когда вы складываете их, вы не получаете ноль.

Вы можете использовать упомянутый вами подход "использование общего знаменателя" (формулы H2 и H3, которые вы разместили), но это просто пойдет по пути, и вы получите неточные результаты, как только сумма продуктов срок начинает округляться.

Почему вы берете гармоническое среднее чисел, которые в любом случае могут быть отрицательными? Это изначально нестабильный расчет (H(-2,3,6+ эпсилон) сильно варьируется для очень маленьких эпсилон).

Источник
Другие вопросы по тегам