2 误差

统计计算的算法要得到正确的结果， 就需要尽可能减少误差。 统计问题中的误差有模型误差、实验误差和数值计算误差， 在统计计算研究中主要解决的是如何减少数值计算误差的问题。

统计计算的算法通常是用来求解某种统计模型。 任何用来解决实际问题的数学模型都或多或少地简化了实际问题， 忽略掉一些细节，从而模型误差不可避免。 如果模型不合适，其它误差控制得再完美， 问题也不能得到解决； 更糟的是， 良好的计算结果会给使用者以错误的信心。 比如，我们使用的回归模型要求观测是独立的， 而实际数据观测有不可忽略的序列相关性， 尽管我们用软件算出了很完美的结果， 这个结果也是错误的。 我们应当仔细选择模型， 尽可能减少模型误差。

建立统计模型所需的数据来自实验、观测、抽样调查等过程， 在这样的过程中会出现实验误差， 包括随机误差、系统误差、过失误差。

随机误差是试验过程中由一系列随机因素引起的不易控制的误差， 可以通过多次重复试验或改进模型设计来减小随机误差。 随机误差可能来自物理量本身的波动， 比如测量风速，就是在测量一个随时变化的物理量， 不可避免地受到随机误差影响。 随机误差可能来自不可控制的随机因素影响， 比如，在用雷达测量飞机的方位和速度时， 可能受到地磁、气温、地形的影响。 由于测量仪器精度的限制也会产生随机误差， 比如用最小刻度是1度的温度计测量温度， 只能把不足1度的值四舍五入或者估计小数点后一位数字。 随机误差也可能来自特定条件下才发生的程序错误。

系统误差是多次测量持续偏高或偏低的误差， 多次重复测量不能消除或减少系统误差。 系统误差可能来自仪器本身的误差， 比如用不锈钢直尺测量家具高度， 直尺本身在温度不同时长度有细微变化。 系统误差也可能来自仪器使用不当， 比如用天平测量质量时天平没有配准。 当发现有系统误差时， 必须找出引起误差的原因并消除。

在记录实验数据时由于人的过失可以导致误差发生, 这样的误差称为过失误差。 比如，在记录仪表（如水表、电表）的读数时看错数字， 在记录数值时写错小数点位置， 在上传数据时报告了过时的或错误的数据，等等。 统计数据分析必须甄别并改正这样的过失误差， 否则会对分析结果产生严重影响。 在使用计算机软件处理数据时程序设计的质量问题也会导致误差发生。 比如，当输入条件不满足模型要求而程序中未进行检查时， 可能给出错误的结果。

数值误差是用电子计算机进行数据存储和计算时产生的误差， 设计算法时必须了解并尽可能避免这种误差。

设\(A\)为结果的精确值，\(a\)为\(A\)在计算机中的近似值，则 \(\Delta = a - A\)称为绝对误差, 简称误差。 \(\delta = \frac{a-A}{A}= \frac{\Delta}{A}\) 称为相对误差。 相对误差没有量纲，常常用百分数表示。 绝对误差和相对误差常常仅考虑其绝对值。 实际中如果能估计绝对误差和相对误差的取值区间或绝对值的上限， 则可以确定误差大小。 如果绝对误差\(\Delta\)的绝对值上限估计为\(\tilde \Delta\), 则相对误差的绝对值上限可以用\(\delta = \frac{\tilde\Delta}{a}\)估计， 因为\(A\)的真实值一般是未知的。 有\(p\)位有效数字的一个十进制近似数的相对误差控制在\(5 \times 10^{-p}\)以下。

没有数值计算经验的读者往往会有一个错误认识， 认为计算机得到的结果总是准确的、可信的。 即使不考虑模型误差、随机误差和系统误差的影响， 用计算机进行数值计算也不可避免地受到误差影响， 只要我们把误差控制在可接受的范围就可以了。

为了解数值计算误差， 我们首先需要了解计算机中数值的存储方式。 计算机中的数都用二进制表示，包括定点表示和浮点表示两种。

定点数主要用于保存整数，是精确表示的， 定点数的四则运算可以得到精确结果， 但整数除法需要特别注意余数问题， 另外定点表示的整数范围有限， 比如用32个二进制位可以表示\(-2^{31} \sim 2^{31}-1\)之间的整数（约10位有效数字）， 定点数的计算结果有可能溢出， 即超出能表示的整数范围。

