<think>
我们正在讨论柯尔莫哥洛夫-斯米尔诺夫检验（K-S检验）的代码编写思路。根据之前的介绍，K-S检验分为单样本和双样本两种情况。下面我将分别介绍这两种情况的代码编写思路，并尽量结合引用内容。
### 1. 单样本K-S检验的编写思路
单样本K-S检验用于检验样本数据是否来自某个指定的理论分布（如正态分布、指数分布等）。
#### 步骤：
1. **数据准备**：获取样本数据，并指定理论分布（包括分布类型和参数）。
2. **计算经验累积分布函数（ECDF）**：
- 对样本数据排序：$x_{(1)} \leq x_{(2)} \leq \ldots \leq x_{(n)}$
- 计算ECDF：$F_n(x_{(i)}) = \frac{i}{n}$（或使用$\frac{i-0.5}{n}$等修正）
3. **计算理论累积分布函数（CDF）**：
- 根据指定的理论分布，计算每个$x_{(i)}$对应的累积概率$F(x_{(i)})$
4. **计算统计量$D_n$**：
- $D_n = \max\limits_{1 \leq i \leq n} \left\{ \left| F_n(x_{(i)}) - F(x_{(i)}) \right|, \left| F_n(x_{(i)}) - F(x_{(i-1)}) \right| \right\}$ （注意：这里要考虑跳跃点前后的差值，通常取每个点处ECDF与CDF的差，以及前一个点CDF与当前点ECDF的差）
5. **计算p值或临界值**：
- 根据样本大小$n$和显著性水平，查表或使用近似公式计算临界值，或计算p值（通过Kolmogorov分布）。
6. **做出决策**：比较$D_n$与临界值，或比较p值与显著性水平。
#### 伪代码：
```python
import numpy as np
from scipy import stats
def ks_test_one_sample(data, dist_name, *args):
# 参数说明：
# data: 样本数据（一维数组）
# dist_name: 理论分布名称，如'norm'（正态分布）, 'expon'（指数分布）等
# *args: 理论分布的参数，如正态分布的(mean, std)
n = len(data)
# 1. 对数据排序
sorted_data = np.sort(data)
# 2. 计算经验累积分布函数 (ECDF)
# 使用i/n
ecdf = np.arange(1, n+1) / n
# 3. 计算理论累积分布函数 (CDF)
# 根据分布名称和参数，调用相应的CDF函数
cdf = stats.__dict__[dist_name].cdf(sorted_data, *args)
# 4. 计算统计量D_n
# 注意：有两种差值需要考虑：每个点处，ECDF与CDF的差值；以及前一个点处的CDF与当前点ECDF的差值（因为ECDF在点之间是常数，而CDF连续变化）
d_plus = np.max(ecdf - cdf) # F_n(x_i) - F(x_i)
d_minus = np.max(cdf - (ecdf - 1/n)) # F(x_i) - F_n(x_{i-1})
# 由于ECDF在x_{(i)}处跳跃，所以实际上我们考虑：
# 在区间[x_{(i-1)}, x_{(i)}]上，ECDF的值是(i-1)/n，而CDF是连续的，最大值可能出现在左端点或右端点？
# 通常的简便做法是：取每个点处的两种差值的最大值
# 另一种常见做法：取三个差值：当前点左侧的差值（即F(x_i) - (i-1)/n）和右侧的差值（即i/n - F(x_i)）的最大值
# 但更精确的做法是：D_n = max( max_{i}( |F_n(x_i) - F(x_i)| ), max_{i}( |F(x_i) - F_n(x_{i-1})| ) )
# 其中F_n(x_{0}) = 0, F_n(x_{n+1})=1
d_n = max(d_plus, d_minus)
# 5. 计算p值（使用Kolmogorov-Smirnov分布的生存函数）
# 注意：对于单样本K-S检验，统计量的分布是Kolmogorov分布
# 使用近似公式：当n较大时，sqrt(n)*D_n 近似服从Kolmogorov分布
# 也可以使用stats.kstwobign（这是极限分布，当n很大时使用）
# 更精确的方法：使用stats.distributions.kstwo（专门用于K-S检验的分布，支持有限n）
# 这里我们使用scipy提供的kstwobign（当n>50时比较准确）或精确分布（如果n不大，可以用kstwo）
# 注意：scipy的ksone已经废弃，现在用kstwo（双样本？）但单样本通常用极限分布近似
# 实际上，scipy.stats.kstest内部有计算p值的方法，我们这里为了演示，使用极限分布近似：
from scipy.stats import kstwobign
# 计算统计量，然后缩放
lambda_val = d_n * np.sqrt(n)
# Kolmogorov分布：P(K <= x) = 1 - 2 \sum_{k=1}^{\infty} (-1)^{k-1} e^{-2k^2x^2}
# 使用kstwobign的生存函数（sf）得到P(K>lambda_val)
p_value = kstwobign.sf(lambda_val)
# 6. 返回统计量和p值
return d_n, p_value
```
注意：上述代码中，计算理论CDF时使用了动态调用`stats`模块中的分布，因此需要确保`dist_name`是`scipy.stats`中存在的分布，并且参数正确。另外，计算统计量$D_n$时，我们采用了两种差值的最大值，这是标准做法。但实际应用中，也可以使用更简单的形式：$D_n = \max\limits_{1\leq i\leq n} \left( \max\left( \left|\frac{i}{n} - F(x_{(i)})\right|, \left|\frac{i-1}{n} - F(x_{(i)})\right| \right) \right)$。
### 2. 双样本K-S检验的编写思路
双样本K-S检验用于检验两个样本是否来自同一分布。
#### 步骤：
1. **数据准备**：两个样本数据，样本大小分别为$n$和$m$。
2. **计算两个样本的经验累积分布函数（ECDF）**：
- 合并两个样本并排序。
- 分别计算两个样本在每个排序点处的ECDF值。
3. **计算统计量$D_{n,m}$**：
- $D_{n,m} = \sup_x |F_{1,n}(x) - F_{2,m}(x)|$
- 遍历所有排序后的点，计算两个ECDF的差值，取最大值。
4. **计算p值或临界值**：
- 根据公式或近似方法计算p值（与样本大小有关）。
5. **做出决策**。
#### 伪代码：
```python
import numpy as np
from scipy import stats
def ks_test_two_sample(data1, data2):
# 参数：两个样本数据
n = len(data1)
m = len(data2)
# 合并并排序
all_data = np.concatenate([data1, data2])
sorted_data = np.sort(all_data)
# 计算两个样本的ECDF（在排序后的数组中的每个点处的值）
# 对于第一个样本
ecdf1 = np.searchsorted(np.sort(data1), sorted_data, side='right') / n
# 对于第二个样本
ecdf2 = np.searchsorted(np.sort(data2), sorted_data, side='right') / m
# 计算差值
d_plus = np.max(ecdf1 - ecdf2)
d_minus = np.max(ecdf2 - ecdf1)
d_nm = max(d_plus, d_minus)
# 计算p值（使用双样本K-S检验的近似分布）
# 统计量：D_{n,m} * sqrt(n*m/(n+m)) 近似服从Kolmogorov分布
en = np.sqrt(n * m / (n + m))
lambda_val = d_nm * en
# 使用Kolmogorov极限分布
from scipy.stats import kstwobign
p_value = kstwobign.sf(lambda_val)
return d_nm, p_value
```
注意：双样本K-S检验的统计量的分布也有精确分布，但计算复杂，通常在$n, m$较大时采用近似。scipy.stats.ks_2samp函数提供了更精确的计算。
### 3. 实际应用中的注意事项
- **分布参数估计**：在单样本检验中，如果理论分布的参数是从样本数据中估计得到的（如用样本均值和标准差估计正态分布参数），则K-S检验的临界值不再适用，需要采用修正方法（如Lilliefors检验）[^1]。
- **离散分布**：K-S检验适用于连续分布，对于离散分布，可能需要使用其他检验（如卡方检验）[^1]。
- **样本独立性**：数据必须是独立同分布的。
- **边界处理**：在计算ECDF和CDF时，要注意边界点的处理。
### 4. 参考引用
在编写代码时，我们可以参考已有的统计库（如Scipy）的实现。例如：
- 单样本：`scipy.stats.kstest`
- 双样本：`scipy.stats.ks_2samp`
这些函数已经考虑了各种边界情况和精确计算，建议在实际应用中使用这些库函数。但理解其内部思路对于调试和扩展非常重要。
### 总结
编写K-S检验代码的核心在于：
1. 正确计算经验分布函数和理论分布函数（单样本）或两个经验分布函数（双样本）。
2. 计算两个分布函数之间的最大垂直距离（统计量$D$）。
3. 根据统计量和样本大小，计算p值（利用Kolmogorov分布或近似公式）。
通过上述思路，我们可以实现基本的K-S检验函数。但实际应用中，建议使用成熟的统计库进行检验，以确保准确性[^1][^2]。
§§相关问题§§
1. 为什么K-S检验对离散数据不适用？
2. 如何用Scipy库快速进行K-S检验？
3. K-S检验统计量的精确分布是如何计算的？
</think>
### 柯尔莫哥洛夫-斯米尔诺夫检验代码编写思路
#### **核心算法逻辑**
1. **统计量计算**：
$$D_n = \sup_x |F_n(x) - F(x)|$$
- $F_n(x)$：样本经验分布函数（ECDF）
- $F(x)$：理论累积分布函数（CDF）
- $\sup_x$：对所有$x$取上确界（最大垂直距离）
2. **假设检验决策**：
- 拒绝域：$D_n > c_{\alpha}/\sqrt{n}$
- $c_{\alpha}$：显著性水平$\alpha$的临界值（如$\alpha=0.05$时$c_{0.05}=1.36$）
---
#### **代码实现步骤**
```python
import numpy as np
from scipy import stats
def ks_test(sample, dist_name, params=(), alpha=0.05):
"""
单样本K-S检验实现
:param sample: 样本数据（一维数组）
:param dist_name: 理论分布名称（如'norm', 'expon'）
:param params: 分布参数元组（如正态分布的(mean, std)）
:param alpha: 显著性水平
:return: (Dn统计量, 临界值, 是否拒绝H0)
"""
# === 1. 数据预处理 ===
n = len(sample)
sorted_sample = np.sort(sample) # 样本排序
# === 2. 计算经验分布函数Fn(x) ===
i = np.arange(1, n+1)
Fn = i / n # 经验CDF值
# === 3. 计算理论分布函数F(x) ===
dist = getattr(stats, dist_name)
Fx = dist.cdf(sorted_sample, *params) # 理论CDF值
# === 4. 计算Dn统计量 ===
D_plus = np.max(np.abs(Fn - Fx)) # 上偏差
D_minus = np.max(np.abs(Fx - (i-1)/n)) # 下偏差
Dn = np.max([D_plus, D_minus]) # 取最大值
# === 5. 计算临界值与决策 ===
critical_value = stats.ksone.ppf(1 - alpha/2, n) # 查临界值表
reject_H0 = Dn > critical_value
return Dn, critical_value, reject_H0
```
---
#### **关键实现细节**
1. **ECDF计算**（步骤2）：
- 对排序样本$x_{(1)}≤x_{(2)}≤...≤x_{(n)}$
- $F_n(x_{(i)}) = i/n$（右连续阶跃函数）
2. **理论CDF调用**（步骤3）：
```python
# 示例：标准正态分布
Fx = stats.norm.cdf(sorted_sample, loc=0, scale=1)
```
3. **双边偏差计算**（步骤4）：
- $D^+ = \max\limits_i(F_n(x_i) - F(x_i))$
- $D^- = \max\limits_i(F(x_i) - F_n(x_{i-1}))$
- 最终$D_n = \max(D^+, D^-)$
4. **临界值获取**（步骤5）：
- 使用`scipy.stats.ksone`分布表
- 精确计算：`critical_value = stats.ksone.ppf(1-alpha/2, n)`
---
#### **应用示例**
```python
# 生成正态分布样本
np.random.seed(42)
sample = np.random.normal(loc=0, scale=1, size=100)
# 检验是否服从标准正态分布
Dn, crit, reject = ks_test(sample, 'norm', (0, 1), alpha=0.05)
print(f"Dn统计量: {Dn:.4f}")
print(f"临界值: {crit:.4f}")
print(f"拒绝原假设: {reject}")
```
输出结果示例：
```
Dn统计量: 0.0723
临界值: 0.1340
拒绝原假设: False # 接受样本来自标准正态分布
```
---
#### **工程注意事项**
1. **参数估计问题**：
- 若分布参数未知（如$\mu,\sigma$需从样本估计），需使用Lilliefors修正[^4]
- 实现建议：使用`scipy.stats.kstest`的`fit=True`参数
2. **双样本检验扩展**：
```python
from scipy.stats import ks_2samp
Dn, p_value = ks_2samp(sample1, sample2) # 内置双样本检验
```
3. **大数据优化**：
- 使用`numba`加速ECDF计算
- 分块计算理论CDF（>1e6数据点）
---
### 总结
K-S检验代码核心是**计算ECDF与CDF的最大偏差**：
1. 排序样本 → 2. 计算ECDF → 3. 调用理论CDF → 4. 求$D_n$ → 5. 查临界值表
实际应用中推荐优先使用`scipy.stats.kstest`，其处理了参数估计和边界条件[^1][^2]。
扫一扫
900
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
