33、环境数据分析中的统计检验方法

环境数据分析中的统计检验方法

1. 引言

在环境数据分析中，统计检验是评估数据质量、模型拟合效果以及判断数据特征显著性的重要手段。本文将介绍多种统计检验方法，包括粒子大小数据的统计检验、广义最小二乘法的卡方检验、拟合改进的检验、频谱峰值显著性检验以及相关时间序列频谱峰值的显著性评估。

2. 粒子大小数据的统计检验

2.1 MATLAB 和 Python 代码实现

在 MATLAB 和 Python 中，对于粒子大小数据的统计检验有如下代码实现：

MATLAB 代码：

% edama_12_07, statistical tests of particle size data
...
if( F<1 )
    F=1/F;
end
PF = 1 - (fcdf(F,NA,NB)-fcdf(1/F,NA,NB));

Python 代码：

# edapy_12_07: statistical tests of particle size data
...
if( F<1 ):
    F=1/F;
PF = 1 - (stats.f.cdf(F,NA,NB)-stats.f.cdf(1/F,NA,NB));

2.2 仪器相关问题的检验结果

对于仪器的相关问题，通过统计检验得到以下结果：
| 问题 | 答案 |
| — | — |
| 校准是否正确？ | 我们不能以 95% 的确定性拒绝原假设，即估计校准值与制造商声明校准值之间的差异是由随机变化引起的。 |
| 方差是否在规格范围内？ | 我们不能以 95% 的确定性拒绝原假设，即估计方差与制造商声明方差之间的差异是由随机变化引起的。 |
| 两次测试之间校准是否发生变化？ | 我们不能以 95% 的确定性拒绝原假设，即两次校准之间的差异是由随机变化引起的。 |
| 测试之间方差是否发生变化？ | 我们不能以 95% 的确定性拒绝原假设，即两次估计方差之间的差异是由随机变化引起的。 |

在每种情况下，答案都是基于特定的原假设给出的，且并非确定性的回答，而是基于概率。

3. 广义最小二乘法的卡方检验

3.1 原假设

原假设为：广义误差与其期望值的差异是由于随机变化引起的。

3.2 广义误差的组成

广义误差 $E_{est}^T$ 由两部分组成：
$E_{est}^T = E_{est} + E_{est}^p$

其中，第一部分是数据的估计误差 $E_{est}$：
$E_{est} = \sum_{i=1}^{N} \frac{(d_{obs}^i - d_{pre}^i)^2}{\sigma_{d_i}^2}$
其中 $d_{pre} = Gm_{est}$

第二部分是先验信息的估计误差 $E_{est}^p$：
$E_{est}^p = \sum_{i=1}^{K} \frac{(h_{pri}^i - h_{pre}^i)^2}{\sigma_{h_i}^2}$
其中 $h_{pri} = Hm_{est}$

数据向量 $d_{pre}$ 的长度为 $N$，先验信息向量 $h_{pri}$ 的长度为 $K$，模型参数向量 $m$ 的长度为 $M$。数据的方差 $\sigma_{d_i}^2$ 和先验信息的方差 $\sigma_{h_i}^2$ 是先验方差。

3.3 总误差的分布和自由度

总误差 $E_{est}^T$ 是正态分布随机变量的平方和，因此服从卡方分布。自由度 $\nu_T$ 为 $(N + K) - M$。总误差的均值为 $\nu_T$，方差为 $2\nu_T$，95% 置信区间近似为：
$\nu_T - 2\sqrt{2\nu_T} < E_{est}^T < \nu_T + 2\sqrt{2\nu_T}$

当广义最小二乘法误差落在该区间之外时，与先验方差不兼容，原假设不太可能成立。误差在区间右侧对应模型拟合数据比随机变化预期的更差，可能是模型过于简单；误差在区间左侧对应模型过度拟合数据，可能是模型过于复杂。

3.4 单个误差的兼容性检验

数据误差 $E_{est}$ 和先验信息误差 $E_{est}^p$ 也服从卡方分布，但自由度比总误差少。使用 Welch–Satterthwaite 近似，$E_{est}$ 分配的自由度为 $\nu = \frac{\nu_T N}{N + K}$，$E_{est}^p$ 分配的自由度为 $\nu_p = \frac{\nu_T K}{N + K}$。95% 置信区间近似为：
$\nu - 2\sqrt{2\nu} < E_{est} < \nu + 2\sqrt{2\nu}$
$\nu_p - 2\sqrt{2\nu_p} < E_{est}^p < \nu_p + 2\sqrt{2\nu_p}$

4. 拟合改进的检验

4.1 问题背景

