Python实战：8种统计分布可视化代码全解析（附避坑指南）

周晓农

172人浏览 · 2026-03-05 00:31:58

周晓农 · 2026-03-05 00:31:58 发布

Python实战：8种统计分布可视化代码全解析（附避坑指南）

很多刚接触数据科学的朋友，拿到一份数据，第一步往往就是画个直方图看看分布。但你是否想过，你看到的那个“形状”，背后可能对应着统计学中一个经典的数学模型？理解这些分布，不仅仅是应付考试，更是数据建模、假设检验、甚至机器学习特征工程的基础。然而，光知道理论公式，一到用Python画图时就容易掉坑：参数设错了图形诡异，分布选错了结论全偏。这篇文章，我就从一个实践者的角度，带你手把手用Python玩转八种核心统计分布的可视化。我们不只画图，更要讲清楚每个参数的实际意义、常见的使用误区，以及如何根据你的数据特点选择合适的分布进行拟合和检验。准备好了吗？让我们从最熟悉的钟形曲线开始。

1. 正态分布：无处不在的“标准模型”

提起统计分布，正态分布（也叫高斯分布）绝对是第一个蹦进脑海的。它描述了大量独立随机因素共同作用下的结果，从人的身高、考试分数到测量误差，随处可见它的身影。其概率密度函数（PDF）那个经典的钟形曲线，由两个参数决定：均值（μ）和标准差（σ）。均值决定了曲线的中心位置，标准差则控制了曲线的“胖瘦”。

在Python中，我们常用 numpy 生成随机样本，用 scipy.stats 获取理论分布，并用 matplotlib 进行可视化。一个常见的坑是，初学者容易混淆 np.random.normal 的 scale 参数（标准差）和 scipy.stats.norm 的 scale 参数（也是标准差），但有时其他函数里会用方差，这点需要特别注意。

下面是一个对比真实数据直方图与理论正态分布曲线的完整示例，其中包含了参数估计：

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

# 生成模拟数据：假设我们有一组用户完成某任务的时间数据（单位：秒）
np.random.seed(42)  # 确保结果可重现
true_mean = 180  # 真实均值：3分钟
true_std = 30    # 真实标准差：30秒
sample_data = np.random.normal(loc=true_mean, scale=true_std, size=1000)

# 计算样本的均值和标准差，作为对总体参数的估计
sample_mean = np.mean(sample_data)
sample_std = np.std(sample_data, ddof=1)  # 注意：ddof=1计算样本标准差

# 绘制直方图（密度形式）与理论PDF曲线
fig, ax = plt.subplots(figsize=(10, 6))
# 直方图
ax.hist(sample_data, bins=50, density=True, alpha=0.6, color='skyblue', edgecolor='black', label='样本直方图')
# 理论正态曲线
x_range = np.linspace(sample_data.min(), sample_data.max(), 1000)
pdf_values = stats.norm.pdf(x_range, loc=sample_mean, scale=sample_std)
ax.plot(x_range, pdf_values, 'r-', lw=3, label=f'正态分布拟合\n(μ={sample_mean:.1f}, σ={sample_std:.1f})')

ax.set_xlabel('任务完成时间 (秒)', fontsize=12)
ax.set_ylabel('概率密度', fontsize=12)
ax.set_title('用户任务完成时间分布与正态拟合', fontsize=14)
ax.legend()
ax.grid(True, linestyle='--', alpha=0.5)
plt.tight_layout()
plt.show()

# 进行正态性检验（Shapiro-Wilk检验）
stat, p_value = stats.shapiro(sample_data)
print(f"Shapiro-Wilk检验统计量: {stat:.4f}")
print(f"P值: {p_value:.4e}")
if p_value > 0.05:
    print("在0.05显著性水平下，无法拒绝数据来自正态分布的原假设。")
else:
    print("在0.05显著性水平下，数据不服从正态分布。")

注意：shapiro 检验在样本量很大时（如>5000），即使分布轻微偏离正态，也可能得到显著的P值。此时应结合Q-Q图等图形方法综合判断。

避坑指南1：参数混淆与图形解读

