理解计算机中的数字

谢益辉 2008-10-04

很多人对统计计算或数值计算不以为然,表现在他们充分相信计算机、在遇到计算问题时就是一句话,“交给软件去做就可以了”。不过我觉得稍微了解一下数字的常识还是很有必要的,一些表达式在数学上成立,在计算上未必成立。据R-FAQ介绍,在"The Elements of Programming Style"一书中有这么一句话:

10.0 times 0.1 is hardly ever 1.0. (10.0乘以0.1通常不是1.0。)

类似地,在R里面还有这样的例子:

.17 - .6 + .43
## 5.551115e-17
.17 + .43 - .6
## 0
.6 - .17 - .43
## -5.551115e-17

奇怪么?不奇怪(而且注意这不是R的问题),想象与现实总是有差距。以前我见有人用R求极大似然估计,直接使用了密度函数值相乘再最大化,却不知成百上千个小数乘起来在计算机中会是怎样难以表达,我换作取对数求和的形式之后,估计就稳定多了。编程的人(尤其是统计编程的人)千万莫偷懒,写代码之前一定要想好把数学问题转换一下,让计算机能够适应你要用的数字。还有一个最简单的例子就是求组合数,理论上可以用排列数(函数factorial())乘乘除除得出来,但是当n取好几百的大数时,想想n的排列是多大的数、计算机能否表达?这种情况下就应该用组合数本身的函数choose()——理论上等价,计算上不等价。万一choose()本身也很难计算了,那只能考虑先求对数再求和最后求幂,这样也许损失精度,但计算上可实现。

最近R-help里面出现过好几例类似的案子,甚至有人在优化时初始值取的是Inf(无穷大),让人哭笑不得、无可奈何。昨天又有人问1 - cumsum(rep(0.1, 10))的最后一个数字为什么不是0。