当为单个数据集提出两个替代模型时，虽然一个模型的总误差 $E$ 比另一个小，但由于误差是随机变量，即使从均值相同的概率密度函数中抽取，不同实现的误差大小也会有所不同。因此，需要检验原假设：两个拟合误差的差异仅由随机变化引起。

4.2 F 检验

在没有先验信息且 $E_{est}^T = E_{est}$ 的最简单情况下，使用 F 检验评估两个拟合估计方差比的显著性：
$F_{K_1,K_2} = \frac{\sigma_{d_1}^2_{est}}{\sigma_{d_2}^2_{est}} = \frac{E_{est}^1 / K_1}{E_{est}^2 / K_2}$
其中 $K_1 = N_1 - M_1$，$K_2 = N_2 - M_2$。第一个模型有 $M_1$ 个模型参数、$N_1$ 个数据和总误差 $E_{est}^1$，第二个模型有 $M_2$ 个模型参数、$N_2$ 个数据和总误差 $E_{est}^2$。当计算 $E$ 时使用的单个误差的先验方差都相等时，它们会在 $E_{est}^1 / E_{est}^2$ 中抵消，此时统计量 $F$ 不依赖于数据的先验方差，可以使用未加权误差 $E = \sum_{i}(d_{obs}^i - d_{pre}^i)^2$。

4.3 示例

以 Black Rock Forest 数据集的 51 小时温度上升片段为例，直线拟合（$N_1 = 2$）数据较好，立方函数拟合（$N_2 = 4$）更好，$F_{est} = 1.112$。原假设是两个函数拟合数据效果相同，差异由随机变化引起。双侧检验得到 $P(F \leq 1/F_{est} \text{ 或 } F \geq F_{est}) = 0.71$，远大于 0.05，因此不能拒绝原假设。

5. 频谱峰值显著性检验

5.1 问题提出

复杂时间序列的频谱通常有许多峰谷，我们想知道高振幅峰是否显著。但确定“显著”的含义需要仔细思考。

5.2 常见情况和原假设

常见情况是长时间序列由“噪声”主导，没有明显的正弦模式。原假设是频谱峰值可以由仅包含随机噪声的时间序列中的随机变化解释。

5.3 分析前提

在分析前，需要明确傅里叶系数是在 $(0, f_{ny})$ 还是 $(-f_{ny}, f_{ny})$ 频率范围内定义，这里选择 $(0, f_{ny})$。

5.4 时间序列的处理和功率谱密度的分布

假设长度为 $N$ 的时间序列 $d_i$ 由零均值、方差为 $\sigma_d^2$ 的不相关随机噪声组成。在计算傅里叶变换前，乘以窗函数 $w_i$，得到 $w_id_i$。其方差近似为 $f_f \sigma_d^2$，其中 $f_f = \frac{1}{N} \sum_{i=1}^{N} w_i^2$ 是窗函数的方差。

当数据不相关且正态分布时，傅里叶系数也不相关且正态分布。功率谱密度 $s^2(f)$ 的每个元素与两个傅里叶系数的平方和成正比。将傅里叶系数归一化为单位方差后，功率谱密度服从自由度为 $p = 2$ 的卡方分布。归一化因子 $c$ 为：
$c = \frac{f_f \sigma_d^2}{2N_f \Delta f}$
其中 $N_f = \frac{N}{2} + 1$。功率谱密度的均值为 $s^2 = 2c$，方差为 $\sigma_{s^2}^2 = 4c^2$。

5.5 示例

考虑一个长度为 $N = 1024$ 的时间序列，由 5 Hz 余弦波和不相关噪声组成，余弦波振幅为 $0.25\sigma_d$。功率谱密度在 5 Hz 处有一个突出的峰值，高度约为平均水平的 10 倍。计算得到 $P(s_0^2/c) = 0.99994$，即功率谱密度预计 99.994% 的时间小于该水平。

如果专门寻找 5 Hz 振荡，原假设是 5 Hz 处的峰值仅由随机变化引起，概率为 $1 - P(s_0^2/c) = 0.00006$，可以以很高的置信度拒绝原假设。但如果没有特定频率的先验知识，功率谱密度有 $N/2 + 1 = 513$ 个元素，峰值在这些可能性中出现的概率为 $(0.99994)^{513} = 0.97$，仍有 3% 的可能性由随机变化引起，虽然可以以大于 95% 的置信度拒绝原假设，但置信度远低于前者。

5.6 频谱峰值显著性检验流程