loc 与 scale：在 scipy.stats 的分布家族中，loc 通常是位置参数（如均值），scale 是尺度参数（如标准差）。但在 np.random.normal 中，参数名直接是 loc 和 scale。务必查阅文档，保持一致性。
直方图分箱（bins）：bins 数量选择不当会扭曲数据印象。太少会丢失细节，太多则会产生噪音。可以尝试 'auto'、'fd'（Freedman-Diaconis规则）或 'sturges' 等规则，并结合业务理解调整。
密度估计（density=True）：绘制直方图时设置 density=True 至关重要，它会让直方图的总面积归一化为1，从而能够与理论概率密度曲线进行直接、准确的比较。忘记设置会导致纵坐标尺度不同，比较失去意义。

2. 二项分布与泊松分布：离散事件的计数法则

当我们关心的是“成功次数”或“事件发生次数”时，就进入了离散分布的领域。其中，二项分布和泊松分布是最重要的两个代表。

二项分布描述的是在固定次数（n）的独立试验中，成功次数（k）的概率分布，每次试验的成功概率为p。比如，抛10次硬币得到正面的次数，或者100个用户中完成注册的人数。

泊松分布则常用于描述在固定时间或空间区间内，某个稀有事件发生次数的概率。它只有一个参数λ（lambda），表示单位区间内的平均发生次数。例如，一个客服中心每小时接到的电话数，或者一个网站每分钟的访问量。

两者的联系在于，当二项分布的n很大而p很小时，它可以由泊松分布近似（λ = n*p）。可视化离散分布时，我们通常使用茎叶图（stem plot）或条形图（bar plot），而不是连接数据点的折线图，以强调其离散性。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom, poisson

# 案例1：二项分布 - 新产品上线，历史点击转化率p=0.02，今日有n=500次独立曝光
n_trials = 500
p_success = 0.02
lambda_poisson = n_trials * p_success  # 泊松近似参数

# 计算二项分布和泊松分布的概率质量函数（PMF）
k_values = np.arange(0, 21)  # 观察成功次数从0到20
pmf_binom = binom.pmf(k_values, n_trials, p_success)
pmf_poisson_approx = poisson.pmf(k_values, lambda_poisson)

# 绘图对比
fig, axes = plt.subplots(1, 2, figsize=(14, 5))

# 二项分布图
axes[0].stem(k_values, pmf_binom, linefmt='b-', markerfmt='bo', basefmt=' ', label=f'二项分布 (n={n_trials}, p={p_success})')
axes[0].set_xlabel('成功次数 (k)', fontsize=12)
axes[0].set_ylabel('概率 P(X=k)', fontsize=12)
axes[0].set_title('二项分布概率质量函数', fontsize=14)
axes[0].legend()
axes[0].grid(True, linestyle='--', alpha=0.5)

# 泊松分布（近似）图
axes[1].stem(k_values, pmf_poisson_approx, linefmt='r-', markerfmt='rs', basefmt=' ', label=f'泊松分布 (λ={lambda_poisson:.1f})')
axes[1].set_xlabel('事件发生次数 (k)', fontsize=12)
axes[1].set_ylabel('概率 P(X=k)', fontsize=12)
axes[1].set_title('泊松分布概率质量函数（用于近似）', fontsize=14)
axes[1].legend()
axes[1].grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

# 案例2：模拟泊松过程 - 模拟一小时内网站访问量
np.random.seed(123)
hourly_visits = np.random.poisson(lam=15, size=30)  # 模拟30小时，每小时平均15次访问
print(f"模拟的30小时每小时访问量：\n{hourly_visits}")
print(f"平均访问量：{np.mean(hourly_visits):.2f}, 方差：{np.var(hourly_visits):.2f}")
# 泊松分布的特性：均值 ≈ 方差

避坑指南2：离散与连续的混淆

PMF vs PDF：对于离散分布（二项、泊松），我们计算和绘制的是概率质量函数（PMF），它给出的是某个具体值k的概率。对于连续分布（正态、指数），我们处理的是概率密度函数（PDF），其在某一点的值不是概率，概率需要通过积分在区间上获得。使用 scipy.stats 时，对应的方法是 .pmf() 和 .pdf()，千万别用错。
可视化选择：用 plt.plot() 画离散分布的PMF点，然后用直线连接，会误导观众认为中间值也有意义。plt.stem() 或 plt.bar() 是更专业的选择。
泊松分布的前提：事件的发生需要是独立的，且单位时间内的平均发生率（λ）是恒定的。如果数据表现出明显的周期性或趋势，则可能不适用泊松分布。

