未记录事实核验 · 状态：维护中
记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。
warning 适用边界与维护说明
本文适合学习卡方分布及其概率计算、模拟、拟合优度检验、独立性检验和总体方差区间估计，示例数据不构成真实药物疗效或人群研究结论。代码宜按正文顺序运行，以复用前文导入的库。卡方可加性要求变量独立；频数检验的卡方近似依赖独立观测和足够的期望频数，拟合优度自由度还须考虑样本估计的参数；总体方差的卡方区间限于独立同分布正态样本。效应量分级仅供示例解释，Pearson残差和样本方差比不能单独判定显著性，不拒绝原假设不等于证明独立或方差相等。正文将班级A的Q-Q图称为“正态性检验”，应理解为图形诊断而非正式检验；可加性示例末尾的结论为固定输出，未随KS检验结果判断，模拟也不能证明分布相同或代替理论证明，两组独立样本的逐项差值不代表生成方法的误差。
🟧 卡方分布（Chi-square Distribution）
🤔 什么是卡方分布？
note 📖 **卡方分布的定义**
卡方分布是由多个独立标准正态分布的平方和构成的连续概率分布。如果X₁, X₂, ..., Xₖ 是独立的标准正态随机变量，那么它们的平方和 χ² = X₁² + X₂² + ... + Xₖ² 服从自由度为k的卡方分布。
数学表达式：
- 概率密度函数：$$f(x) = \frac{1}{2^{k/2}\Gamma(k/2)} x^{k/2-1} e^{-x/2}, \quad x > 0$$
- 记号：$\chi^2 \sim \chi^2(k)$，其中 $k$ 是自由度
重要特性：
- 只取非负值（$x \geq 0$）
- 右偏分布，随自由度增大趋于对称
- 均值 = 自由度 $k$
- 方差 = $2k$
🏗️ 生活中的卡方分布
想象你在做一个"投掷色子"的实验：
example 🎲 **色子公平性检验**
你怀疑一个六面色子是否公平。投掷600次，每个面理论上应该出现100次。实际结果：
- 面1：95次，面2：110次，面3：85次
- 面4：120次，面5：90次，面6：100次
卡方统计量 = Σ(观察值-期望值)²/期望值；在公平色子的原假设成立、各次投掷独立且期望频数足够大时，该统计量的分布近似为自由度5的卡方分布。
其他应用场景：
- 检验药物疗效是否与性别独立
- 分析不同地区销售数据是否存在显著差异
- 验证抽样调查结果是否符合总体分布
- 质量控制中的方差检验
📊 卡方分布的性质和特征
1. 不同自由度的分布形状
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
def plot_chi2_properties():
"""展示卡方分布的基本性质"""
# 设置中文字体
plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. 不同自由度的PDF
x = np.linspace(0, 25, 1000)
degrees_of_freedom = [1, 2, 3, 5, 10, 15]
colors = ['red', 'blue', 'green', 'orange', 'purple', 'brown']
for df, color in zip(degrees_of_freedom, colors):
pdf = stats.chi2.pdf(x, df)
ax1.plot(x, pdf, color=color, linewidth=2, label=f'df = {df}')
ax1.set_xlabel('x')
ax1.set_ylabel('概率密度')
ax1.set_title('不同自由度的卡方分布')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 25)
# 2. 均值和方差的变化
df_range = np.arange(1, 21)
means = df_range # 均值 = df
variances = 2 * df_range # 方差 = 2*df
std_devs = np.sqrt(variances)
ax2.plot(df_range, means, 'b-o', label='均值 = df', linewidth=2, markersize=6)
ax2.plot(df_range, std_devs, 'r-s', label='标准差 = √(2df)', linewidth=2, markersize=6)
ax2.set_xlabel('自由度 (df)')
ax2.set_ylabel('数值')
ax2.set_title('卡方分布的均值与标准差')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 3. 累积分布函数
x_cdf = np.linspace(0, 20, 1000)
for df, color in zip([1, 3, 5, 10], colors[:4]):
cdf = stats.chi2.cdf(x_cdf, df)
ax3.plot(x_cdf, cdf, color=color, linewidth=2, label=f'df = {df}')
ax3.set_xlabel('x')
ax3.set_ylabel('累积概率')
ax3.set_title('卡方分布的累积分布函数')
ax3.legend()
ax3.grid(True, alpha=0.3)
ax3.set_xlim(0, 20)
# 4. 分位数图
percentiles = np.array([0.9, 0.95, 0.975, 0.99, 0.999])
df_values = [1, 2, 3, 5, 10]
for df in df_values:
quantiles = stats.chi2.ppf(percentiles, df)
ax4.plot(percentiles * 100, quantiles, 'o-', label=f'df = {df}', linewidth=2, markersize=6)
ax4.set_xlabel('百分位数 (%)')
ax4.set_ylabel('分位数值')
ax4.set_title('卡方分布的分位数')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 输出重要性质
print("卡方分布的重要性质:")
print("=" * 40)
print("1. 定义域: x ≥ 0 (只取非负值)")
print("2. 右偏分布，df增大时趋于对称")
print("3. 均值 E(X) = df")
print("4. 方差 Var(X) = 2·df")
print("5. 标准差 SD(X) = √(2·df)")
print("6. 模式 Mode = max(0, df-2)")
print("7. 可加性: 若X₁ ~ χ²(df₁)、X₂ ~ χ²(df₂)且相互独立，则X₁ + X₂ ~ χ²(df₁+df₂)")
plot_chi2_properties()
2. 数值特征验证
def verify_chi2_properties():
"""验证卡方分布的数值特征"""
df_values = [1, 3, 5, 10, 20]
n_samples = 100000
print("卡方分布数值特征验证:")
print("=" * 50)
print("自由度\t理论均值\t模拟均值\t理论方差\t模拟方差\t理论标准差\t模拟标准差")
print("-" * 80)
np.random.seed(42)
for df in df_values:
# 生成随机样本
samples = stats.chi2.rvs(df, size=n_samples)
# 理论值
theory_mean = df
theory_var = 2 * df
theory_std = np.sqrt(2 * df)
# 模拟值
sim_mean = np.mean(samples)
sim_var = np.var(samples, ddof=1)
sim_std = np.std(samples, ddof=1)
print(f"{df}\t{theory_mean:.2f}\t\t{sim_mean:.2f}\t\t{theory_var:.2f}\t\t{sim_var:.2f}\t\t{theory_std:.2f}\t\t{sim_std:.2f}")
# 可加性验证
print(f"\n可加性验证:")
print("-" * 30)
df1, df2 = 3, 5
df_sum = df1 + df2
# 独立生成两个卡方分布
chi2_1 = stats.chi2.rvs(df1, size=n_samples)
chi2_2 = stats.chi2.rvs(df2, size=n_samples)
chi2_sum = chi2_1 + chi2_2
# 直接生成df1+df2的卡方分布
chi2_direct = stats.chi2.rvs(df_sum, size=n_samples)
print(f"χ²({df1}) + χ²({df2}) 的均值: {np.mean(chi2_sum):.3f}")
print(f"χ²({df_sum}) 的均值: {np.mean(chi2_direct):.3f}")
print(f"理论均值: {df_sum}")
print(f"\nχ²({df1}) + χ²({df2}) 的方差: {np.var(chi2_sum, ddof=1):.3f}")
print(f"χ²({df_sum}) 的方差: {np.var(chi2_direct, ddof=1):.3f}")
print(f"理论方差: {2 * df_sum}")
verify_chi2_properties()
💻 Python实现详解
1. 基本计算
from scipy import stats
import numpy as np
def chi2_calculations():
"""卡方分布的各种计算"""
df = 5 # 自由度为5
print(f"卡方分布 χ²({df}) 的计算:")
print("=" * 30)
# 1. 概率密度函数值
x_values = [1, 3, 5, 7, 10, 15]
print("概率密度函数值:")
for x in x_values:
pdf_value = stats.chi2.pdf(x, df)
print(f"f({x}) = {pdf_value:.4f}")
# 2. 累积分布函数值
print(f"\n累积分布函数值:")
for x in [2, 5, 7, 10, 15]:
cdf_value = stats.chi2.cdf(x, df)
print(f"P(X ≤ {x}) = {cdf_value:.4f}")
# 3. 上尾概率（生存函数）
print(f"\n上尾概率:")
for x in [5, 7, 10, 15]:
sf_value = stats.chi2.sf(x, df)
print(f"P(X > {x}) = {sf_value:.4f}")
# 4. 分位数
print(f"\n重要分位数:")
percentiles = [0.9, 0.95, 0.975, 0.99, 0.995, 0.999]
for p in percentiles:
quantile = stats.chi2.ppf(p, df)
print(f"{p:.1%}分位数: {quantile:.4f}")
# 5. 区间概率
print(f"\n区间概率:")
intervals = [(2, 8), (3, 10), (5, 15)]
for a, b in intervals:
prob = stats.chi2.cdf(b, df) - stats.chi2.cdf(a, df)
print(f"P({a} ≤ X ≤ {b}) = {prob:.4f}")
# 6. 统计特征
print(f"\n统计特征:")
mean = stats.chi2.mean(df)
var = stats.chi2.var(df)
std = stats.chi2.std(df)
print(f"均值: {mean:.2f}")
print(f"方差: {var:.2f}")
print(f"标准差: {std:.2f}")
return df
df = chi2_calculations()
2. 随机数生成和可视化
def chi2_simulation_and_visualization():
"""卡方分布模拟和可视化"""
df = 8
n_samples = 10000
print(f"卡方分布 χ²({df}) 模拟分析:")
print("=" * 35)
# 生成随机数
np.random.seed(42)
samples = stats.chi2.rvs(df, size=n_samples)
# 基本统计
print(f"样本量: {n_samples:,}")
print(f"样本均值: {np.mean(samples):.3f} (理论值: {df})")
print(f"样本方差: {np.var(samples, ddof=1):.3f} (理论值: {2*df})")
print(f"样本标准差: {np.std(samples, ddof=1):.3f} (理论值: {np.sqrt(2*df):.3f})")
print(f"样本最小值: {np.min(samples):.3f}")
print(f"样本最大值: {np.max(samples):.3f}")
# 分位数比较
print(f"\n分位数比较:")
print("百分位\t样本分位数\t理论分位数\t误差")
print("-" * 40)
percentiles_check = [0.25, 0.5, 0.75, 0.9, 0.95]
for p in percentiles_check:
sample_quantile = np.percentile(samples, p * 100)
theory_quantile = stats.chi2.ppf(p, df)
error = abs(sample_quantile - theory_quantile)
print(f"{p:.0%}\t{sample_quantile:.3f}\t\t{theory_quantile:.3f}\t\t{error:.3f}")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. 直方图vs理论PDF
ax1.hist(samples, bins=50, density=True, alpha=0.7, color='skyblue', edgecolor='black')
x_theory = np.linspace(0, max(samples), 1000)
pdf_theory = stats.chi2.pdf(x_theory, df)
ax1.plot(x_theory, pdf_theory, 'r-', linewidth=2, label=f'理论PDF χ²({df})')
ax1.set_xlabel('x')
ax1.set_ylabel('概率密度')
ax1.set_title('样本分布 vs 理论分布')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. Q-Q图
from scipy.stats import probplot
probplot(samples, dist=stats.chi2(df), plot=ax2)
ax2.set_title('Q-Q图: 样本 vs 理论χ²分布')
ax2.grid(True, alpha=0.3)
# 3. 经验CDF vs 理论CDF
x_sorted = np.sort(samples)
empirical_cdf = np.arange(1, len(x_sorted) + 1) / len(x_sorted)
theoretical_cdf = stats.chi2.cdf(x_sorted, df)
ax3.plot(x_sorted, empirical_cdf, label='经验CDF', linewidth=1)
ax3.plot(x_sorted, theoretical_cdf, 'r--', label='理论CDF', linewidth=2)
ax3.set_xlabel('x')
ax3.set_ylabel('累积概率')
ax3.set_title('经验CDF vs 理论CDF')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 4. 箱线图分析
quartiles = np.percentile(samples, [25, 50, 75])
iqr = quartiles[2] - quartiles[0]
ax4.boxplot(samples, vert=True, patch_artist=True,
boxprops=dict(facecolor='lightblue', alpha=0.7))
ax4.set_ylabel('值')
ax4.set_title('样本箱线图')
ax4.grid(True, alpha=0.3)
# 添加统计信息
ax4.text(1.1, quartiles[1], f'中位数: {quartiles[1]:.2f}', fontsize=10)
ax4.text(1.1, quartiles[1] - iqr*0.3, f'IQR: {iqr:.2f}', fontsize=10)
plt.tight_layout()
plt.show()
return samples
samples = chi2_simulation_and_visualization()
🎯 卡方分布的实际应用
1. 拟合优度检验
def goodness_of_fit_test():
"""拟合优度检验案例"""
print("案例1: 色子公平性检验")
print("=" * 30)
# 投掷600次色子的观察结果
observed = np.array([95, 110, 85, 120, 90, 100]) # 各面出现次数
expected = np.array([100, 100, 100, 100, 100, 100]) # 期望次数（公平色子）
faces = ['面1', '面2', '面3', '面4', '面5', '面6']
print("观察结果:")
for face, obs, exp in zip(faces, observed, expected):
print(f"{face}: 观察{obs}次, 期望{exp}次, 偏差{obs-exp:+d}")
# 手动计算卡方统计量
chi2_manual = np.sum((observed - expected)**2 / expected)
df = len(observed) - 1 # 本例概率预先指定且未估计额外参数，自由度 = 类别数 - 1
# 使用scipy计算
chi2_stat, p_value = stats.chisquare(observed, expected)
print(f"\n检验结果:")
print(f"卡方统计量 χ² = {chi2_stat:.4f}")
print(f"自由度 df = {df}")
print(f"p值 = {p_value:.4f}")
# 临界值
alpha = 0.05
critical_value = stats.chi2.ppf(1 - alpha, df)
print(f"\n决策 (α = {alpha}):")
print(f"临界值 χ²_{alpha}({df}) = {critical_value:.4f}")
if chi2_stat > critical_value:
print(f"拒绝原假设: 色子不公平 (χ² = {chi2_stat:.4f} > {critical_value:.4f})")
else:
print(f"不拒绝原假设: 没有足够证据表明色子不公平")
if p_value < alpha:
print(f"p值判断: 拒绝原假设 (p = {p_value:.4f} < {alpha})")
else:
print(f"p值判断: 不拒绝原假设 (p = {p_value:.4f} ≥ {alpha})")
# 可视化
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# 观察值vs期望值
x_pos = np.arange(len(faces))
width = 0.35
bars1 = ax1.bar(x_pos - width/2, observed, width, label='观察值', alpha=0.7, color='skyblue')
bars2 = ax1.bar(x_pos + width/2, expected, width, label='期望值', alpha=0.7, color='lightcoral')
ax1.set_xlabel('色子面')
ax1.set_ylabel('频次')
ax1.set_title('观察频次 vs 期望频次')
ax1.set_xticks(x_pos)
ax1.set_xticklabels(faces)
ax1.legend()
ax1.grid(True, alpha=0.3)
# 在柱子上标注数值
for bar1, bar2, obs, exp in zip(bars1, bars2, observed, expected):
ax1.text(bar1.get_x() + bar1.get_width()/2, bar1.get_height() + 2,
str(obs), ha='center', va='bottom', fontsize=10)
ax1.text(bar2.get_x() + bar2.get_width()/2, bar2.get_height() + 2,
str(exp), ha='center', va='bottom', fontsize=10)
# 卡方分布和检验统计量
x_chi2 = np.linspace(0, 20, 1000)
pdf_chi2 = stats.chi2.pdf(x_chi2, df)
ax2.plot(x_chi2, pdf_chi2, 'b-', linewidth=2, label=f'χ²({df})分布')
ax2.axvline(chi2_stat, color='red', linestyle='--', linewidth=2,
label=f'检验统计量 = {chi2_stat:.3f}')
ax2.axvline(critical_value, color='green', linestyle='--', linewidth=2,
label=f'临界值 = {critical_value:.3f}')
# 填充拒绝域
x_reject = x_chi2[x_chi2 >= critical_value]
pdf_reject = stats.chi2.pdf(x_reject, df)
ax2.fill_between(x_reject, 0, pdf_reject, alpha=0.3, color='red', label='拒绝域')
ax2.set_xlabel('χ² 值')
ax2.set_ylabel('概率密度')
ax2.set_title('卡方检验图示')
ax2.legend()
ax2.grid(True, alpha=0.3)
ax2.set_xlim(0, 20)
plt.tight_layout()
plt.show()
goodness_of_fit_test()
2. 独立性检验
def independence_test():
"""独立性检验案例"""
print("\n\n案例2: 药物疗效与性别独立性检验")
print("=" * 40)
# 2x2列联表：药物疗效 vs 性别
# 行：疗效（有效/无效），列：性别（男/女）
observed_table = np.array([
[60, 40], # 有效：男60人，女40人
[30, 70] # 无效：男30人，女70人
])
# 标签
effect_labels = ['有效', '无效']
gender_labels = ['男性', '女性']
print("观察数据 (列联表):")
print("\t\t男性\t女性\t总计")
row_totals = observed_table.sum(axis=1)
col_totals = observed_table.sum(axis=0)
grand_total = observed_table.sum()
for i, effect in enumerate(effect_labels):
print(f"{effect}\t\t{observed_table[i,0]}\t{observed_table[i,1]}\t{row_totals[i]}")
print(f"总计\t\t{col_totals[0]}\t{col_totals[1]}\t{grand_total}")
# 计算期望频数
expected_table = np.outer(row_totals, col_totals) / grand_total
print(f"\n期望频数:")
print("\t\t男性\t女性")
for i, effect in enumerate(effect_labels):
print(f"{effect}\t\t{expected_table[i,0]:.1f}\t{expected_table[i,1]:.1f}")
# 卡方检验
chi2_stat, p_value, dof, expected = stats.chi2_contingency(observed_table, correction=False) # 使用未作连续性校正的Pearson统计量，与下方标准Cramér's V计算一致
print(f"\n独立性检验结果:")
print(f"卡方统计量 χ² = {chi2_stat:.4f}")
print(f"自由度 df = {dof}")
print(f"p值 = {p_value:.4f}")
# 效应量 - Cramér's V
n = grand_total
cramers_v = np.sqrt(chi2_stat / (n * (min(observed_table.shape) - 1)))
print(f"Cramér's V = {cramers_v:.4f}")
# 判断效应大小
if cramers_v < 0.1:
effect_size = "微弱"
elif cramers_v < 0.3:
effect_size = "小"
elif cramers_v < 0.5:
effect_size = "中等"
else:
effect_size = "大"
print(f"效应大小: {effect_size}")
# 决策
alpha = 0.05
critical_value = stats.chi2.ppf(1 - alpha, dof)
print(f"\n决策 (α = {alpha}):")
if p_value < alpha:
print("拒绝原假设: 药物疗效与性别存在关联")
print("结论: 药物疗效可能因性别而异")
else:
print("不拒绝原假设: 没有足够证据表明药物疗效与性别有关，但不能据此证明二者独立")
print("结论: 没有足够证据表明疗效与性别有关")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. 堆叠柱状图
bottom_values = [0, 0]
colors = ['lightblue', 'lightcoral']
for i, effect in enumerate(effect_labels):
ax1.bar(gender_labels, observed_table[i], bottom=bottom_values,
label=effect, alpha=0.7, color=colors[i])
# 更新底部值用于堆叠
if i == 0:
bottom_values = observed_table[i]
ax1.set_ylabel('人数')
ax1.set_title('药物疗效按性别分布')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. 热图显示观察值
import matplotlib.colors as mcolors
im1 = ax2.imshow(observed_table, cmap='Blues', aspect='auto')
ax2.set_xticks(range(len(gender_labels)))
ax2.set_xticklabels(gender_labels)
ax2.set_yticks(range(len(effect_labels)))
ax2.set_yticklabels(effect_labels)
ax2.set_title('观察频数热图')
# 在格子中添加数值
for i in range(len(effect_labels)):
for j in range(len(gender_labels)):
ax2.text(j, i, f'{observed_table[i, j]}',
ha='center', va='center', fontsize=14, fontweight='bold')
# 3. 残差分析
residuals = (observed_table - expected_table) / np.sqrt(expected_table)
im2 = ax3.imshow(residuals, cmap='RdBu_r', aspect='auto', vmin=-3, vmax=3)
ax3.set_xticks(range(len(gender_labels)))
ax3.set_xticklabels(gender_labels)
ax3.set_yticks(range(len(effect_labels)))
ax3.set_yticklabels(effect_labels)
ax3.set_title('Pearson残差（用于探索偏离，不直接判定单元格显著性）')
# 添加颜色条
plt.colorbar(im2, ax=ax3, shrink=0.8)
# 在格子中添加残差值
for i in range(len(effect_labels)):
for j in range(len(gender_labels)):
color = 'white' if abs(residuals[i, j]) > 1.5 else 'black'
ax3.text(j, i, f'{residuals[i, j]:.2f}',
ha='center', va='center', fontsize=12, color=color, fontweight='bold')
# 4. 比例图
# 计算行百分比（每种疗效中的性别比例）
row_props = observed_table / row_totals[:, np.newaxis] * 100
x_pos = np.arange(len(effect_labels))
width = 0.35
bars1 = ax4.bar(x_pos - width/2, row_props[:, 0], width,
label='男性%', alpha=0.7, color='lightblue')
bars2 = ax4.bar(x_pos + width/2, row_props[:, 1], width,
label='女性%', alpha=0.7, color='lightcoral')
ax4.set_xlabel('疗效')
ax4.set_ylabel('百分比 (%)')
ax4.set_title('各疗效组中的性别比例')
ax4.set_xticks(x_pos)
ax4.set_xticklabels(effect_labels)
ax4.legend()
ax4.grid(True, alpha=0.3)
# 在柱子上标注百分比
for bar1, bar2, male_pct, female_pct in zip(bars1, bars2, row_props[:, 0], row_props[:, 1]):
ax4.text(bar1.get_x() + bar1.get_width()/2, bar1.get_height() + 1,
f'{male_pct:.1f}%', ha='center', va='bottom', fontsize=10)
ax4.text(bar2.get_x() + bar2.get_width()/2, bar2.get_height() + 1,
f'{female_pct:.1f}%', ha='center', va='bottom', fontsize=10)
plt.tight_layout()
plt.show()
# 详细解读
print(f"\n详细解读:")
print("1. 观察与期望的对比表明治疗效果可能存在性别差异")
print("2. 男性有效率: {:.1f}%".format(observed_table[0,0]/col_totals[0]*100))
print("3. 女性有效率: {:.1f}%".format(observed_table[0,1]/col_totals[1]*100))
if abs(residuals[0,0]) > 2:
print("4. 男性有效人数的Pearson残差绝对值超过2，提示值得进一步检查的偏离，不单独构成显著性结论")
if abs(residuals[0,1]) > 2:
print("5. 女性有效人数的Pearson残差绝对值超过2，提示值得进一步检查的偏离，不单独构成显著性结论")
independence_test()
3. 方差齐性检验
def variance_homogeneity_test():
"""方差齐性检验案例"""
print("\n\n案例3: 不同班级成绩方差齐性检验")
print("=" * 40)
# 模拟三个班级的考试成绩
np.random.seed(42)
# 假设班级A方差较小，班级B方差适中，班级C方差较大
class_A = np.random.normal(85, 8, 30) # 均值85，标准差8
class_B = np.random.normal(82, 12, 30) # 均值82，标准差12
class_C = np.random.normal(80, 15, 30) # 均值80，标准差15
classes = {'班级A': class_A, '班级B': class_B, '班级C': class_C}
print("各班级成绩统计:")
print("-" * 30)
print("班级\t样本量\t均值\t标准差\t方差")
class_stats = {}
for name, scores in classes.items():
n = len(scores)
mean_score = np.mean(scores)
std_score = np.std(scores, ddof=1)
var_score = np.var(scores, ddof=1)
class_stats[name] = {'n': n, 'mean': mean_score, 'std': std_score, 'var': var_score}
print(f"{name}\t{n}\t{mean_score:.1f}\t{std_score:.2f}\t\t{var_score:.2f}")
# Bartlett检验：用于独立正态总体样本，原假设下统计量渐近服从自由度为组数减一的卡方分布
from scipy.stats import bartlett
bartlett_stat, bartlett_p = bartlett(class_A, class_B, class_C)
print(f"\nBartlett方差齐性检验:")
print(f"检验统计量: {bartlett_stat:.4f}")
print(f"p值: {bartlett_p:.4f}")
# Levene检验：对非正态数据通常更稳健，此处使用中位数中心；其统计量采用F分布参照，不是卡方分布参照
from scipy.stats import levene
levene_stat, levene_p = levene(class_A, class_B, class_C, center='median')
print(f"\nLevene方差齐性检验:")
print(f"检验统计量: {levene_stat:.4f}")
print(f"p值: {levene_p:.4f}")
# 决策
alpha = 0.05
print(f"\n决策 (α = {alpha}):")
if bartlett_p < alpha:
print(f"Bartlett检验: 拒绝方差齐性假设 (p = {bartlett_p:.4f} < {alpha})")
else:
print(f"Bartlett检验: 不拒绝方差齐性假设 (p = {bartlett_p:.4f} ≥ {alpha})")
if levene_p < alpha:
print(f"Levene检验: 拒绝方差齐性假设 (p = {levene_p:.4f} < {alpha})")
else:
print(f"Levene检验: 不拒绝方差齐性假设 (p = {levene_p:.4f} ≥ {alpha})")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. 箱线图比较
data_for_boxplot = [class_A, class_B, class_C]
labels = list(classes.keys())
bp = ax1.boxplot(data_for_boxplot, labels=labels, patch_artist=True)
colors = ['lightblue', 'lightgreen', 'lightcoral']
for patch, color in zip(bp['boxes'], colors):
patch.set_facecolor(color)
patch.set_alpha(0.7)
ax1.set_ylabel('成绩')
ax1.set_title('各班级成绩分布比较')
ax1.grid(True, alpha=0.3)
# 2. 方差比较柱状图
class_names = list(class_stats.keys())
variances = [class_stats[name]['var'] for name in class_names]
bars = ax2.bar(class_names, variances, alpha=0.7, color=colors)
ax2.set_ylabel('方差')
ax2.set_title('各班级成绩方差比较')
ax2.grid(True, alpha=0.3)
# 在柱子上标注数值
for bar, var in zip(bars, variances):
ax2.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 5,
f'{var:.1f}', ha='center', va='bottom', fontsize=12, fontweight='bold')
# 3. 残差图（用于检查正态性假设）
# 计算标准化残差
all_residuals = []
class_labels_expanded = []
for i, (name, scores) in enumerate(classes.items()):
mean_score = class_stats[name]['mean']
std_score = class_stats[name]['std']
residuals = (scores - mean_score) / std_score
all_residuals.extend(residuals)
class_labels_expanded.extend([i] * len(scores))
# 按班级绘制残差
for i, (name, color) in enumerate(zip(class_names, colors)):
class_residuals = [res for j, res in enumerate(all_residuals) if class_labels_expanded[j] == i]
ax3.scatter([i] * len(class_residuals), class_residuals,
alpha=0.6, color=color, s=30, label=name)
ax3.set_xticks(range(len(class_names)))
ax3.set_xticklabels(class_names)
ax3.set_ylabel('标准化残差')
ax3.set_title('标准化残差分布')
ax3.axhline(y=0, color='black', linestyle='-', alpha=0.5)
ax3.axhline(y=2, color='red', linestyle='--', alpha=0.7)
ax3.axhline(y=-2, color='red', linestyle='--', alpha=0.7)
ax3.grid(True, alpha=0.3)
ax3.legend()
# 4. Q-Q图检验正态性（以班级A为例）
from scipy import stats as scipy_stats
scipy_stats.probplot(class_A, dist="norm", plot=ax4)
ax4.set_title('班级A成绩的Q-Q图 (正态性检验)')
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 进一步分析
print(f"\n进一步分析:")
print("1. 方差差异程度:")
max_var = max(variances)
min_var = min(variances)
variance_ratio = max_var / min_var
print(f" 最大方差/最小方差 = {variance_ratio:.2f}")
if variance_ratio > 4:
print(" 样本方差比超过4，差异较大；是否拒绝方差齐性假设需结合检验结果及其适用条件")
elif variance_ratio > 2:
print(" 方差差异中等，需要谨慎解释")
else:
print(" 样本方差比不超过2，差异相对较小；该比值不能单独证明方差齐性")
print(f"\n2. 实际意义:")
print(" - 如果方差不齐，说明不同班级学生成绩的离散程度不同")
print(" - 这可能反映了教学方法、学生基础等因素的影响")
print(" - 在进行方差分析前需要考虑方差齐性假设是否满足")
variance_homogeneity_test()
🧮 卡方分布的理论应用
1. 样本方差的分布
def sample_variance_distribution():
"""演示样本方差的卡方分布性质"""
print("样本方差的卡方分布性质:")
print("=" * 35)
# 从正态总体N(μ, σ²)独立同分布抽样；下述样本方差卡方分布及总体方差置信区间依赖这一前提
true_mean = 10
true_std = 3
true_var = true_std**2
sample_sizes = [5, 10, 20, 30]
n_experiments = 10000
print(f"总体: N({true_mean}, {true_var})")
print(f"实验次数: {n_experiments:,}")
fig, axes = plt.subplots(2, 2, figsize=(15, 12))
axes = axes.flatten()
for idx, n in enumerate(sample_sizes):
print(f"\n样本量 n = {n}:")
# 进行多次抽样实验
np.random.seed(42)
sample_variances = []
for _ in range(n_experiments):
sample = np.random.normal(true_mean, true_std, n)
sample_var = np.var(sample, ddof=1) # 使用无偏估计
sample_variances.append(sample_var)
sample_variances = np.array(sample_variances)
# 理论：(n-1)S²/σ² ~ χ²(n-1)
chi2_statistics = (n - 1) * sample_variances / true_var
# 统计结果
print(f" 样本方差均值: {np.mean(sample_variances):.3f} (理论: {true_var:.3f})")
print(f" (n-1)S²/σ²的均值: {np.mean(chi2_statistics):.3f} (理论: {n-1})")
print(f" (n-1)S²/σ²的方差: {np.var(chi2_statistics, ddof=1):.3f} (理论: {2*(n-1)})")
# 可视化
ax = axes[idx]
# 直方图
ax.hist(chi2_statistics, bins=50, density=True, alpha=0.7,
color='skyblue', edgecolor='black', label='样本统计量')
# 理论卡方分布
x_theory = np.linspace(0, max(chi2_statistics), 1000)
pdf_theory = stats.chi2.pdf(x_theory, n-1)
ax.plot(x_theory, pdf_theory, 'r-', linewidth=2,
label=f'χ²({n-1})理论分布')
ax.set_title(f'n = {n}, (n-1)S²/σ² ~ χ²({n-1})')
ax.set_xlabel('(n-1)S²/σ²')
ax.set_ylabel('概率密度')
ax.legend()
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 置信区间应用
print(f"\n置信区间应用:")
print("-" * 20)
n = 20
confidence_levels = [0.90, 0.95, 0.99]
print(f"当样本量 n = {n} 时，总体方差σ²的置信区间:")
for conf_level in confidence_levels:
alpha = 1 - conf_level
# 卡方分布的临界值
chi2_lower = stats.chi2.ppf(alpha/2, n-1)
chi2_upper = stats.chi2.ppf(1-alpha/2, n-1)
print(f"\n{conf_level:.0%}置信区间:")
print(f" χ²_{alpha/2:.3f}({n-1}) = {chi2_lower:.4f}")
print(f" χ²_{1-alpha/2:.3f}({n-1}) = {chi2_upper:.4f}")
print(f" σ²的置信区间: [(n-1)S²/{chi2_upper:.3f}, (n-1)S²/{chi2_lower:.3f}]")
# 用一个具体的样本方差来计算
sample_var_example = 10.5
ci_lower = (n-1) * sample_var_example / chi2_upper
ci_upper = (n-1) * sample_var_example / chi2_lower
print(f" 若S² = {sample_var_example}, 则σ²的{conf_level:.0%}置信区间: [{ci_lower:.3f}, {ci_upper:.3f}]")
sample_variance_distribution()
2. 卡方分布的可加性
def chi2_additivity():
"""验证卡方分布的可加性"""
print("\n\n卡方分布的可加性验证:")
print("=" * 30)
print("定理: 若 X₁ ~ χ²(df₁), X₂ ~ χ²(df₂) 且独立")
print("则 X₁ + X₂ ~ χ²(df₁ + df₂)")
# 设定参数
df1, df2 = 3, 7
df_sum = df1 + df2
n_samples = 50000
print(f"\n验证: χ²({df1}) + χ²({df2}) = χ²({df_sum})")
print(f"样本量: {n_samples:,}")
np.random.seed(42)
# 方法1: 分别生成两个卡方分布然后相加
chi2_1 = stats.chi2.rvs(df1, size=n_samples)
chi2_2 = stats.chi2.rvs(df2, size=n_samples)
chi2_sum = chi2_1 + chi2_2
# 方法2: 直接生成df1+df2的卡方分布
chi2_direct = stats.chi2.rvs(df_sum, size=n_samples)
# 方法3: 通过正态分布构造
# χ²(k) = Z₁² + Z₂² + ... + Zₖ²，其中 Zᵢ ~ N(0,1)
normals_sum = np.sum(np.random.standard_normal((n_samples, df_sum))**2, axis=1)
# 统计比较
methods = {
'χ²(3) + χ²(7)': chi2_sum,
'χ²(10)直接生成': chi2_direct,
'10个N(0,1)²之和': normals_sum
}
print("\n统计特征比较:")
print("方法\t\t\t均值\t方差\t标准差\tQ1\tQ2\tQ3")
print("-" * 80)
theoretical_mean = df_sum
theoretical_var = 2 * df_sum
theoretical_std = np.sqrt(theoretical_var)
print(f"理论值\t\t\t{theoretical_mean:.2f}\t{theoretical_var:.2f}\t{theoretical_std:.2f}\t-\t-\t-")
for method_name, data in methods.items():
mean_val = np.mean(data)
var_val = np.var(data, ddof=1)
std_val = np.std(data, ddof=1)
q1, q2, q3 = np.percentile(data, [25, 50, 75])
print(f"{method_name}\t{mean_val:.2f}\t{var_val:.2f}\t{std_val:.2f}\t{q1:.2f}\t{q2:.2f}\t{q3:.2f}")
# KS检验比较分布相似性
from scipy.stats import kstest, ks_2samp
print(f"\n分布相似性检验 (KS检验):")
print("-" * 35)
# 与理论分布比较
for method_name, data in methods.items():
ks_stat, p_value = kstest(data, lambda x: stats.chi2.cdf(x, df_sum))
print(f"{method_name} vs χ²({df_sum})理论分布:")
print(f" KS统计量: {ks_stat:.6f}, p值: {p_value:.6f}")
# 两两比较
print(f"\n方法间比较:")
method_pairs = [
('χ²(3) + χ²(7)', 'χ²(10)直接生成'),
('χ²(3) + χ²(7)', '10个N(0,1)²之和'),
('χ²(10)直接生成', '10个N(0,1)²之和')
]
for method1, method2 in method_pairs:
ks_stat, p_value = ks_2samp(methods[method1], methods[method2])
print(f"{method1} vs {method2}:")
print(f" KS统计量: {ks_stat:.6f}, p值: {p_value:.6f}")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. PDF比较
x = np.linspace(0, 30, 1000)
theoretical_pdf = stats.chi2.pdf(x, df_sum)
ax1.plot(x, theoretical_pdf, 'k-', linewidth=3, label=f'理论χ²({df_sum})')
colors = ['red', 'blue', 'green']
alphas = [0.7, 0.7, 0.7]
for i, (method_name, data) in enumerate(methods.items()):
ax1.hist(data, bins=60, density=True, alpha=alphas[i],
color=colors[i], label=method_name, edgecolor='white', linewidth=0.5)
ax1.set_xlabel('x')
ax1.set_ylabel('概率密度')
ax1.set_title('卡方分布可加性验证 - PDF比较')
ax1.legend()
ax1.grid(True, alpha=0.3)
ax1.set_xlim(0, 30)
# 2. CDF比较
x_cdf = np.sort(chi2_sum)
theoretical_cdf = stats.chi2.cdf(x_cdf, df_sum)
empirical_cdf = np.arange(1, len(x_cdf) + 1) / len(x_cdf)
ax2.plot(x_cdf, theoretical_cdf, 'k-', linewidth=2, label='理论CDF')
ax2.plot(x_cdf, empirical_cdf, 'r--', alpha=0.7, linewidth=2, label='经验CDF')
ax2.set_xlabel('x')
ax2.set_ylabel('累积概率')
ax2.set_title('累积分布函数比较')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 3. Q-Q图
from scipy.stats import probplot
probplot(chi2_sum, dist=stats.chi2(df_sum), plot=ax3)
ax3.set_title('Q-Q图: χ²(3) + χ²(7) vs χ²(10)')
ax3.grid(True, alpha=0.3)
# 4. 残差分析
# 比较χ²(3) + χ²(7)和χ²(10)直接生成的差异
differences = chi2_sum - chi2_direct
ax4.hist(differences, bins=50, alpha=0.7, color='purple', edgecolor='black')
ax4.axvline(np.mean(differences), color='red', linestyle='--', linewidth=2,
label=f'均值 = {np.mean(differences):.4f}')
ax4.set_xlabel('χ²(3) + χ²(7) - χ²(10)')
ax4.set_ylabel('频数')
ax4.set_title('两种方法的差异分布')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
print(f"\n结论:")
print("1. 三种方法生成的分布在统计上几乎相同")
print("2. 验证了卡方分布的可加性: χ²(df₁) + χ²(df₂) = χ²(df₁+df₂)")
print("3. 也验证了卡方分布的定义: k个独立标准正态分布的平方和")
chi2_additivity()
🔗 相关知识
- 正态分布
- 假设检验
- t分布
- Hadoop/HDFS分布式文件系统
- 拟合优度检验
- 独立性检验
- 方差
#统计学 #卡方分布 #概率分布 #连续分布 #基础