拟合的卡方检验

实验中一个常见的任务是，手头有一组数据，要拟合一条曲线。然后要检验拟合的优度。在使用卡方(\(\chi^2\))或者约化卡方(reduced chi-squares, \(\chi^2_\mathrm{red}\)​​)检验时，会遇到自由度到底等于几的问题。本文先参考[1-2]介绍了测量数据为何服从正态分布，再参考[3]介绍了线性回归的概念和方法，最后参考[4]解释了自由度的问题。

整篇文章不涉及高深的数学知识，也没有数学意义上的严格证明，只有直观解释和物理上的推导，是为理工科实验数据处理而总结的。

测量的物理量的均值

如果我们每次测得的物理量的值服从某均值为\(\mu\)的正态分布，对这样的一组测量结果取均值，视该均值为一个新的随机变量，则其期望是\(\mu\)​，方差是\(S^2/n\)​，其中\(S^2\)​是该组测量结果的样本方差。这是由正态分布的性质决定的。

当物理量的测量值并不服从正态分布时，我们一样可以在\(n\to \infin\)时得到类似结果，推导如下。

仅限定独立同分布，总体的均值为\(\mu\)，方差为\(\sigma^2\)。记

\[z=\frac{\sum_{i=1}^nx_i-n\mu}{\sqrt{n}\sigma} \]

这时中心极限定理给出[2]

\[\lim_{n\to \infin}P\left(z\leqslant z_0\right)=\Phi(z_0) \]

其中\(\Phi(z_0)\)​​为标准正态分布\(N(0,1)\)​​的累积分布函数。换言之，当\(n\)​​很大时，随机变量\(z\)​​趋于标准正态分布\(N(0,1)\)​​，即

\[\bar{x}\equiv\frac{1}{n}\sum_{i=1}^nx_i\sim N(\mu,\sigma^2/n), n\to\infin \]

如果在上式右侧用样本方差\(S^2\)代替总体方差\(\sigma^2\)，用样本均值\(\bar{x}\)代替总体均值\(\mu\) (注意左侧的\(\bar{x}\)理解为新构建的一个随机变量，右侧的\(\bar{x}\)理解为某一组测量的均值，它作为总体均值\(\mu\)的一个估计)，则可以说一组重复测量(\(n\gg 1\))的平均值，作为一个新的随机变量，它满足正态分布，均值为\(\bar{x}\)，方差为\(S^2/n\)。

线性拟合

现在考虑线性拟合，先简述线性拟合的概念。包含\(P\)个待拟合参数\(b_i\)的线性模型表达式如下：

\[f(\mathbf{x};\mathbf{b})=\sum_{i=1}^Pb_ig_i(\mathbf{x}) \]

其中\(g_i(\mathbf{x})\)是\(P\)个线性无关的连续函数。注意线性拟合只要求表达式关于待拟合参数是线性的，而不要求关于\(\mathbf{x}\)是线性的。一个线性模型的例子为

\[f(\mathbf{x};\mathbf{b})=b_1x_1+b_2x_2^2+b_3x_3^3 \]

而非线性模型的例子为

\[f(\mathbf{x};\mathbf{b})=b_3x_2\cos(b_1x_1+b_2) \]

现有一组数据点\((\mathbf{x}_i,y_i),i=1,2,\cdots,N\)，其中\(y_i\)是在条件\(\mathbf{x}_i\)下独立多次重复测量取均值的结果，且\(y_i\)有对应的误差\(\sigma_i\)，该误差\(\sigma_i\)就是上一节所说的正态分布的标准差，也即该组重复测量结果的标准误。

我们有数据点在拟合时，要最小化如下目标函数\(\chi^2\)，待优化的参数是\(\mathbf{b}\)​

\[\chi^2=\sum_{i=1}^N\left(\frac{y_i-f(\mathbf{x}_i;\mathbf{b})}{\sigma_i}\right)^2 \]

其中线性模型\(f\)的表达式如前。优化过程就是求解方程组