3. 均匀、指数与伽马分布：时间与等待的模型

这一组分布常与“时间”、“间隔”、“等待”相关。均匀分布描述的是在一个区间内所有结果等可能出现的场景，比如抽奖、随机采样。指数分布是泊松过程的“另一面”，它描述的是连续独立事件发生的时间间隔，比如两次客服来电的间隔时间。伽马分布则可以看作是多个独立指数分布变量之和的分布，用于描述直到第k个事件发生所需的等待时间。

理解它们的参数是关键：

均匀分布 Uniform(a, b)：a 是下限，b 是上限。
指数分布 Exp(λ)：λ 是速率参数（单位时间发生次数），均值是 1/λ。在 scipy.stats.expon 中，我们用 scale 参数，其值为 1/λ。
伽马分布 Gamma(α, β)：α (形状参数) 可以理解为事件发生的次数，β (尺度参数) 与指数分布的尺度参数意义相同。当 α=1 时，伽马分布退化为指数分布。

让我们通过一个模拟网站用户访问的场景，将三者联系起来：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform, expon, gamma

# 模拟场景：分析用户访问行为
np.random.seed(2024)

# 1. 均匀分布：模拟用户在一天24小时内随机开始访问的时间（小时）
visit_start_hour = uniform.rvs(loc=0, scale=24, size=1000)

# 2. 指数分布：模拟用户会话之间的间隔时间（分钟），假设平均间隔30分钟（λ=1/30）
mean_interval = 30  # 平均间隔30分钟
rate_lambda = 1.0 / mean_interval  # 速率参数
session_intervals = expon.rvs(scale=mean_interval, size=1000)  # scale = 1/λ

# 3. 伽马分布：模拟用户完成一个需要多个步骤（比如3步）的任务总耗时（分钟）
# 假设每个步骤耗时独立且服从指数分布（均值10分钟），总耗时即为形状参数α=3的伽马分布
step_mean_time = 10
shape_alpha = 3  # 步骤数
scale_beta = step_mean_time  # 尺度参数，这里等于指数分布的均值
task_total_time = gamma.rvs(a=shape_alpha, scale=scale_beta, size=1000)

# 绘制三个子图
fig, axes = plt.subplots(1, 3, figsize=(16, 5))

# 子图1：均匀分布 - 访问时间分布
axes[0].hist(visit_start_hour, bins=24, density=True, alpha=0.7, color='lightgreen', edgecolor='darkgreen')
x_unif = np.linspace(0, 24, 500)
pdf_unif = uniform.pdf(x_unif, loc=0, scale=24)
axes[0].plot(x_unif, pdf_unif, 'g-', lw=2, label='均匀分布PDF')
axes[0].set_xlabel('一天中的时间 (小时)', fontsize=12)
axes[0].set_ylabel('概率密度', fontsize=12)
axes[0].set_title('用户访问开始时间（均匀分布）', fontsize=14)
axes[0].legend()
axes[0].grid(True, linestyle='--', alpha=0.5)

# 子图2：指数分布 - 会话间隔时间
axes[1].hist(session_intervals, bins=50, density=True, alpha=0.7, color='lightcoral', edgecolor='darkred', label='模拟数据')
x_exp = np.linspace(0, session_intervals.max(), 500)
pdf_exp = expon.pdf(x_exp, scale=mean_interval)
axes[1].plot(x_exp, pdf_exp, 'r-', lw=2, label=f'指数分布PDF (λ={rate_lambda:.3f})')
axes[1].set_xlabel('会话间隔时间 (分钟)', fontsize=12)
axes[1].set_ylabel('概率密度', fontsize=12)
axes[1].set_title('用户会话间隔时间（指数分布）', fontsize=14)
axes[1].legend()
axes[1].grid(True, linestyle='--', alpha=0.5)

