编程计算机器误差的值,2 误差 | 统计计算

2.2.2 计算误差

由于计算机存储的有限性，

浮点数的四则运算不符合结合律，计算结果与计算次序有关。

不同的算法可能导致不同的计算误差，

应该尽可能选用计算精度高的数学公式，

另外在设计算法时需要注意避免一些损失精度的做法。

例2.3 (结合律的反例)浮点数的四则运算不符合结合律，计算结果与计算次序有关。

\(1.1 + 0.1 - 1.2\)的精确结果是0，

但是在R中计算，得

而

两个计算次序不同，结果都不等于零而且不同计算次序结果不同。

Julia计算结果相同。

※※※※※

例2.4 (累加的简化例子)

计算

\[

\sum_{n=1}^{999} \frac{1}{n(n+1)}

\]

可以直接累加计算，造成很大的累积误差。 只要把公式变换成

\[\begin{aligned}

\sum_{n=1}^{999} \frac{1}{n(n+1)}

= \sum_{n=1}^{999}\left( \frac{1}{n} - \frac{1}{n+1} \right)

= 1 - \frac{1}{1000} = 0.999

\end{aligned}

\]

就只要计算一个除法和一个减法。

在R中计算，

后一算法与精确值0.999的差等于0，

前一算法与精确值0.999的差等于6.66133814775094e-16。

Julia中计算：

累加计算的结果有\(10^{-16}\)级别的误差。

※※※※※

在进行浮点数加减时，两个绝对值相差很大的数的加减严重损失精度。

设\(|a| \gg |b|\)(\(\gg\)表示“远大于”), 则\(a+b\)会被舍入为\(a\)。

比如，计算

\[\begin{align*}

0.1234 + 0.00001234

\end{align*}\]

如果仅允许保留4位有效数字，则结果为\(0.1234\)。

为避免这样的问题，应该只对大小相近的数相加减，

比如，可以把小的数组合后分别相加，再加起来。

两个相近数相减损失有效数字。如:

\[\begin{align*}

0.8522 - 0.8511 = 0.0011

\end{align*}\]

有效数字个数从4位减低到2位。

统计中如下的方差计算公式：

\[\begin{align*}

\sum_{i=1}^n (x_i - \bar x)^2 = \sum_{i=1}^n x_i^2 - \frac{1}{n} (\sum_{i=1}^n x_i)^2

\end{align*}\]

就存在这样的问题。

在对分布函数\(F(x)\)计算右尾部概率\(1 - F(x)\)时，

如果\(x\)很大则\(F(x)\)很接近1，计算\(1 - F(x)\)会造成结果有效数字损失。

为此，统计软件中计算分布函数时常常可以直接计算右尾概率，

以及概率的对数。

例2.5

逻辑斯蒂分布函数为\(F(x) = \frac{1}{1 + e^{-x}}\),

当\(x\)很大时如果直接计算

\(1 - \frac{1}{1 + e^{-x}}\)，就会出现有效位数损失； 只要改用

\[

\begin{aligned}

1 - F(x) = \frac{1}{1 + e^x}

\end{aligned}

\]

计算，就可以避免两个相近数相减，提高精度。

例如，

用\(1 - \frac{1}{1 + e^{-x}}\)计算\(1 - F(8)\)结果为:

用\(\frac{1}{1 + e^x}\)计算\(1 - F(8)\)结果为

第一种计算公式损失了4位以上的有效数字:

用Julia语言进行上述计算：

2-element Array{String,1}:

"0.00033535013046637197"

"0.0003353501304664781"

第一种计算公式损失了4位以上的有效数字。

※※※※※

例2.6计算\(\sqrt{x^2 + 1} - |x|\)，

直接计算在\(|x|\)很大时会有较大误差，

改成\(1/(\sqrt{x^2 + 1} + |x|)\)则不损失精度。

例2.7 (数值微分中的浮点运算误差)

如果某函数参数相对变动在机器精度U附近时函数值也有变化，

则函数值的计算结果被舍入误差严重影响，不能得到有效计算结果。

对一元可微函数\(f(x)\),

