Поддержание точности с плавающей точкой с текущим средним значением

Мне нужно вычислить среднеквадратичную ошибку 16-разрядной операции для произвольного количества точек данных (более 100 миллионов). Я решил пойти со средним значением, поэтому мне не пришлось бы беспокоиться о переполнении из добавления большого числа квадратов ошибок. У 100 миллионов образцов у меня были проблемы с точностью с плавающей запятой (неточные результаты), поэтому я переместился в double.

Вот мой код

int iDifference = getIdeal() - getValue();

m_iCycles++;


// calculate the running MSE as

// http://en.wikipedia.org/wiki/Moving_average

// MSE(i + 1) = MSE(i) + (E^2 - MSE(i))/(i + 1)

m_dMSE = m_dMSE + ((pow((double)iDifference,2) - m_dMSE) / (double)m_iCycles);

Есть ли лучший способ реализовать это для поддержания точности? Я считал нормализацию MSE одним и просто держал сумму с окончательным делением на завершение, чтобы вычислить среднее значение.

Ответ 1

Числа с плавающей запятой не переполняются в подобной ситуации, они только теряют точность. Таким образом, здесь нет преимуществ от среднего значения по текущей сумме. Следствием является то же, растет ли текущая сумма или знаменатель.

Чтобы поддерживать точность в общей сумме, сохраняйте промежуточные итоги вместо одного итога. Просто продолжайте добавлять к промежуточному итогу, пока добавление еще одного приведет к переполнению. Затем переходите к следующему промежуточному итогу. Поскольку они все одинаковы по порядку величины (в базе 2), оптимальная точность может быть достигнута путем преобразования в плавающую точку и использования парного накопления в одну конечную сумму.

// first = errors, second = counter
typedef pair< vector< uint32_t >, uint32_t > running_subtotals;

void accumulate_error( uint32_t error, running_subtotals &acc ) {
    ( numeric_limits< uint32_t >::max() - error < acc.first.back()?
        * acc.first.insert( acc.first.end(), 0 ) : acc.first.back() )
        += error; // add error to current subtotal, or new one if needed
    ++ acc.second; // increment counter
}

double get_average_error( running_subtotals const &total ) {
    vector< double > acc( total.first.begin(), total.first.end() );
    while ( acc.size() != 1 ) {
        if ( acc.size() % 2 ) acc.push_back( 0 );
        for ( size_t index = 0; index < acc.size() / 2; ++ index ) {
            acc[ index ] = acc[ index * 2 ] + acc[ index * 2 + 1 ];
        }
        acc.resize( acc.size() / 2 );
    }
    return acc.front() / total.second;
}

Ответ 2

Возможно, вы захотите посмотреть Kahan Summation Algorithm - это не совсем то, что вам нужно здесь, но оно решает очень похожую проблему, и вы может быть в состоянии адаптировать его к вашим потребностям.

Ответ 3

Если ваши другие решения не работают, вы можете исследовать библиотеку Bignum

"GMP является бесплатной библиотекой для произвольной арифметики точности, работающей со знаком целых чисел, рациональных чисел и чисел с плавающей запятой. Практически не предел точности, кроме тех, которые подразумеваются доступной памятью в машине GMP. GMP имеет богатый набор функций, а функции имеют обычный интерфейс".

Ответ 4

То, что у вас есть, является экспоненциальной скользящей средней. Это более тяжело переносит ошибки в данных, чем раньше. Вам нужно линейное среднее. Средние ваши данные в блоках, скажем, 1 миллион, затем возьмите среднее из этих блоков. Вы даже можете сделать это на нескольких уровнях. Это будет равномерно взвешивать все точки ошибки.