# 子图3：伽马分布 - 任务总耗时
axes[2].hist(task_total_time, bins=50, density=True, alpha=0.7, color='lightblue', edgecolor='darkblue', label='模拟数据')
x_gamma = np.linspace(0, task_total_time.max(), 500)
pdf_gamma = gamma.pdf(x_gamma, a=shape_alpha, scale=scale_beta)
axes[2].plot(x_gamma, pdf_gamma, 'b-', lw=2, label=f'伽马分布PDF (α={shape_alpha}, β={scale_beta})')
axes[2].set_xlabel('任务总耗时 (分钟)', fontsize=12)
axes[2].set_ylabel('概率密度', fontsize=12)
axes[2].set_title('多步骤任务总耗时（伽马分布）', fontsize=14)
axes[2].legend()
axes[2].grid(True, linestyle='--', alpha=0.5)

plt.tight_layout()
plt.show()

# 关键参数解读表格
print("\n分布参数与实际意义解读：")
print("| 分布类型 | 关键参数 | 参数含义 | 模拟案例中的意义 |")
print("| :--- | :--- | :--- | :--- |")
print(f"| 均匀分布 | a=0, b=24 | 取值范围下限与上限 | 一天中的小时数（0-24点） |")
print(f"| 指数分布 | λ={rate_lambda:.3f} | 事件发生速率 | 平均每{mean_interval}分钟发生一次会话 |")
print(f"| 伽马分布 | α={shape_alpha}, β={scale_beta} | 形状与尺度参数 | 完成{shape_alpha}个独立步骤，每步平均耗时{scale_beta}分钟 |")

避坑指南3：尺度参数与速率参数

指数分布的参数化：这是最大的混淆点。指数分布有两种常见参数化方式：速率参数λ（expon(scale=1/λ)）和尺度参数β（expon(scale=β)），其中 β = 1/λ。scipy.stats.expon 使用的是尺度参数β。务必明确你手中的参数是平均间隔时间（β）还是单位时间发生率（λ）。
伽马分布参数：scipy.stats.gamma 使用形状参数 a（即α）和尺度参数 scale（即β）。另一个常见的参数化是使用形状参数α和速率参数β（rate = 1/scale）。在调用函数前，确认你的参数定义。
数据范围：指数分布和伽马分布的定义域是 [0, +∞)。如果你的数据包含负值或零值，需要检查数据生成过程或考虑其他分布（如对数正态分布）。

4. t分布与F分布：统计推断的基石

t分布和F分布在经典的频率派统计推断中扮演着核心角色，它们通常不是直接对原始数据建模，而是作为检验统计量的抽样分布出现。

t分布：当我们用样本均值估计总体均值，并用样本标准差代替未知的总体标准差时，标准化后的统计量服从t分布。其形状类似正态分布，但尾部更厚，自由度（df）越小，尾部越厚。随着自由度增大，t分布趋近于标准正态分布。它主要用于单样本/双样本t检验、回归系数的显著性检验。
F分布：定义为两个独立的卡方分布变量除以各自自由度后的比值。它主要用于方差分析（ANOVA）、比较两个总体方差是否相等（F检验）、以及多元回归的整体显著性检验。

可视化这些分布，有助于我们理解假设检验中P值和临界值的由来。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import t, f, norm
import pandas as pd

# 比较不同自由度的t分布与标准正态分布
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))

# 子图1：t分布曲线族
x = np.linspace(-4, 4, 1000)
ax1.plot(x, norm.pdf(x), 'k--', lw=2, label='标准正态分布 (df=∞)')
for df in [1, 2, 5, 10]:
    ax1.plot(x, t.pdf(x, df), lw=2, label=f't分布 (df={df})')
ax1.set_xlabel('x', fontsize=12)
ax1.set_ylabel('概率密度', fontsize=12)
ax1.set_title('t分布与标准正态分布对比', fontsize=14)
ax1.legend()
ax1.grid(True, linestyle='--', alpha=0.5)
ax1.set_ylim(0, 0.42)

# 子图2：F分布曲线族（固定分母自由度）
x_f = np.linspace(0, 5, 1000)
dfd = 10  # 固定分母自由度
for dfn in [1, 2, 5, 10]:  # 变化分子自由度
    ax2.plot(x_f, f.pdf(x_f, dfn, dfd), lw=2, label=f'F分布 (dfn={dfn}, dfd={dfd})')
