高精度计算平均值的最佳策略答案

【问题标题】：Best strategy to compute average with high precision高精度计算平均值的最佳策略
【发布时间】：2016-05-25 19:42:28
【问题描述】：

我正在比较两种计算随机数平均值的算法。

第一个算法将所有数字相加，然后除以最后的项目数
第二种算法计算每次迭代的平均值，并在收到新数据时重复使用结果

我想这里没有什么革命性的东西，而且我不是数学家，所以我不能为这两种算法命名。

这是我的代码：

#include <iostream>
#include <iomanip>
#include <cstdlib>

class Average1
{
public:
    Average1() : total( 0 ), count( 0 ) {}

    void add( double value )
    {
        total += value;
        count++;
    }

    double average()
    {
        return total/count;
    }

private:
    double total;
    size_t count;
};

class Average2
{
public:
    Average2() : av( 0 ), count( 0 ) {}

    void add( double value )
    {
        av = (av*count + value)/(count+1);
        count++;
    }

    double average()
    {
        return av;
    }

private:
    double av;
    size_t count;
};

void compare()
{
    Average1 av1;
    Average2 av2;
    double temp;
    for ( size_t i = 0; i != 100000000; ++i )
    {
        temp = static_cast<double>(std::rand()) / static_cast<double>(RAND_MAX);
        av1.add( temp );
        av2.add( temp );
    }

    std::cout << std::setprecision(20) << av1.average() << std::endl;
    std::cout << std::setprecision(20) << av2.average() << std::endl;
}

int main()
{
    compare();
    return 0;
}

输出是：

0.50001084285722707801
0.50001084285744978875

差异肯定是由于double 类型精度。

到底哪一个是好方法？哪一个给出了真正的数学平均值（或最接近...）？

【问题讨论】：

标签： c++ algorithm average average-precision

【解决方案1】：

如果你真的想要高精度：

考虑任意精度算术（例如使用GMP）
考虑Kahan summation algorithm （可能的编译器问题）
考虑Shewchuk's-algorithm（在Python 中可以使用math.fsum）

编辑： math.fsum 中的 python-docs 也链接到 this Overview of approaches

【讨论】：

【解决方案2】：

我的猜测是第一类给出了更可靠的结果。在第二种情况下，在每次迭代中，由于除以计数，您会进行一些近似，最终所有这些近似加起来就是您看到的结果差异。相反，在第一种情况下，您只需在计算最终除法时进行近似计算。

【讨论】：

不一定正确。如果在某些时候累积和非常接近下一个数字的负数，则添加数字列表可能会损失很多精度。
是的，但是如果 av*count 恰好接近 -value，那么这也适用于第二类的情况。在这种特定情况下，我会选择头等舱，因为它避免了在每次添加新数字时进行除法。但是，我相信也有一些方法可以防止您提到的问题。
是的，第二个也好不到哪里去。经典的稳定算法是newmean = oldmean + (newvalue - oldmean)/newcount。尽管它涉及每次迭代的除法，但它比朴素求和算法更稳定。（我相信与 Kahan 的求和大致相当。）

【解决方案3】：

John D. Cook 给出了他推荐的精彩分析：

av = av + (value - av)/count;

他的帖子以Comparing three methods of computing standard deviation开头。

然后Theoretical explanation for numerical results

最后一个Accurately computing running variance

【讨论】：

【解决方案4】：

我自己的想法是，在除以之前，两者都会计算 count 乘以一个很大的数字，这就解释了为什么你的结果是近似的。我会这样做：

class Average3
{
public:
    Average3() : av( 0 ), count( 0 ) {}

    void add( double value )
    {
        count++;
        av +=  (value - av)/count;
    }
...

但是在添加最后一个数字时，您仍然会失去精度，因为与平均值相比，添加值/计数很小。我很高兴知道我的直觉是否正确

【讨论】：