\[f'(x) = \lim_{h\to 0} \frac{f(x+h) - f(x)}{h}\]

当\(h\)很小时可以用\(\frac{f(x+h) - f(x)}{h}\)近似计算\(f'(x)\)。

取\(f(x) = e^x\), 用这种方法近似计算\(x=0\)处\(f'(x)\)的值。真值等于1。

令\(g(h) = \frac{e^h - 1}{h} - 1\)为近似计算的误差。

用浮点数计算，分别取\(h=10^{-1}, 10^{-2}, \dots, 10^{-20}\)。

发现\(g(0.1) = 0.052\)有较大误差，逐步减小\(h\)的值， 同时误差变小，

\(g(10^{-8}) = -6.1 \times 10^{-9}\)， 这时误差已经比较小。

为了进一步减小误差， 取更小的\(h\)值，发现并不能减小误差，

\(g(10^{-11}) = 8.3\times 10^{-8}\)与\(h=10^{-8}\)时误差相近。

继续减小\(h\)的值，发现\(g(10^{-15}) = 0.11\)， 误差反而变得很大，

这是因为\(\exp(10^{-15}) - 1\)是两个十分接近的数相减，

有效位数损失很多，除以\(10^{-15}\)又放大了误差。

进一步减小\(h\)，\(g(10^{-16})=-1\),

这是因为\(\exp(10^{-16}) - 1\)的浮点结果等于0。

此例的Julia语言版本：

20×2 Array{Float64,2}:

0.1 0.0517092

0.01 0.00501671

0.001 0.000500167

0.0001 5.00017e-5

1.0e-5 5.00001e-6

1.0e-6 4.99962e-7

1.0e-7 4.94337e-8

1.0e-8 -6.07747e-9

1.0e-9 8.27404e-8

1.0e-10 8.27404e-8

1.0e-11 8.27404e-8

1.0e-12 8.89006e-5

1.0e-13 -0.000799278

1.0e-14 -0.000799278

1.0e-15 0.110223

1.0e-16 -1.0

1.0e-17 -1.0

1.0e-18 -1.0

1.0e-19 -1.0

1.0e-20 -1.0

※※※※※

例2.8误差的累积会放大误差。

把0.1累加一千万次，与真实值1e6比较。

误差约\(10^{-9}\)。

注意双精度数\(\varepsilon_m\)在\(10^{-16}\)左右，

1e6的表示误差在\(10^{-10}\)左右，

这个累加误差比表示误差大了一个数量级。

Julia语言版本：

误差约\(10^{-4}\)。

注意双精度数\(\varepsilon_m\)在\(10^{-16}\)左右，

1e6的表示误差在\(10^{-10}\)左右，

这个累加误差比表示误差大了6个数量级。

Julia语言结果与R语言的计算结果不同。

※※※※※

浮点数做乘法运算时，结果有效位主要取决于精度较低的数。

做除法运算时，如果分母绝对值很小则会产生较大的绝对误差。

在比较两个浮点数是否相等时，

因为浮点数表示和计算的不精确性，

程序中不应该直接判断两个浮点数是否相等，

而应该判断两者差的绝对值是否足够小。

设计计算机数值计算的算法时一定要注意浮点算法的局限。

要了解能保存的最大和最小实数。例如，浮点数的范围有时不够用。

超过最大值会发生向上溢出,

绝对值小于最小正浮点数会发生向下溢出。

向上溢出一般作为错误；向下溢出经常作为结果零。

还要了解相对误差的范围，

用机器精度\(U\)或\(\varepsilon_m\)表示。

例2.9 (计算同生日概率)

设一个班有\(n\)个人，

则至少有两人的生日的月、日重合的概率为

\[\begin{aligned}

p_n = 1 - \frac{P_{365}^n}{365^n}

\end{aligned}\]

其中\(P_{365}^n = 365\times 364 \times\cdots \times (365-n+1)\)。

如果直接计算分式的分子和分母， 则当\(n\)达到121时分母的计算溢出。

只要改用如下计算次序就可以解决这个问题：

\[\begin{aligned}

P_n = 1 - \prod_{j=0}^{n-1} \frac{365-j}{365}

\end{aligned}\]

事实上，以上的对\(n\)的循环还可以用cumprod()函数简化为：

Julia语言版本：

365-element Array{Float64,1}:

0.0

0.00273973

0.00820417

0.0163559

0.0271356

0.0404625

0.0562357

0.0743353

0.0946238

0.116948

0.141141

0.167025

0.19441

⋮

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

1.0

注意Julia语言的循环不损失效率，

不需要费尽心思改成向量化的程序风格；

而R如果用循环则会慢一到二个数量级。

※※※※※