ax2.set_xlabel('F值', fontsize=12)
ax2.set_ylabel('概率密度', fontsize=12)
ax2.set_title('F分布（分母自由度dfd=10）', fontsize=14)
ax2.legend()
ax2.grid(True, linestyle='--', alpha=0.5)
ax2.set_ylim(0, 1.2)

plt.tight_layout()
plt.show()

# 实战演示：模拟一个双样本t检验，并可视化其t统计量的分布
np.random.seed(88)
# 模拟两组数据：A组（新方法），B组（旧方法）
group_a = np.random.normal(loc=105, scale=15, size=25)  # 均值105，标准差15，样本量25
group_b = np.random.normal(loc=100, scale=15, size=30)  # 均值100，标准差15，样本量30

# 执行独立双样本t检验（假设方差相等）
from scipy.stats import ttest_ind
t_stat, p_val = ttest_ind(group_a, group_b, equal_var=True)
print(f"\n双样本t检验结果：")
print(f"t统计量 = {t_stat:.4f}")
print(f"P值 = {p_val:.4f}")

# 在t分布上标出计算得到的t统计量
fig, ax = plt.subplots(figsize=(10, 6))
df = len(group_a) + len(group_b) - 2  # 自由度 = n1 + n2 - 2
x_t = np.linspace(-4, 4, 1000)
y_t = t.pdf(x_t, df)
ax.plot(x_t, y_t, 'b-', lw=2, label=f't分布 (df={df})')
# 绘制拒绝域（双侧检验，α=0.05）
alpha = 0.05
t_critical = t.ppf(1 - alpha/2, df)  # 临界值
x_fill_left = np.linspace(-4, -t_critical, 200)
x_fill_right = np.linspace(t_critical, 4, 200)
ax.fill_between(x_fill_left, t.pdf(x_fill_left, df), color='red', alpha=0.3, label=f'拒绝域 (α={alpha})')
ax.fill_between(x_fill_right, t.pdf(x_fill_right, df), color='red', alpha=0.3)
# 标记计算出的t统计量
ax.axvline(x=t_stat, color='green', linestyle='--', lw=2, label=f'观测t值 ({t_stat:.2f})')
ax.set_xlabel('t值', fontsize=12)
ax.set_ylabel('概率密度', fontsize=12)
ax.set_title(f'双样本t检验：t统计量抽样分布 (df={df})', fontsize=14)
ax.legend()
ax.grid(True, linestyle='--', alpha=0.5)
plt.show()

避坑指南4：检验的前提与自由度

t检验的前提：独立双样本t检验通常要求数据独立性、正态性（或大样本）和方差齐性。使用 ttest_ind 时，equal_var 参数需要根据方差是否相等来设置（可用Levene检验或F检验判断）。前提不满足可能导致结论错误。
自由度的计算：不同场景下t分布的自由度计算方式不同。对于单样本t检验，df = n - 1。对于独立双样本t检验（方差相等），df = n1 + n2 - 2。对于配对样本t检验，df = n - 1（n为对数）。错误的自度会影响P值的计算。
F检验的用途：F检验常用于方差分析（ANOVA）比较多个组均值，或检验两个总体方差是否相等（此时统计量是样本方差之比）。scipy.stats.f_oneway 用于前者，scipy.stats.levene 或 scipy.stats.bartlett 更稳健地用于方差齐性检验。明确你的检验目的，选择正确的函数。

5. 高级可视化与拟合优度评估

掌握了单个分布的可视化后，我们需要更进一步：如何判断一份数据到底服从哪种分布？以及如何同时对比多个分布的拟合效果？这里介绍两种强大的工具：Q-Q图和概率图。

Q-Q图（分位数-分位数图） 是一种直观比较数据分位数与理论分布分位数的图形方法。如果数据服从该理论分布，点会大致落在一条直线上。

概率图 是另一种变体，特别适用于检验数据是否服从特定的位置-尺度分布族（如正态、指数、威布尔等）。scipy.stats 提供了 probplot 函数来方便地绘制概率图。

让我们用一组模拟的偏态数据，来演示如何用这些工具评估不同分布的拟合优度。

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import warnings
warnings.filterwarnings('ignore')  # 忽略部分拟合警告

