下图是sas中proc mi对缺失数据填补的基本原理，本文对这段文字进行解读。

1、采用X1-Xk中非缺失数据建立线性回归模型，得到β0-βk等参数的估计值：β̂ 0-β̂ k，σ̂²为残差的方差。有关β̂ 的估计在我的另一篇文章《回归参数估计和协方差矩阵的推导》中有所详述。

2、通过极大似然估计推导出参数σ̂²的后验分布。其中σ̂²的公式为σ²* = σ̂²(nj-k-1)/g。

该公式原理如下：

(1) 在线性回归中，残差平方和(RSS)除以σ²服从卡方分布：RSS/σ² ~ χ²(nj-k-1)，因为：

         首先，在线性回归中，我们假设误差项ε服从正态分布：ε ~ N(0, σ²)，标准化后的误差项ε/σ ~ N(0, 1)。对于n个独立的标准正态随机变量，它们的平方和服从卡方分布，自由度为n（Z₁² + Z₂² + ... + Zₙ² ~ χ²(n)，其中每个Z都是标准正态分布N(0,1)），而RSS是残差的平方和RSS = Σe² = Σ(y - ŷ)²，故而服从卡方分布。

(2) 最大似然估计的残差方差为：σ̂² = RSS/(nj-k-1)，因为：

        首先，由于误差项服从正态分布，单个观测值yi的概率密度函数为f(yi|xi,β,σ²) = (1/√(2πσ²)) * exp(-(yi - xiᵀβ)²/(2σ²))。

        由于观测值独立，联合概率密度（似然函数）是各个观测值概率密度的乘积（换成人话就是：n个受试者在经过试验后恰好得到现在的观测结果的概率），L(β,σ²|y) = ∏(i=1 to n) f(yi|xi,β,σ²)，由于x和y是已知结果，所以这个函数只有β和σ²是未知数。

        将以上公式展开得：L(β,σ²|y) = (2πσ²)^(-n/2) * exp(-(1/(2σ²)) * ∑(yi - xiᵀβ)²)，，由于其中∑(yi - xiᵀβ)²就是RSS（残差平方和），故而得到L(β,σ²|y) = (2πσ²)^(-n/2) * exp(-RSS/(2σ²))。

        将该公式取对数得到：ln L = -n/2 * ln(2π) - n/2 * ln(σ²) - RSS/(2σ²)。由于极大似然法假设参数得取值为概率最大得结果，所以为求该函数最大值，对两边求导并使得结果为0：

        ∂(ln L)/∂(σ²) = -n/(2σ²) + RSS/(2σ⁴) = 0，进而得到n/σ² = RSS/σ⁴，整理得到σ² = RSS/n。最后由于需要调整自由度，故而将n改为(nj-k-1)，，其中nj为总样本数，k为自变量个数。

        最终得到公式：σ̂² = RSS/(nj-k-1)

(3) 将这两个关系σ̂² = RSS/(nj-k-1)，RSS/σ² ~ χ²(nj-k-1)结合：(nj-k-1)σ̂²/σ² ~ χ²(nj-k-1)

        令 g 为从χ²(nj-k-1)分布中抽取的随机数，则：g = (nj-k-1)σ̂²/σ²。解出σ²：σ² = σ̂²(nj-k-1)/g

3、通过极大似然估计推导出参数β̂ 的后验分布。其中β̂ 的公式为β̂ * = β̂ +σV*(1/2)Z，V*(1/2)表示V的开方，Z为k+1维独立正太分布统计量。

        在我的另一篇文章《回归参数估计和协方差矩阵的推导》中有详述Var(β̂) = σ²(X'X)⁻¹ ，其中V=(X'X)⁻¹。

4、将β和 σ²都换成后验分布的值以后，得到新的线性回归方程：

Y=β̂ *X+Zσ*²

5、采用新的线性回归方程对缺失值进行填补。