高精度计算平均值的最佳策略
Best strategy to compute average with high precision
我正在比较两种计算随机数平均值的算法。
- 第一个算法将所有数字相加并除以最后的项目数
- 第二种算法计算每次迭代的平均值,并在收到新数据时重新使用结果
我想这里没有什么革命性的东西,而且我不是数学家所以我不能给这两个算法命名。
这是我的代码:
#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
类型精度。
到底哪种方法好呢?哪一个给出了真实的数学平均值(或最接近...)?
我的猜测是第一个 class 给出了更可靠的结果。在第二种情况下,在每次迭代中,由于除以计数,你会做一些近似,最终所有这些近似值加起来就是你看到的结果差异。相反,在第一种情况下,您只是在计算最终除法时进行近似。
如果你真的想要高精度:
- 考虑任意精度算术(例如 GMP)
- 考虑Kahan求和算法
(可能是编译器问题)
- 考虑 Shewchuk's-algorithm (which is available in Python as math.fsum)
编辑: math.fsum 中的 python-docs 也链接到 this Overview of approaches
我自己的想法是,在除以它之前,两者都计算了一个很大的数字,这就解释了为什么你的结果是近似值。我会这样做:
class Average3
{
public:
Average3() : av( 0 ), count( 0 ) {}
void add( double value )
{
count++;
av += (value - av)/count;
}
...
但是在添加最后一个数字时您仍然会失去精度,因为添加 value/count 与平均值相比较小。我很高兴知道我的直觉是否正确
John D. Cook 给出了一个很好的分析,他推荐:
av = av + (value - av)/count;
他的帖子以 Comparing three methods of computing standard deviation 开头。
我正在比较两种计算随机数平均值的算法。
- 第一个算法将所有数字相加并除以最后的项目数
- 第二种算法计算每次迭代的平均值,并在收到新数据时重新使用结果
我想这里没有什么革命性的东西,而且我不是数学家所以我不能给这两个算法命名。
这是我的代码:
#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
类型精度。
到底哪种方法好呢?哪一个给出了真实的数学平均值(或最接近...)?
我的猜测是第一个 class 给出了更可靠的结果。在第二种情况下,在每次迭代中,由于除以计数,你会做一些近似,最终所有这些近似值加起来就是你看到的结果差异。相反,在第一种情况下,您只是在计算最终除法时进行近似。
如果你真的想要高精度:
- 考虑任意精度算术(例如 GMP)
- 考虑Kahan求和算法 (可能是编译器问题)
- 考虑 Shewchuk's-algorithm (which is available in Python as math.fsum)
编辑: math.fsum 中的 python-docs 也链接到 this Overview of approaches
我自己的想法是,在除以它之前,两者都计算了一个很大的数字,这就解释了为什么你的结果是近似值。我会这样做:
class Average3
{
public:
Average3() : av( 0 ), count( 0 ) {}
void add( double value )
{
count++;
av += (value - av)/count;
}
...
但是在添加最后一个数字时您仍然会失去精度,因为添加 value/count 与平均值相比较小。我很高兴知道我的直觉是否正确
John D. Cook 给出了一个很好的分析,他推荐:
av = av + (value - av)/count;
他的帖子以 Comparing three methods of computing standard deviation 开头。