# 生成一组右偏（正偏）的模拟数据，例如网站用户停留时间（分钟）
np.random.seed(369)
# 使用伽马分布生成右偏数据
true_shape, true_scale = 2.5, 2.0
skewed_data = np.random.gamma(shape=true_shape, scale=true_scale, size=500)

# 尝试用三种分布去拟合：正态分布、指数分布、伽马分布
distributions_to_fit = [
    ('正态分布', stats.norm, None),  # None表示使用MLE自动估计参数
    ('指数分布', stats.expon, None),
    ('伽马分布', stats.gamma, (true_shape,))  # 固定形状参数？不，我们让MLE估计
]

# 创建图形
fig, axes = plt.subplots(2, 3, figsize=(16, 10))
axes = axes.flatten()

for idx, (dist_name, dist_obj, fit_params) in enumerate(distributions_to_fit):
    ax_hist = axes[idx]
    ax_qq = axes[idx + 3]

    # --- 子图1：直方图与拟合的PDF叠加 ---
    ax_hist.hist(skewed_data, bins=40, density=True, alpha=0.6, color='gray', edgecolor='black', label='数据直方图')
    # 拟合分布参数（使用最大似然估计MLE）
    if fit_params:
        params = dist_obj.fit(skewed_data, f0=fit_params[0])  # 固定部分参数
    else:
        params = dist_obj.fit(skewed_data)  # 完全MLE估计
    # 生成拟合的PDF曲线
    xmin, xmax = skewed_data.min(), skewed_data.max()
    x = np.linspace(xmin, xmax, 1000)
    pdf_fitted = dist_obj.pdf(x, *params)
    ax_hist.plot(x, pdf_fitted, 'r-', lw=3, label=f'{dist_name}拟合')
    ax_hist.set_xlabel('用户停留时间 (分钟)', fontsize=11)
    ax_hist.set_ylabel('概率密度', fontsize=11)
    ax_hist.set_title(f'{dist_name}拟合对比', fontsize=13)
    ax_hist.legend()
    ax_hist.grid(True, linestyle='--', alpha=0.5)

    # --- 子图2：Q-Q图 ---
    if dist_name == '伽马分布':
        # 对于伽马分布，probplot需要形状参数
        (osm, osr), (slope, intercept, r) = stats.probplot(skewed_data, dist=dist_name, sparams=params, plot=ax_qq)
        fit_label = f'Q-Q图 (形状={params[0]:.2f})'
    else:
        (osm, osr), (slope, intercept, r) = stats.probplot(skewed_data, dist=dist_name, plot=ax_qq)
        fit_label = f'Q-Q图'
    ax_qq.plot(osm, osm * slope + intercept, 'r-', lw=2, label=fit_label)
    ax_qq.set_xlabel('理论分位数', fontsize=11)
    ax_qq.set_ylabel('样本分位数', fontsize=11)
    ax_qq.set_title(f'{dist_name} Q-Q图', fontsize=13)
    ax_qq.legend()
    ax_qq.grid(True, linestyle='--', alpha=0.5)
    # 计算并显示R^2，衡量线性拟合程度
    ax_qq.text(0.05, 0.85, f'$R^2$ = {r**2:.4f}', transform=ax_qq.transAxes, fontsize=12,
               bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))

plt.tight_layout()
plt.show()

# 使用Kolmogorov-Smirnov检验进行定量评估
print("\nKolmogorov-Smirnov拟合优度检验结果（P值越大，越支持数据来自该分布）：")
ks_results = []
for dist_name, dist_obj, _ in distributions_to_fit:
    params = dist_obj.fit(skewed_data)
    # KS检验：比较经验分布函数与拟合的理论分布函数
    ks_stat, ks_pvalue = stats.kstest(skewed_data, dist_obj.cdf, args=params)
    ks_results.append((dist_name, ks_stat, ks_pvalue))
    print(f"{dist_name:10} KS统计量={ks_stat:.4f}, P值={ks_pvalue:.4f}")

# 找出P值最大的分布（最优拟合）
best_fit = max(ks_results, key=lambda x: x[2])
print(f"\n根据KS检验，最优拟合分布是：{best_fit[0]} (P值={best_fit[2]:.4f})")

避坑指南5：拟合优度检验的局限性