\[\frac{\partial \chi^2}{\partial b_i}=0, i=1,\cdots,P\tag{$*$} \]

根据(9)和(12)式，可知线性模型的优点就在于\(\chi^2\)​关于诸\(b_i\)​是二次的，求导后得到一组线性方程组，可以直接求解。为了将此事说清楚，记

\[\mathbf{y}=(y_1,y_2,\cdots,y_N)^T\\ \mathbf{b}=(b_1,b_2,\cdots,b_P)^T\\ \mathbf{X}_{n,p}=g_p(\mathbf{x}_n)\\ \mathbf{\Sigma}=\mathrm{diag}\{\sigma_1^2,\sigma_2^2,\cdots,\sigma_N^2\} \]

其中\(\mathbf{y}\)​是\(N\times1\)​列向量，\(\mathbf{b}\)​是\(P\times 1\)​的列向量，\(\mathbf{X}\)​是\(N\times P\)​矩阵（通常\(N\gg P\)​即数据点数目远远多于待拟合参数的数目），\(\mathbf{\Sigma}\)​是\(N\times N\)​的测量误差的协方差矩阵，因为独立测量，协方差矩阵对角。利用上面定义的矩阵，可以将\(\chi^2\)​改写如下

\[\chi^2=(\mathbf{y}-\mathbf{X}\mathbf{b})^T\mathbf{\Sigma}^{-1}(\mathbf{y}-\mathbf{X}\mathbf{b}) \]

而求解线性方程组得到的解表示为(加帽表示估计值)

\[\hat{\mathbf{b}}=(\mathbf{X}^T\mathbf{\Sigma}^{-1}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Sigma}^{-1}\mathbf{y} \]

于是\(\mathbf{y}\)的预测值\(\hat{\mathbf{y}}\)可以表示为

\[\hat{\mathbf{y}}=\mathbf{X}\hat{\mathbf{b}}=\mathbf{X}(\mathbf{X}^T\mathbf{\Sigma}^{-1}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Sigma}^{-1}\mathbf{y}=\mathbf{H}\mathbf{y} \]

其中

\[\mathbf{H}=\mathbf{X}(\mathbf{X}^T\mathbf{\Sigma}^{-1}\mathbf{X})^{-1}\mathbf{X}^T\mathbf{\Sigma}^{-1} \]

它把测量的\(\mathbf{y}\)​​转换成预测的\(\hat{\mathbf{y}}\)​​. 拟合参数的有效数目等于\(\mathrm{tr}(\mathbf{H})\)​​即\(\mathbf{H}\)​​的迹，它也等于\(\mathbf{X}\)​​的秩[3]。前面提到，我们限制诸\(g_i(x)\)​​线性无关，就是为了避免\(\mathbf{X}\)​​非满秩（注意到\(N\gg P\)​​），或者换句话说，是为了保证求解的线性方程组\((*)\)​式中的方程都是线性独立的。

在下面讨论中，始终默认诸\(g_i(x)\)​​线性无关。

线性拟合的\(\chi^2\)自由度

拟合的过程就是最小化\(\chi^2\)​​​的过程。其物理意义是，在考虑到测量误差的情形下，正确的参数值应该给出最小的残差。我们令

\[z_i=\frac{y_i-f(\mathbf{x}_i;\mathbf{b})}{\sigma_i} \]

显然\(z_i\sim N(0,1)\)​，且\(\chi^2\)​可改写为\(\chi^2=\sum_{i=1}^N z_i^2\)​​。根据统计学常识，\(n\)个服从标准正态分布的随机变量的平方和服从自由度为\(n\)的\(\chi^2\)​分布。但事实是，如果(18)式中的\(\mathbf{b}\)是最小二乘回归给出的，则自由度并非\(N\)，而是\(N-P\)​​​。下面将给出一个直观解释和一个数学推导。

直观解释