graph TD;
    A[获取时间序列数据] --> B[乘以窗函数];
    B --> C[计算傅里叶变换];
    C --> D[计算功率谱密度];
    D --> E[确定峰值高度和频率];
    E --> F[计算归一化因子 c];
    F --> G[计算 P(speak/c)];
    G --> H{是否有特定频率先验知识};
    H -- 是 --> I[计算特定频率峰值随机变化概率];
    H -- 否 --> J[计算任意频率峰值随机变化概率];
    I --> K{是否拒绝原假设};
    J --> K;
    K -- 是 --> L[峰值显著];
    K -- 否 --> M[峰值不显著];

6. 相关时间序列频谱峰值的显著性评估

6.1 问题背景

许多观测时间序列的相邻样本强相关，与不相关随机时间序列比较不现实。这种时间序列的功率谱密度不均匀，呈红色，大部分功率在低频。

6.2 相关时间序列的构建

用于比较的随机时间序列 $y_i$ 由不相关时间序列 $x_i$ 构建：
$y_i = \begin{cases}
x_1, & i = 1 \
\beta y_{i-1} + x_i, & i = 2, 3, 4, \cdots
\end{cases}$
其中 $x_i$ 是零均值、方差为 $\sigma^2$ 的不相关随机数，$\beta$ 是 0 到 1 之间的常数。$y_i$ 是自回归 AR1 过程的结果，比 $x_i$ 更平滑。

6.3 功率谱密度的计算

通过 z 变换技术计算 $y_i$ 的功率谱密度。将 $y_i - \beta y_{i-1} = x_i$ 转化为 IIR 滤波器形式，得到 $v^ y = u^ x$，其中 $v = \begin{bmatrix} 1 \ -\beta \end{bmatrix}$，$u = [1]$。

$v(z) = 1 - \beta z$，其傅里叶变换为 $v(f) = 1 - \beta \exp(-\frac{\pi if}{f_{ny}})$。$y(f)$ 的频谱为：
$|y(f)|^2 = \frac{|x(f)|^2}{|v(f)|^2}$
其中 $|v(f)|^2 = 1 + \beta^2 - 2\beta \cos(\frac{\pi f}{f_{ny}})$。

从频谱统计的角度，$\frac{1}{|v(f)|^2}$ 是一个乘法常数。$y(f)$ 的功率谱密度的均值和方差与 $x(f)$ 的关系为：
$\text{mean}(|y(f)|^2) = \frac{\text{mean}(|x(f)|^2)}{|v(f)|^2}$
$\text{var}(|y(f)|^2) = \frac{\text{var}(|x(f)|^2)}{|v(f)|^4}$

已知 $\text{mean}(|x(f)|^2) = 2c$，$\text{var}(|x(f)|^2) = 4c^2$，则 $|y(f)|^2|v(f)|^2/c$ 服从自由度为 2 的卡方分布。$y(t)$ 的功率谱密度的均值为 $s^2 = \frac{2c}{|v(f)|^2}$，方差为 $\sigma_{s^2}^2 = \frac{4c^2}{|v(f)|^4}$。

6.4 相关时间序列频谱峰值显著性评估流程

graph TD;
    A[获取不相关时间序列 xi] --> B[构建相关时间序列 yi];
    B --> C[计算 v(f) 和 |v(f)|^2];
    C --> D[计算 xi 的功率谱密度 |x(f)|^2];
    D --> E[计算 yi 的功率谱密度 |y(f)|^2];
    E --> F[确定 yi 频谱峰值高度和频率];
    F --> G[计算 |y(f)|^2|v(f)|^2/c];
    G --> H[计算峰值随机变化概率];
    H --> I{是否拒绝原假设};
    I -- 是 --> J[峰值显著];
    I -- 否 --> K[峰值不显著];

7. 总结

本文介绍了环境数据分析中多种统计检验方法，包括粒子大小数据的统计检验、广义最小二乘法的卡方检验、拟合改进的检验、频谱峰值显著性检验以及相关时间序列频谱峰值的显著性评估。这些方法在评估数据质量、模型拟合效果和判断数据特征显著性方面具有重要作用。在实际应用中，需要根据具体问题选择合适的检验方法，并仔细考虑原假设的设定和结果的解释。

8. 各检验方法的对比与应用场景

8.1 检验方法对比

检验方法	适用场景	原假设	检验统计量	关键参数
粒子大小数据的统计检验	评估仪器校准和方差是否符合要求	差异由随机变化引起	F 统计量	NA, NB, F
广义最小二乘法的卡方检验	评估广义误差与期望值差异是否由随机变化引起	差异由随机变化引起	卡方统计量	$\nu_T$, $\nu$, $\nu_p$
拟合改进的检验	比较两个替代模型的拟合效果	两个拟合误差的差异仅由随机变化引起	F 统计量	$K_1$, $K_2$, $F_{K_1,K_2}$
频谱峰值显著性检验	判断频谱中高振幅峰值是否显著	频谱峰值可以由仅包含随机噪声的时间序列中的随机变化解释	卡方统计量	$c$, $p$, $s^2$
相关时间序列频谱峰值的显著性评估	评估相关时间序列中频谱峰值的显著性	频谱峰值可以由相关随机时间序列中的随机变化解释	卡方统计量	$c$, $