检验的效力：像KS检验这样的拟合优度检验，在大样本情况下非常敏感，即使分布只有轻微偏离也可能拒绝原假设。因此，不能仅凭P值小于0.05就断定分布不合适，必须结合图形（如Q-Q图、直方图拟合图）进行综合判断。
参数估计的影响：上述KS检验中，我们使用了基于数据估计的参数（MLE）。这实际上是一种“复合”假设检验，其临界值与简单假设检验不同。scipy.stats.kstest 在提供 args 参数时，其P值在某些情况下可能需要调整或使用更专门的检验（如Lilliefors检验用于正态性）。
多分布比较：当比较多个候选分布时，除了看拟合优度检验的P值，还可以比较赤池信息准则（AIC） 或贝叶斯信息准则（BIC），这些准则在考虑拟合优度的同时，惩罚了模型复杂度（参数数量）。对于嵌套模型（如指数分布是伽马分布的特例），还可以使用似然比检验。

6. 实战案例：从数据到分布选择的完整流程

理论讲得再多，不如一个完整的实战案例。假设你是一家电商公司的数据分析师，拿到了一份“用户首次购买前浏览商品页面次数”的数据。你的任务是：探索其分布，并选择一个合适的模型，用于后续的库存预测或用户行为模拟。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns

# 模拟业务数据：用户首次购买前的浏览页面数（显然是非负整数，且可能具有大量低值和长尾）
np.random.seed(999)
# 使用负二项分布（Negative Binomial）模拟，它常用于描述比泊松分布方差更大的计数数据。
# 负二项分布可以理解为，获得r次“成功”前，经历“失败”次数的分布。
# 这里我们用它模拟浏览页面数（“失败”）直到决定购买（“成功”）。
n = 1000  # 样本量
# 参数：r（成功次数，决定分布的峰位置），p（每次试验的成功概率，决定离散程度）
r, p = 5, 0.3
page_views_data = stats.nbinom.rvs(r, p, size=n)  # 生成数据

print("数据描述性统计：")
print(pd.Series(page_views_data).describe())
print(f"方差/均值比（离散指数）：{np.var(page_views_data)/np.mean(page_views_data):.3f}")
# 泊松分布方差等于均值，比值约等于1。若远大于1，则存在过离散（over-dispersion）。

# 步骤1：初步可视化 - 直方图与经验CDF
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 5))

# 直方图（注意计数数据用整数分箱）
ax1.hist(page_views_data, bins=range(0, page_views_data.max()+2, 2), density=True, alpha=0.7, color='steelblue', edgecolor='navy')
ax1.set_xlabel('浏览页面数', fontsize=12)
ax1.set_ylabel('概率密度（估计）', fontsize=12)
ax1.set_title('用户首次购买前浏览页面数分布（直方图）', fontsize=14)
ax1.grid(True, linestyle='--', alpha=0.5)

# 经验累积分布函数（ECDF）
sorted_data = np.sort(page_views_data)
y_ecdf = np.arange(1, len(sorted_data)+1) / len(sorted_data)
ax2.step(sorted_data, y_ecdf, where='post', lw=2, color='darkorange', label='经验CDF')
ax2.set_xlabel('浏览页面数', fontsize=12)
ax2.set_ylabel('累积概率', fontsize=12)
ax2.set_title('经验累积分布函数（ECDF）', fontsize=14)
ax2.grid(True, linestyle='--', alpha=0.5)
ax2.legend()

plt.tight_layout()
plt.show()

# 步骤2：尝试用几种常见离散分布拟合，并计算AIC/BIC
candidate_distributions = {
    '泊松分布': (stats.poisson, (), 'mle'),  # 一个参数λ
    '负二项分布': (stats.nbinom, (), 'mm'),   # 两个参数n, p，使用矩估计初始值更稳定
    '几何分布': (stats.geom, (), 'mle'),      # 负二项分布r=1时的特例
}

