简介
统计计算在统计建模中扮演着非常关键的角色。虽然我们常常着眼于统计模型、不同的统计推断方法等,但当我们实际将统计理论付诸实践时,我们就需要统计算法及其计算机实现,因此统计计算是统计建模中不可或缺的一环。本书聚焦于如何将统计模型与数据相契合,以及我们用于实现这一契合的统计算法。
统计理论与计算机实现的距离
对数概率
在统计计算中,经常将概率 \(P(E)\) 求对数,原因有二:(1)计算机无法表示极小的数字,(2)求对数可以把乘法转化为加法,而计算机在加法方面要快很多。
表示极小概率
虽然计算机能够完成所有的计算,但其浮点表示法意味着计算机不能以完美的精度表示小数。.Machine$double.xmin
给出可以满足IEEE 754浮点计算技术标准的要求的最小正数:
$double.xmin
.Machine#> [1] 2.225e-308
但事实上,R 能处理的最小正数为 \(5e-324\)。
在统计建模过程中,我们通常会涉及计算某概率分布密度函数,以一元正态分布为例,其密度函数可以写为:\[f(x | \mu,\sigma^2) = \frac{1}{\sqrt{2\pi}\sigma}\exp\{\frac{1}{2\sigma^2}(x-\mu)^2\}.\]给定 \(x\),\(\mu\) 和 \(\sigma\),我们可以通过 dnorm()
函数很容易的计算对应的密度值:
dnorm(1, mean = 0 , sd = 1)
#> [1] 0.242
但因为密度函数的指数形式,密度值会很快变小,以至于计算机无法表示(下溢),如:
dnorm(40, mean = 0 , sd = 1)
#> [1] 0
如果我们做一个对数变换,则可避免计算机无法表示的问题。且对数变换后的数值更容易在计算机中存储:
dnorm(40, mean = 0 , sd = 1, log = TRUE)
#> [1] -800.9
此外在统计计算中,经常会遇见连个密度函数相除,如 \(\frac{f(x)}{g(x)}\),此时在计算过程中很容易遇到计算机上溢或下溢的问题,通过对数变换,\(\exp(\log(f(x)) - \log(g(x)))\) 在数值上会更加稳定。
线性回归模型
数学表达式的形式和在实际计算机计算该表达式的方式可能截然不同,一个非常明显的例子就是线性回归模型的估计。对一个线性模型:\[y = X\beta + \epsilon,\]其中 \(y\) 是 \(n \times 1\) 的响应变量, \(X\) 是 \(n \times p\) 的自变量矩阵,\(\beta\) 是 \(p \times 1\) 的回归系数,\(\epsilon\) 是 \(n \times 1\) 的误差项。
在统计学相关的教科书中我们可以得出该模型参数的最小二乘或者最大似然估计为 \(\hat{\beta} = (X^TX)^{-1}X^Ty\),这并不意味着计算机实现时我们就需要直接计算 \((X^TX)^{-1}X^Ty\),其主要原因是对矩阵 \(X^TX\) 直接求逆在计算上是非常昂贵的,而且如果遇到多重共线性问题,求逆运算在计算上并不稳定。那怎么办呢?