例2.1 (Julia中的定点整数(*)) 显示Julia语言中整数的界限与二进制表示。

整数1的表示：

上面的x是64位定点整数，首位的0是符号位，表示非负， 后面有62个0和1个1。

上面的x是8位定点整数， 首位的0是符号位，表示非负，后面有6个0和1个1。

上面的x是8位无符号定点整数，有7个0和1个1, 没有符号位。

上面的x是8位定点整数，首位1表示负号，后面有7个1。

Julia中各种整数类型的最大可表示值：

Int8约有2-3位十进制有效数字。

Int16约有4-5位十进制有效数字。

Int32约有9-10位十进制有效数字。

Int64约有18-19位十进制有效数字。

Int64的加法超出存储界限时发生溢出， 而且不提供错误信息：

变成了负数。

但是Int8, Int16, Int32可以自动提升到更大储存的类型:

注意Julia中整数加法、减法、乘法结果还是整数， 除法按实数除法计算， 整数之间的乘法也是整数。 两个整数之间进行乘方计算时结果也是整数， 这是Julia语言的特殊规定， 其它语言中一般不这样规定， R语言中的四则运算和乘方运算基本都是按浮点数计算。

Julia程序：

溢出了。写成实数的乘方：

※※※※※

数值计算时主要使用浮点表示, 比如用64个二进制位能表示绝对值在\(10^{-308} \sim 10^{308}\)之间的实数， 有效数字大约有\(16 \sim 17\)位。 二进制浮点数的表示包括\((S,E,F)\)三个部分， 即正负号、指数部分、小数部分。 小数部分的位数\(d\)决定了精度。一般存储\(2^{-1} \leq F < 1\), 这样\(F\)中首位一般是1，有些计算机不保存这一位。 例如，二进制数\(101.101_2=0.101101_2 \times 2^3\) 可以表示为\((+, +11_2, .101101_2)\)。

不同的计算机硬件可能采用不同的浮点数表示， 在现代常用的IEEE 64位浮点数的表示中， 一个浮点数\(x\)被表示成 \[ x = \pm 2^{q - 1023} \times (1.b_1 b_2 \dots b_{52}) \] 其中符号\(\pm\)用64位的首位表示， 随后11位存储指数部分， 看成一个11位无符号定点整数\(q\)， 指数部分\(E\)代表\(\times 2^{q-1023}\)， 但是\(q=0\)和\(q=2^{11}-1\)有特殊解释。 最后的52位是二进制小数部分\(F\)， 约定小数点前面总是1并且不保存这个1。 当\(q=0\)时，小数部分约定小数点前面是0， 保存所有的二进制小数位， 这时小数部分的有效位数比正常情况少。 \(q > 0\)时绝对值最小的浮点数约为2.2E-308。 当\(q=2^{11}-1\)(所有的11位指数位都为1)时， 代表某些不正常的实数， 这时如果小数位都是零， 用来代表Inf(无穷大); 如果小数位有1， 用来代表NaN(比如0.0/0.0的结果)。 \(q < 2^{11}-1\)时绝对值最大的浮点数约为1.80E308。

由于\(E\)和\(F\)两部分的位数限制， 二进制浮点数只有有限个， 这些数的集合记为\(\mathcal F\)， 在这个集合中可以制定一个次序， \((S,E,F)\)的下一个数是\((S, E, F+2^{-d})\)。 对于实数域\(\mathbb R\)中的实数\(x\)， 当\(|x|\)超出浮点数范围时，我们不能表示； 当\(|x|\)没有超出浮点数范围时， 一般也只能用\(\mathcal F\)中的数来近似表示。 比如，十进制数\(0.1\)如果表示成二进制数， 是一个无限循环小数\(0.000\dot{1}10\dot{0}_2\), 只能存储其中有限位。 在近似表示时对于不能保存的位数计算机中可能用了舍入法或者截断法处理, 这样造成的误差叫做舍入误差。

例2.2 (Julia中的浮点数的二进制表示(*)) 用Julia语言演示64位浮点数\(0.1\)的二进制存储。