fit_results = []
for dist_name, (dist_obj, fit_args, fit_method) in candidate_distributions.items():
    try:
        # 拟合参数
        if dist_name == '负二项分布':
            # 负二项分布参数估计有时需要好的初始值，这里用矩估计
            mean = np.mean(page_views_data)
            var = np.var(page_views_data)
            p_est = mean / var if var > mean else 0.9  # 估计p
            n_est = mean * p_est / (1 - p_est) if p_est < 1 else 1  # 估计n
            params = (n_est, p_est)
        else:
            params = dist_obj.fit(page_views_data, *fit_args, method=fit_method)
        # 计算对数似然
        loglik = np.sum(dist_obj.logpmf(page_views_data, *params))
        # 计算AIC和BIC
        k = len(params)  # 参数个数
        n_obs = len(page_views_data)
        aic = 2 * k - 2 * loglik
        bic = k * np.log(n_obs) - 2 * loglik
        fit_results.append({
            '分布': dist_name,
            '参数': params,
            '对数似然': loglik,
            'AIC': aic,
            'BIC': bic
        })
    except Exception as e:
        print(f"拟合 {dist_name} 时出错: {e}")
        continue

# 展示拟合结果
df_results = pd.DataFrame(fit_results).sort_values('AIC')
print("\n候选分布拟合优度比较（AIC/BIC越小越好）：")
print(df_results.to_string(index=False))

# 步骤3：可视化最佳拟合分布
best_dist_name = df_results.iloc[0]['分布']
best_dist_info = candidate_distributions[best_dist_name]
best_dist_obj = best_dist_info[0]
best_params = df_results.iloc[0]['参数']

# 计算最佳分布的理论PMF
x_vals = np.arange(0, page_views_data.max() + 5)
if best_dist_name == '负二项分布':
    # scipy的nbinom参数化是 (n, p)，其中n是成功次数，p是成功概率
    pmf_best = best_dist_obj.pmf(x_vals, *best_params)
elif best_dist_name == '泊松分布':
    pmf_best = best_dist_obj.pmf(x_vals, *best_params)
else:
    pmf_best = best_dist_obj.pmf(x_vals, *best_params)

# 绘制对比图
fig, ax = plt.subplots(figsize=(10, 6))
# 绘制数据直方图（归一化）
ax.hist(page_views_data, bins=range(0, page_views_data.max()+2), density=True, alpha=0.7, color='lightgray', edgecolor='black', label='观测数据（直方图）')
# 绘制最佳拟合分布的PMF
ax.stem(x_vals, pmf_best, linefmt='r-', markerfmt='ro', basefmt=' ', label=f'最佳拟合：{best_dist_name}')
ax.set_xlabel('浏览页面数', fontsize=12)
ax.set_ylabel('概率', fontsize=12)
ax.set_title(f'用户浏览页面数分布与最佳拟合 ({best_dist_name})', fontsize=14)
ax.legend()
ax.grid(True, linestyle='--', alpha=0.5)
# 在图上标注参数
param_text = f"拟合参数: {best_params}"
ax.text(0.65, 0.95, param_text, transform=ax.transAxes, fontsize=11,
        verticalalignment='top', bbox=dict(boxstyle='round', facecolor='wheat', alpha=0.8))
plt.tight_layout()
plt.show()

print(f"\n业务解读与建议：")
print(f"1. 数据探索：浏览页面数均值为{np.mean(page_views_data):.1f}，方差为{np.var(page_views_data):.1f}，方差/均值比为{np.var(page_views_data)/np.mean(page_views_data):.2f}，明显大于1，存在'过离散'现象，简单的泊松分布可能不适用。")
print(f"2. 模型选择：根据AIC/BIC准则，{best_dist_name} 提供了最好的拟合。")
print(f"3. 应用方向：基于此分布模型，可以：")
print("   - 预测有多少比例的用户在浏览超过X个页面后才会购买。")
print("   - 模拟生成符合真实用户行为的浏览数据，用于A/B测试或系统压力测试。")
print("   - 识别异常用户（例如浏览页面数远超出模型预测范围的用户）。")

这个案例展示了从数据探索、描述性统计、图形化分析，到模型拟合、评估和选择的完整闭环。关键在于理解数据的业务背景（非负整数、计数过程），然后选择合适的候选分布族（这里是离散分布），最后使用统计准则（如AIC）和图形化工具（直方图叠加、Q-Q图）做出决策。记住，没有“唯一正确”的分布，只有“在当前业务目标和数据证据下更合适”的模型。