考虑3个数据点\((x_i,y_i),i=1,2,3\)的直线拟合\(f(x)=kx+b\)，待拟合参数有\(2\)个. 作为最简单的尝试，我们连接两点\((x_1,y_1)\)和\((x_3,y_3)\)得到一条直线作为拟合直线，然后使用\((x_2,y_2)\)去评估拟合得好不好。于是在\(\chi^2\)的计算中， 只有点\((x_2,y_2)\)对求和有贡献，因为另外两个点都在直线上，贡献为零。这时候，很显然\(\chi^2\)的自由度为\(1\)，正好等于\(3-2\)​.

数学推导

为了稍微严格地解释这一问题（物理上的严格），我们需要从\(\chi^2\)​分布的定义出发[4]。先讨论另一个问题，即前面定义的统计量\(\chi^2\)​服从何种概率分布，不妨记为\(P(\chi^2)\)​​​。同样，​记

\[z_i=\frac{y_i-f(\mathbf{x}_i;\mathbf{b})}{\sigma_i} \]

显然根据如前所述，设有\(N\)​​​个独立的\(z_i\sim N(0,1)\)​​​​​，\(\chi^2\)​​​改写为\(\chi^2=\sum_{i=1}^N z_i^2\)​​​. 考虑\(N\)​​​维空间中的随机点\((z_1,z_2,\cdots,z_N)\)​​​，该随机点的概率密度显然正比于\(N\)​​​个独立的\(z_i\)​​​的概率密度的乘积，即正比于\(\exp(-\chi^2/2)\)​​​​​​。我们现在要求\(\chi^2\)​​​的概率密度，这可以通过两步实现：先通过\(N\)​​​维空间中随机点的概率密度求出随机变量\(\chi\)​​​的概率密度\(P(\chi)\)​​​，然后转换为\(\chi^2\)​​​的概率密度\(P(\chi^2)\)​​​。考虑这\(N\)​​​维空间中所有\(\chi\)​​​相等的随机点构成的曲面，显然是一个半径为\(\chi\)​​的\(N\)​​维球面，该球面的表面积显然正比于\(\chi^{N-1}\)​​。于是随机点落在该半径为\(\chi\)​​、厚度为\(\mathrm{d}\chi\)​​的球壳上的概率，就正比于球壳的体积乘以随机点的概率密度，即\(\exp(-\chi^2/2)\chi^{N-1}\mathrm{d}\chi\)​​，另一方面，该概率也等于随机变量\(\chi\)​​落在区间\(\chi\)​​到\(\chi+\mathrm{d}\chi\)​​上的概率\(P(\chi)\mathrm{d}\chi\)​​​​。于是

\[P(\chi)\propto \exp(-\chi^2/2)\chi^{N-1} \]

从而转换到\(\chi^2\)的概率密度为

\[P(\chi^2)\propto \exp(-\chi^2/2)\chi^{N-2} \]

至于比例常数，只需将\(P(\chi^2)\)归一化即可算出，它和\(N\)有关，结果为

\[P(\chi^2,N)=\frac{2^{-N/2}}{\Gamma(N/2)}\chi^{N-2}\mathrm{e}^{-\chi^2/2} \]

上式就是所谓的\(\chi^2\)分布。

如果现在\(N\)个\(z_i\)依然服从\(N(0,1)\)，但相互不独立，而是满足如下线性约束，

\[\sum_{i=1}^N C_iz_i=0 \]

现在重新考虑\(\chi^2\)​的概率分布。显然，这时随机点将位于原先的\(N\)​维空间的一个\(N-1\)​​维子空间中，剩余的推导照旧，可得出现在的\(\chi^2\)​服从自由度为\(N-1\)​的\(\chi^2\)​分布\(P(\chi^2,N-1)\)​。推而广之，如果我们有\(P\)​个不线性相关的线性约束

\[\sum_{i=1}^N C^{(n)}_iz_i=0,n=1,\cdots,P \]

则\(\chi^2\)服从自由度为\(N-P\)的\(\chi^2\)​分布。