这是Julia语言中Float64类型的浮点数，按IEEE 64位浮点数格式存储。 首位0是符号位， 后面11位是指数位，存储\(q\)的无符号整数值，恰好等于1019， 表示\(2^{1019 - 1023} = 2^{-4}\)， 最后52位是小数位， 都在二进制小数点后面， 而小数点前面有一个1， 所以近似表示\(1.\dot 100 \dot 1\)， 因此表示为 \[\begin{aligned} 0.1 \approx & + 2^{-4} \times \left[ 1 + (1 \times 2^{-1} + 0 \times 2^{-2} + 0 \times 2^{-3} + 1 \times 2^{-4}) \right. \\ & \left. + (1 \times 2^{-5} + 0 \times 2^{-6} + 0 \times 2^{-7} + 1 \times 2^{-8}) + \cdots \right] . \end{aligned}\]

※※※※※

设实数\(z\)在浮点数可表示范围内， 将其浮点表示记为\(\mbox{fl}(z)\)， 则绝对舍入误差满足 \[\begin{aligned} |\mbox{fl}(z) - z| \leq U |z|, \forall z. \end{aligned}\] 其中\(U\)称为机器单位(machine unit)， \(\mbox{fl}\)用截断近似时\(U=2^{-(d-1)}\), \(\mbox{fl}\)用四舍五入时\(U = 2^{-d}\)。 一个浮点数\(z\)用浮点表示\(\mbox{fl}(z)\)后误差范围如下 \[\begin{align*} \mbox{fl}(z) = z(1+u), \ |u| \leq U. \end{align*}\]

和\(U\)类似的一个数是机器\(\varepsilon_{\text{m}}\)， \(\varepsilon_{\text{m}}\)是使得\(1+\varepsilon_{\text{m}}\) 和\(1\)的表示不相同的最小正浮点数。 可以认为\(U \approx \varepsilon_{\text{m}}/2\)。 典型的双精度计算\(\varepsilon_{\text{m}}\)约为\(10^{-16}\)数量级， 而单精度的相应值与双精度相差约\(10^9\)倍, 即单精度计算与双精度计算的精度相差约9位有效数字。 使用双精度实数可以减少表示误差和计算误差， 但是不能消除这两种误差， 所以不能盲目相信使用了双精度后就总是能准确地进行数值计算。

在R软件中用变量.Machine$double.eps保存双精度计算的机器\(\varepsilon_{\text{m}}\)值。 在Julia语言中用eps(Float64)获得双精度浮点数的机器\(\varepsilon_{\text{m}}\)值。

由于计算机存储的有限性， 浮点数的四则运算不符合结合律，计算结果与计算次序有关。 不同的算法可能导致不同的计算误差， 应该尽可能选用计算精度高的数学公式， 另外在设计算法时需要注意避免一些损失精度的做法。

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

\(1.1 + 0.1 - 1.2\)的精确结果是0， 但是在R中计算，得

而

两个计算次序不同，结果都不等于零而且不同计算次序结果不同。 Julia计算结果与R的结果相同。

※※※※※

例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} \] 就只要计算一个除法和一个减法。

\(z_2\)定义为 \[\begin{align} z_2 = z_2(x_1, x_2) = \frac{x_1 - \sqrt{x_1^2 - 4 x_2}}{2} \tag{2.1} \end{align}\] 在\(x_2\)变化时条件数 \[\begin{aligned} C =& \left| x_2 \frac{ \partial z_2(x_1, x_2)}{\partial x_2} / z_2(x_1, x_2) \right| \\ =& x_2 \frac{1}{x_1^2 - 4 x_2} \frac{2}{x_1 - \sqrt{x_1^2 - 4 x_2}} = \frac{z_1}{z_1 - z_2} \end{aligned}\] 当\(z_1\)较大，\(z_2\)很小时\(C\)接近1， 问题适定性很好。 当判别式\(x_1^2 - 4 x_2\)很接近于零时两个根\(z_1\)和\(z_2\)差很小， 条件数很大， 这时\(x_2\)的一个微小变化可能引起根\(z_2\)的很大变化。

即使在\(C=1\)时不适当的算法也可能会造成结果不稳定。 设判别式\(x_1^2 - 4 x_2\)不接近于零， 两个根\(z_1\)和\(z_2\)差距很大， 但是当\(x_2\)很小时公式(2.1)的分子中\(x_1\)和\(\sqrt{x_1^2 - 4 x_2}\)很接近， 相减会造成精度损失， 改用公式 \[\begin{align} z_2 = \frac{2}{x_1 + \sqrt{x_1^2 - 4 x_2}} \tag{2.2} \end{align}\] 则避免了相近数相减， 提高了精度。

下面程序中x1>0, x2>0，方程的左边看成z, x1, x2的函数：