8.2 应用场景选择

• 粒子大小数据的统计检验 ：适用于仪器校准和方差的质量控制，确保仪器测量结果的准确性和稳定性。
• 广义最小二乘法的卡方检验 ：在使用广义最小二乘法进行模型拟合时，用于评估模型误差是否合理，判断模型复杂度是否合适。
• 拟合改进的检验 ：当有多个模型可供选择时，用于确定哪个模型的拟合效果更好，避免因随机误差而误判。
• 频谱峰值显著性检验 ：在分析时间序列频谱时，用于识别有意义的频率成分，排除随机噪声的干扰。
• 相关时间序列频谱峰值的显著性评估 ：对于相邻样本强相关的时间序列，通过与相关随机时间序列比较，准确评估频谱峰值的显著性。

9. 代码实现的注意事项

9.1 MATLAB 代码

在使用 MATLAB 进行统计检验时，需要注意函数的输入参数和返回值。例如，在粒子大小数据的统计检验中， fcdf 函数用于计算 F 分布的累积分布函数，需要正确输入 F 值、自由度等参数。

% edama_12_07, statistical tests of particle size data
...
if( F<1 )
    F=1/F;
end
PF = 1 - (fcdf(F,NA,NB)-fcdf(1/F,NA,NB));

9.2 Python 代码

Python 中使用 scipy.stats 库进行统计检验。在频谱峰值显著性检验中， stats.chi2.cdf 函数用于计算卡方分布的累积分布函数。

# edapy_12_11, confidence of spectral peak
...
ppeak = stats.chi2.cdf(speak/c,p);

在编写代码时，还需要注意数据的预处理，如时间序列的窗函数处理、数据的归一化等，以确保统计检验的准确性。

10. 实际案例分析

10.1 环境温度数据拟合

以 Black Rock Forest 数据集的温度数据为例，使用直线和立方函数进行拟合。通过拟合改进的检验，判断立方函数是否真的比直线拟合效果更好。

模型	模型参数数量	数据数量	总误差	F 统计量	P 值
直线拟合	2	51	$E_{est}^1$	1.112	0.71
立方函数拟合	4	51	$E_{est}^2$	-	-

由于 P 值远大于 0.05，不能拒绝原假设，说明两个函数拟合效果的差异可能是由随机变化引起的。

10.2 时间序列频谱分析

考虑一个长度为 1024 的时间序列，由 5 Hz 余弦波和不相关噪声组成。通过频谱峰值显著性检验，判断 5 Hz 处的峰值是否显著。

• 专门寻找 5 Hz 振荡时，概率为 0.00006，可以以很高的置信度拒绝原假设，认为 5 Hz 处的峰值是显著的。
• 没有特定频率的先验知识时，峰值在所有可能性中出现的概率为 0.97，仍有 3% 的可能性由随机变化引起，虽然可以以大于 95% 的置信度拒绝原假设，但置信度远低于前者。

11. 结论

11.1 统计检验的重要性

环境数据分析中的统计检验方法为评估数据质量、模型拟合效果和判断数据特征显著性提供了有力工具。通过合理运用这些方法，可以避免因随机误差而做出错误的判断，提高分析结果的可靠性。

11.2 方法选择与应用

在实际应用中，需要根据具体问题选择合适的检验方法，并仔细考虑原假设的设定和结果的解释。同时，要注意代码实现的细节，确保数据预处理和统计计算的准确性。

11.3 未来展望

随着环境数据的不断丰富和分析需求的不断提高，统计检验方法也将不断发展和完善。未来可以进一步研究如何结合多种检验方法，提高分析的准确性和可靠性，以及如何将这些方法应用于更复杂的环境问题中。

11.4 实际应用建议

• 在进行统计检验前，充分了解数据的特点和问题的背景，合理设定原假设。
• 对数据进行适当的预处理，如去除异常值、归一化等，以提高检验的准确性。
• 结合多种检验方法进行分析，相互验证结果，减少误判的可能性。
• 在解释检验结果时，要考虑实际问题的意义，避免过度依赖统计数据。

总之，环境数据分析中的统计检验方法是一个复杂而重要的领域，需要不断学习和实践，才能更好地应用这些方法解决实际问题。