总结起来，如果参数\(\hat{\mathbf{b}}\)​​是通过拟合得出，那么最终的\(\chi^2|_{\mathbf{b}=\hat{\mathbf{b}}}\)​​的自由度为\(N-P\)​​，因为这时按\((*)\)​式，存在着\(P\)​个线性约束；如果参数\(\mathbf{b}\)​是通过其他不依赖于数据集\((\mathbf{x}_i,y_i)\)​的方法得出的（比如来自文献经验值），那么通过数据集计算得出的\(\chi^2\)​的自由度就是\(N\)​​。​

对于某些情况，例如非线性模型，就没有形如\((*)\)​​式的线性约束，这时候“自由度”这一概念严格来说是不存在的。对此，文献[4]认为只要模型的非线性行为并不夸张，人们还是可以勉强使用“自由度”这一概念。但是文献[3]坚持对于非线性模型，并不存在“自由度”这一概念。并且文献[3]还认为，带参数上下界约束的线性模型本质上仍然是非线性的，也没有“自由度”这一概念，对此我的理解是只要最后求出的待拟合参数没有落到约束边界上，那么这个模型仍然可以认为是线性的。

拟合优度的\(\chi^2\)​​检验

如前所述，如果\(\chi^2|_{\mathbf{b}=\hat{\mathbf{b}}}\)​太大，则说明拟合太差；如果拟合很好，那么计算出的\(\chi^2|_{\mathbf{b}=\hat{\mathbf{b}}}\)​应该很好得服从相应自由度的\(\chi^2\)​分布，即：不会太大，也不会太小。那么问题就成为，怎样的\(\chi^2\)​​值算大呢？可以用\(\chi^2\)​概率来衡量[4]，

\[\mathrm{Prob}(\chi^2,K)=\int_{\chi^2}^\infty\mathrm{d}\chi^2P(\chi^2,K) \]

上式的物理意义是，一个真正服从相应自由度的\(\chi^2\)​分布的随机变量，大于待检验的\(\chi^2\)的值的概率。如果这个值很小，比方说\(0.001\)，就意味着如果拟合是正确的那么误差未免太大了，误差达到这么大的概率只有\(0.001\)​。那如何判断\(\chi^2\)值是否太小呢？结论是没法判断，因为\(\chi^2\)检验本身是单边检验，拒绝域仅位于一侧。当然，如果计算出来的\(\chi^2\)确实非常小，比如\(\mathrm{Prob}(\chi^2,K)\approx 1\)，那么可能是测量的\(\sigma_i\)​搞错了，或者确实就是个小概率事件。

有时也会使用约化卡方(reduced chi-squares)[5]，即\(\chi^2_\mathrm{red}\equiv\chi^2/K\)​，它等于\(\chi^2\)​除以其自由度。可以验证，自由度为\(K\)​的\(\chi^2\)​分布，其均值\(E(\chi^2)=K\)​，方差\(D(\chi^2)=2K\)​。于是\(E(\chi^2_\mathrm{red})=1\)​，方差\(D(\chi^2_\mathrm{red})=2/K\)​，当自由度\(K\)​很大时，\(\chi^2\)​分布趋于高斯分布，于是有\(\chi^2_\mathrm{red}\sim N(1,\sigma^2=2/K),K\gg1\)​. 这时约化卡方趋于一个以\(1\)​​为中心的较窄的高斯分布，可以检验拟合的好坏。​但仍需注意，尽管约化卡方趋于高斯分布两边对称看上去好像可以用来做双边检验，但仍然和卡方一样，是单边检验。

参考文献

[1] 陈希孺. 概率论与数理统计. P90-91.

[2] 同上, P133-134.

[3] Andrae. Dos and don'ts of reduced chi-squares (2010).

[4] Statistics - A guide and reference to the use of statistical methods in the physical sciences, chapter 6 and 8 (1993).

[5] Data reduction and error analysis for physical sciences, chapter 11 (2003)