做数据分析这几年我越来越觉得在处理“看起来有问题”的数据时与其硬套单一高斯分布去拟合不如换一个更贴合实际的模型。比如某次我在分析一批传感器读数理论上应该是围绕一个标准值波动的结果画出来的直方图明显有两个峰一个集中在正常范围另一个偏到一侧。这种形态如果强行用单个高斯去描述拟合曲线会非常别扭后续的置信区间计算也会跟着失真。这种情况在统计里有个标准称呼双峰高斯分布也叫双峰高斯混合模型。处理这类数据时我通常的做法是先做一遍蒙特卡洛模拟把理论 PDF 和 CDF 画出来再对照真实样本验证。整套流程用 Python 的 NumPy、SciPy 和 Matplotlib 就能完成既能帮我们理解分布形态也能为后续的统计推断提供依据。这篇文章我会从建模思路、采样逻辑、绘图细节到常见坑位完整走一遍适合刚接触概率统计但想动手落地的人也适合日常需要跟异常分布打交道的分析人员和科研工作者。1. 双峰分布到底长什么样为什么不能当普通高斯处理1.1 混合权重、均值、标准差三个参数决定了整体形态双峰高斯分布本质上是两个高斯分布的加权叠加。数学上写成f(x) w1 * N(x; μ1, σ1²) (1 - w1) * N(x; μ2, σ2²)这里的 w1 是第一个峰的权重取值范围在 0 到 1 之间。μ1 和 μ2 是两个分量各自的中心位置σ1 和 σ2 是它们的离散程度。直观来看w1 越接近 0.5两个峰的高度差异越小μ1 和 μ2 离得越远两个峰分得越开σ 越大单个峰就越矮胖两个峰叠在一起后可能看起来就是一个又宽又平的分布。实操中最容易遇到的误判是把这种宽平分布当成方差很大的单一高斯分布但两者的尾部特征和分位数完全不同。举一个现实中很典型的例子网约车平台的订单量在全天时间轴上早晚高峰各出现一个高峰如果按小时统计订单数分布就是典型的双峰形态。还有一些场景比如产品售后赔款金额、生物实验中某种指标的测量值、金融资产收益率的分布都可能出现两个明显的中心。对这些数据建模时双峰高斯分布比单高斯有更强的描述能力。1.2 “硬算公式”和“蒙特卡洛模拟”各自的作用边界有人可能会问PDF 和 CDF 的公式都摆在那里直接用解析式画图不就行了吗为什么还要蒙特卡洛模拟这个疑问很合理。对于一维双峰高斯分布理论 PDF 和 CDF 的确可以直接用公式计算我在后文也会这么做。但蒙特卡洛模拟的价值在于当我们需要生成大量符合该分布的合成样本时只有通过随机抽样才能得到比如用来测试一个聚类算法、验证一个假设检验方法的功效或者做蒙特卡洛积分。当模型从一维扩展到多维、从两个分量扩展到多个分量或者引入截断、混合噪声之后解析 CDF 往往变得非常困难这时候抽样模拟几乎是唯一可操作的路径。即使只需要画 PDF 和 CDF模拟样本也可以用来验证自己的理论计算是否正确。把经验分布和理论曲线放在一起对比是自查代码最直接的方式。所以我的建议是把两者都做。理论曲线用来当基准模拟样本用来做实际验证和后续数据生成这样你在面对任何分布假设时都能有足够底气。2. 蒙特卡洛采样引擎两步抽样的核心逻辑2.1 混合权重如何决定样本来源从一个双峰高斯分布中生成随机样本逻辑上可以分两步走第一步为每个样本随机选择“它来源于哪个峰”。这一步用均匀分布随机数实现生成一个服从 U(0,1) 的随机数 u如果 u w1就认为样本来自第一个高斯分量否则来自第二个。第二步在该分量内部用该分量的均值和标准差生成正态分布样本。这个方法在统计里叫“成分混合抽样”实现非常简单而且时间效率很高。需要注意的一个细节是如果先选择了成分再逐个生成样本最后样本的顺序天然就会按成分“分堆”假如后续要做可视化或混淆排序最好在生成后随机打乱。2.2 完整采样代码与随机种子策略下面这一段是我在实际项目中反复使用的采样函数参数放在字典里方便调整和记录实验条件import numpy as np from scipy import stats import matplotlib.pyplot as plt # 全局随机数生成器保证结果可复现 rng np.random.default_rng(42) # 双峰高斯分布的参数 params { w1: 0.35, mu1: -1.0, sigma1: 0.8, mu2: 3.5, sigma2: 1.6, } def sample_two_gauss(n, params, rng): w1 params[w1] mu1 params[mu1] sigma1 params[sigma1] mu2 params[mu2] sigma2 params[sigma2] # 第一步按权重决定每个样本属于哪个分量 u rng.random(n) n1 int(np.sum(u w1)) # 第二步分别从两个正态分布中抽样 samples np.empty(n) samples[:n1] rng.normal(mu1, sigma1, n1) samples[n1:] rng.normal(mu2, sigma2, n - n1) # 打乱顺序避免样本按分量“排好队” rng.shuffle(samples) return samples这里我建议所有用到随机数的地方统一使用np.random.default_rng(seed)而不是旧的np.random.seed()。前者的优势是生成器对象可以显式传入不同的实验模块不会因为某个模块调用了全局随机函数而影响全局状态。我在跑蒙特卡洛实验时习惯给每类实验分配独立的 rng 实例这样即使某段代码出现随机数消耗不一致的问题也不会污染其他部分的复现结果。2.3 样本量和随机种子对结果的影响你可能已经注意到了我直接生成了 5 万个样本。为什么我不先用 500 个看看效果因为 PDF 绘图对样本量的敏感度很高。n 500 时直方图的每个 bin 里可能只有个位数甚至零个样本密度曲线凹凸不平n 5000 时主体轮廓能看出来但双峰的边界位置仍然会抖动n 50000 时经验分布和理论曲线基本能重合到视觉上很难区分。做展示、写报告或验证算法时我一般直接用 50000 到 100000 个样本。如果只是快速验证代码逻辑1000 个也够用但别指望它画出来的直方图多光滑。随机种子也同样关键。固定住种子之后同样代码跑出来的样本完全一致这个特性在写论文、复现他人实验和排查 bug 时极为重要。我在每次代码文件的头部都写清“seed42”这已经成了我的个人习惯。3. 理论 PDF 和 CDF 的计算不能只靠直方图3.1 用 SciPy 的 norm 分布组装理论曲线采样完之后下一步是画理论 PDF。这里直接用scipy.stats.norm分别计算两个高斯分量的密度再按权重相加。计算网格的范围要覆盖到两个分量的主要区间一般取 μ1-4σ1 到 μ24σ2这样曲线的两端都压到接近零的幅度图看起来更完整。# 计算用于绘图的网格 x np.linspace(-4.5, 9.5, 1000) # 理论PDF两个高斯密度加权相加 pdf_theory ( params[w1] * stats.norm.pdf(x, params[mu1], params[sigma1]) (1 - params[w1]) * stats.norm.pdf(x, params[mu2], params[sigma2]) ) # 理论CDF对应两个正态CDF的加权相加 cdf_theory ( params[w1] * stats.norm.cdf(x, params[mu1], params[sigma1]) (1 - params[w1]) * stats.norm.cdf(x, params[mu2], params[sigma2]) )这里有一个容易出错的点PDF 画出来是曲线CDF 画出来也是曲线但两者的量纲完全不同。PDF 的纵轴是“概率密度”可以大于 1CDF 的纵轴是“累积概率”永远在 0 到 1 之间。如果把它们画在同一个坐标系里CDF 几乎会贴着 x 轴形同一条水平线所以要么分两张图画要么用双 y 轴后面我会详细说。3.2 经验 CDF 的两种实现方法和经验 PDF也就是直方图不同经验 CDF 不需要分箱因此没有 bin 数量这个超参数稳定性更好。它的计算方式非常直观把样本从小到大排序在第 i 个数据点处累积概率就是 i/n。最简单的实现方式sorted_samples np.sort(samples) n_samples len(sorted_samples) ecdf_y np.arange(1, n_samples 1) / n_samples画图时用step而不是直接用散点或折线因为真实 CDF 是阶梯函数。用wherepost表示阶梯从当前 x 开始一直保持到下一个观测值之前这更贴近统计里“右连续”的约定。另一种更严谨的实现方法是使用np.searchsorted统计每个取值之前的样本比例。对于大量重复数据的场景这种写法更准确def empirical_cdf(data, grid): sorted_data np.sort(data) return np.searchsorted(sorted_data, grid, sideright) / len(data)这种方式的好处是可以在任意网格上计算经验 CDF方便和理论 CDF 在同一条 x 轴上比较也方便后续计算 KS 检验里的最大距离。4. 绘图阶段最容易踩的坑从默认样式到能直接放进论文的图4.1 直方图和理论密度怎么叠加才不误导人画 PDF 对比图时我最常看到的问题是直方图的纵轴忘了设置成密度。默认情况下Matplotlib 的hist返回的是频数如果样本量很大、bin 比较宽直方图纵轴数值会变得很大直接和理论 PDF 叠在一起曲线会被压成一条几乎水平的线。所以必须加上一个关键参数fig, ax plt.subplots(figsize(8, 5)) ax.hist(samples, bins80, densityTrue, alpha0.55, label模拟样本直方图) ax.plot(x, pdf_theory, r-, lw2, label理论PDF) ax.set_xlabel(x) ax.set_ylabel(概率密度) ax.legend() ax.set_title(双峰高斯分布模拟直方图与理论PDF)densityTrue会把直方图归一化使得所有 bin 的面积之和为 1这样才能和理论概率密度曲线放在同一个尺度下比较。另一个常见问题是 bin 数量选得太少或太多。选太少双峰可能被“糊”成一个宽峰选太多每个 bin 的随机起伏太大。对于 5 万样本60 到 100 个 bin 是经验安全区间。如果样本量较小可以参考“平均每个 bin 至少 10 个样本”的原则来倒推 bin 数量。4.2 中文标签、坐标刻度与科研绘图风格Matplotlib 的中文支持是老生常谈的问题。默认字体里没有中文字体直接写set_xlabel(样本值)会显示成方框。解决办法是plt.rcParams[font.sans-serif] [SimHei, Noto Sans CJK SC, Microsoft YaHei] plt.rcParams[axes.unicode_minus] False第二行的用处往往被忽略设置了中文字体后坐标轴上的负号有时会显示成乱码因为默认的 Unicode 负号和中文字体不兼容。加上这一行负号就能正常显示。如果要把图放进论文或实验报告我还会进一步微调用figsize(8, 5)控制比例避免默认图形太扁用savefig(pdf_cdf.png, dpi200, bbox_inchestight)保存保证四周不留白将线条颜色改成色盲友好的配色比如深红、深蓝、纯黑图例位置显式指定避免默认位置盖住曲线。4.3 双 y 轴处理 PDF 和 CDF 的组合图如果你想在一张图里同时展示 PDF 和 CDF双 y 轴是最好的方式。PDF 用左侧 y 轴CDF 用右侧 y 轴两者互不干扰fig, ax1 plt.subplots(figsize(8, 5)) ax1.hist(samples, bins80, densityTrue, alpha0.55, label模拟样本直方图) ax1.plot(x, pdf_theory, r-, lw2, label理论PDF) ax1.set_xlabel(x) ax1.set_ylabel(概率密度) ax2 ax1.twinx() ax2.step(sorted_samples, ecdf_y, wherepost, colork, lw1.0, label经验CDF) ax2.plot(x, cdf_theory, b--, lw1.5, label理论CDF) ax2.set_ylabel(累积概率) # 合并两个坐标轴的图例 lines1, labels1 ax1.get_legend_handles_labels() lines2, labels2 ax2.get_legend_handles_labels() ax1.legend(lines1 lines2, labels1 labels2, locbest) plt.tight_layout() plt.savefig(two_gauss_pdf_cdf.png, dpi200, bbox_inchestight)这里要注意twinx()创建的第二个坐标轴只是共用 x 轴y 轴范围是独立的所以即使 PDF 峰值大于 1CDF 也依然能正常显示。5. 参数没选好时图会骗人三个典型陷阱与校准方法5.1 两个峰靠太近时视觉上根本看不出“双峰”如果 μ1 和 μ2 的距离小于 σ1 σ2 的好几倍两个峰叠加后中间不会出现明显谷底整体看起来就是一个扁平的“胖高斯”。这种情况下直接从直方图判断分布形态很容易漏判。我处理这类数据时不会只依赖直方图而是同时观察 CDF 的形态。单高斯分布的 CDF 是一条平滑 S 形曲线而双峰分布尤其是两峰分离较远时CDF 会在中间出现一个相对平缓的“台阶”。所谓台阶就是累积概率在某个区间内增长明显变慢甚至停滞这是双峰结构非常典型的信号。所以画图时我总会把 PDF 和 CDF 一起输出PDF 负责直观展示峰的位置CDF 负责暴露 PDF 中不明显的变化——这是一种很有效的“互查机制”。5.2 权重失衡时小峰容易被淹没当 w1 特别小比如 0.05第一个分量在直方图里可能只会留下一个小小的隆起。这时就算 μ1 和 μ2 离得足够远直方图也可能因为 bin 宽度不合适而看不出这个峰。解决办法有两个。一个是增加样本量让弱势峰积累足够的样本另一个是调整 bin 宽度采用更细的 bin 或者改用核密度估计KDE来补充视觉效果。SciPy 里可以这样快速生成 KDE 曲线kde stats.gaussian_kde(samples, bw_method0.3) pdf_kde kde(x)带宽bw_method越小曲线越细致但噪声也越大。对弱势峰来说0.2 到 0.4 的带宽通常能暴露更多细节。5.3 用最小二乘拟合反推混合参数模拟完了之后有时候还需要用模拟数据反推真实参数。比如你确实有一批真实观测数据猜测它来自双峰高斯分布想验证一下并恢复出两个峰的中心位置和权重。最常用的办法是拟合直方图。from scipy.optimize import curve_fit def two_gauss_pdf(x, w1, mu1, sigma1, mu2, sigma2): return ( w1 * stats.norm.pdf(x, mu1, sigma1) (1 - w1) * stats.norm.pdf(x, mu2, sigma2) ) hist, bin_edges np.histogram(samples, bins80, densityTrue) bin_centers (bin_edges[:-1] bin_edges[1:]) / 2 popt, pcov curve_fit( two_gauss_pdf, bin_centers, hist, p0[0.5, 0.0, 1.0, 3.0, 1.0], bounds([0, -10, 0.01, -10, 0.01], [1, 10, 10, 10, 10]), )拟合完成后把popt里的参数代入理论 PDF和原始直方图叠加对比如果曲线贴合度不错说明双峰高斯的假设是合理的。实操中要特别留意初值p0的设定因为混合模型的似然面经常有多个局部最优解初值离真实值太远容易拟合出负权重或完全错误的峰位置。下表总结了我在调参过程中遇到的现象、可能原因和应对方法现象可能原因常用对策直方图只有一个明显峰两峰间距太小改用 CDF 或 KDE 观察小峰时有时无样本量不足或权重过低增加样本量调细 bin理论 PDF 与直方图整体偏高忘记加 densityTrue直方图归一化CDF 出现“平台阶”两峰分离明显加强双峰模型假设检验拟合参数偏离预期初值太差用分位点估计初值6. 双峰高斯模拟在实际项目里的三种用法6.1 从 CDF 快速反查区间概率模拟并画出 CDF 后我们发现它不仅能直观展示分布形态还能用于概率查询。比如想知道样本落在 0 到 2 之间的概率理论 CDF 直接相减就能得到p_range stats.norm.cdf(2, params[mu1], params[sigma1]) - stats.norm.cdf(0, params[mu1], params[sigma1]) print(p_range)但如果分布的方差或偏态复杂手动用公式不保险最稳妥的办法还是基于大样本模拟做分位数估算。用np.percentile(samples, [2.5, 50, 97.5])可以快速算出 95% 覆盖区间。这个方法在业务里经常用来做基准区间判断比如判断当前指标是否明显偏离正常波动范围。6.2 生成合成数据集做算法验证双峰高斯分布是制造“分类困难样本”的天然工具。当我在测试一个聚类算法时如果样本是两个重叠的高斯簇恰好能检验算法能否正确把混合区域里的样本分开。生成这样的合成数据非常简单前面定义的sample_two_gauss直接就能用而且可以附加标签把来自分量的 0 的样本标为类别 0来自分量的 1 的样本标为类别 1。这种做法在论文实验和模型评估中很实用因为真实标注数据往往稀少而合成数据可以控制簇间距、权重、噪声水平从而系统测试算法的鲁棒性。6.3 后续扩展多维、多峰与核密度对比这篇文章只展示了一维双峰的情况但现实中数据往往更复杂。如果你想继续扩展有三个方向很自然把两个高斯分量替换成多维高斯采样时用rng.multivariate_normalPDF 用scipy.stats.multivariate_normal.pdf把双峰扩展到三峰甚至更多混合权重变成一个概率向量仍然沿用“先按权重选分量再从分量抽样”的思路用核密度估计对比理论 PDF检查模型假设和真实样本的差异。多维情况的蒙特卡洛模拟码会复杂一些但核心理念完全一致先确定模型结构再设计抽样步骤最后用理论和经验分布互相验证。这也是我为什么一直强调要把一维向基础流程彻底搞明白的原因后面的扩展都是在这套流程上做加法。最后分享一个我的习惯每改动一次参数就把参数字典、随机种子、样本量记在代码注释里。双峰分布看起来只是多了三个参数但一旦开始调权重、调均值、调标准差组合爆炸的速度非常快不做好实验记录很容易忘了当前图对应的参数版本。把这些基础工作做扎实整个模拟流程才能真正成为后面分析工作的可靠支撑。
双峰高斯分布与蒙特卡洛模拟：Python实现指南
做数据分析这几年我越来越觉得在处理“看起来有问题”的数据时与其硬套单一高斯分布去拟合不如换一个更贴合实际的模型。比如某次我在分析一批传感器读数理论上应该是围绕一个标准值波动的结果画出来的直方图明显有两个峰一个集中在正常范围另一个偏到一侧。这种形态如果强行用单个高斯去描述拟合曲线会非常别扭后续的置信区间计算也会跟着失真。这种情况在统计里有个标准称呼双峰高斯分布也叫双峰高斯混合模型。处理这类数据时我通常的做法是先做一遍蒙特卡洛模拟把理论 PDF 和 CDF 画出来再对照真实样本验证。整套流程用 Python 的 NumPy、SciPy 和 Matplotlib 就能完成既能帮我们理解分布形态也能为后续的统计推断提供依据。这篇文章我会从建模思路、采样逻辑、绘图细节到常见坑位完整走一遍适合刚接触概率统计但想动手落地的人也适合日常需要跟异常分布打交道的分析人员和科研工作者。1. 双峰分布到底长什么样为什么不能当普通高斯处理1.1 混合权重、均值、标准差三个参数决定了整体形态双峰高斯分布本质上是两个高斯分布的加权叠加。数学上写成f(x) w1 * N(x; μ1, σ1²) (1 - w1) * N(x; μ2, σ2²)这里的 w1 是第一个峰的权重取值范围在 0 到 1 之间。μ1 和 μ2 是两个分量各自的中心位置σ1 和 σ2 是它们的离散程度。直观来看w1 越接近 0.5两个峰的高度差异越小μ1 和 μ2 离得越远两个峰分得越开σ 越大单个峰就越矮胖两个峰叠在一起后可能看起来就是一个又宽又平的分布。实操中最容易遇到的误判是把这种宽平分布当成方差很大的单一高斯分布但两者的尾部特征和分位数完全不同。举一个现实中很典型的例子网约车平台的订单量在全天时间轴上早晚高峰各出现一个高峰如果按小时统计订单数分布就是典型的双峰形态。还有一些场景比如产品售后赔款金额、生物实验中某种指标的测量值、金融资产收益率的分布都可能出现两个明显的中心。对这些数据建模时双峰高斯分布比单高斯有更强的描述能力。1.2 “硬算公式”和“蒙特卡洛模拟”各自的作用边界有人可能会问PDF 和 CDF 的公式都摆在那里直接用解析式画图不就行了吗为什么还要蒙特卡洛模拟这个疑问很合理。对于一维双峰高斯分布理论 PDF 和 CDF 的确可以直接用公式计算我在后文也会这么做。但蒙特卡洛模拟的价值在于当我们需要生成大量符合该分布的合成样本时只有通过随机抽样才能得到比如用来测试一个聚类算法、验证一个假设检验方法的功效或者做蒙特卡洛积分。当模型从一维扩展到多维、从两个分量扩展到多个分量或者引入截断、混合噪声之后解析 CDF 往往变得非常困难这时候抽样模拟几乎是唯一可操作的路径。即使只需要画 PDF 和 CDF模拟样本也可以用来验证自己的理论计算是否正确。把经验分布和理论曲线放在一起对比是自查代码最直接的方式。所以我的建议是把两者都做。理论曲线用来当基准模拟样本用来做实际验证和后续数据生成这样你在面对任何分布假设时都能有足够底气。2. 蒙特卡洛采样引擎两步抽样的核心逻辑2.1 混合权重如何决定样本来源从一个双峰高斯分布中生成随机样本逻辑上可以分两步走第一步为每个样本随机选择“它来源于哪个峰”。这一步用均匀分布随机数实现生成一个服从 U(0,1) 的随机数 u如果 u w1就认为样本来自第一个高斯分量否则来自第二个。第二步在该分量内部用该分量的均值和标准差生成正态分布样本。这个方法在统计里叫“成分混合抽样”实现非常简单而且时间效率很高。需要注意的一个细节是如果先选择了成分再逐个生成样本最后样本的顺序天然就会按成分“分堆”假如后续要做可视化或混淆排序最好在生成后随机打乱。2.2 完整采样代码与随机种子策略下面这一段是我在实际项目中反复使用的采样函数参数放在字典里方便调整和记录实验条件import numpy as np from scipy import stats import matplotlib.pyplot as plt # 全局随机数生成器保证结果可复现 rng np.random.default_rng(42) # 双峰高斯分布的参数 params { w1: 0.35, mu1: -1.0, sigma1: 0.8, mu2: 3.5, sigma2: 1.6, } def sample_two_gauss(n, params, rng): w1 params[w1] mu1 params[mu1] sigma1 params[sigma1] mu2 params[mu2] sigma2 params[sigma2] # 第一步按权重决定每个样本属于哪个分量 u rng.random(n) n1 int(np.sum(u w1)) # 第二步分别从两个正态分布中抽样 samples np.empty(n) samples[:n1] rng.normal(mu1, sigma1, n1) samples[n1:] rng.normal(mu2, sigma2, n - n1) # 打乱顺序避免样本按分量“排好队” rng.shuffle(samples) return samples这里我建议所有用到随机数的地方统一使用np.random.default_rng(seed)而不是旧的np.random.seed()。前者的优势是生成器对象可以显式传入不同的实验模块不会因为某个模块调用了全局随机函数而影响全局状态。我在跑蒙特卡洛实验时习惯给每类实验分配独立的 rng 实例这样即使某段代码出现随机数消耗不一致的问题也不会污染其他部分的复现结果。2.3 样本量和随机种子对结果的影响你可能已经注意到了我直接生成了 5 万个样本。为什么我不先用 500 个看看效果因为 PDF 绘图对样本量的敏感度很高。n 500 时直方图的每个 bin 里可能只有个位数甚至零个样本密度曲线凹凸不平n 5000 时主体轮廓能看出来但双峰的边界位置仍然会抖动n 50000 时经验分布和理论曲线基本能重合到视觉上很难区分。做展示、写报告或验证算法时我一般直接用 50000 到 100000 个样本。如果只是快速验证代码逻辑1000 个也够用但别指望它画出来的直方图多光滑。随机种子也同样关键。固定住种子之后同样代码跑出来的样本完全一致这个特性在写论文、复现他人实验和排查 bug 时极为重要。我在每次代码文件的头部都写清“seed42”这已经成了我的个人习惯。3. 理论 PDF 和 CDF 的计算不能只靠直方图3.1 用 SciPy 的 norm 分布组装理论曲线采样完之后下一步是画理论 PDF。这里直接用scipy.stats.norm分别计算两个高斯分量的密度再按权重相加。计算网格的范围要覆盖到两个分量的主要区间一般取 μ1-4σ1 到 μ24σ2这样曲线的两端都压到接近零的幅度图看起来更完整。# 计算用于绘图的网格 x np.linspace(-4.5, 9.5, 1000) # 理论PDF两个高斯密度加权相加 pdf_theory ( params[w1] * stats.norm.pdf(x, params[mu1], params[sigma1]) (1 - params[w1]) * stats.norm.pdf(x, params[mu2], params[sigma2]) ) # 理论CDF对应两个正态CDF的加权相加 cdf_theory ( params[w1] * stats.norm.cdf(x, params[mu1], params[sigma1]) (1 - params[w1]) * stats.norm.cdf(x, params[mu2], params[sigma2]) )这里有一个容易出错的点PDF 画出来是曲线CDF 画出来也是曲线但两者的量纲完全不同。PDF 的纵轴是“概率密度”可以大于 1CDF 的纵轴是“累积概率”永远在 0 到 1 之间。如果把它们画在同一个坐标系里CDF 几乎会贴着 x 轴形同一条水平线所以要么分两张图画要么用双 y 轴后面我会详细说。3.2 经验 CDF 的两种实现方法和经验 PDF也就是直方图不同经验 CDF 不需要分箱因此没有 bin 数量这个超参数稳定性更好。它的计算方式非常直观把样本从小到大排序在第 i 个数据点处累积概率就是 i/n。最简单的实现方式sorted_samples np.sort(samples) n_samples len(sorted_samples) ecdf_y np.arange(1, n_samples 1) / n_samples画图时用step而不是直接用散点或折线因为真实 CDF 是阶梯函数。用wherepost表示阶梯从当前 x 开始一直保持到下一个观测值之前这更贴近统计里“右连续”的约定。另一种更严谨的实现方法是使用np.searchsorted统计每个取值之前的样本比例。对于大量重复数据的场景这种写法更准确def empirical_cdf(data, grid): sorted_data np.sort(data) return np.searchsorted(sorted_data, grid, sideright) / len(data)这种方式的好处是可以在任意网格上计算经验 CDF方便和理论 CDF 在同一条 x 轴上比较也方便后续计算 KS 检验里的最大距离。4. 绘图阶段最容易踩的坑从默认样式到能直接放进论文的图4.1 直方图和理论密度怎么叠加才不误导人画 PDF 对比图时我最常看到的问题是直方图的纵轴忘了设置成密度。默认情况下Matplotlib 的hist返回的是频数如果样本量很大、bin 比较宽直方图纵轴数值会变得很大直接和理论 PDF 叠在一起曲线会被压成一条几乎水平的线。所以必须加上一个关键参数fig, ax plt.subplots(figsize(8, 5)) ax.hist(samples, bins80, densityTrue, alpha0.55, label模拟样本直方图) ax.plot(x, pdf_theory, r-, lw2, label理论PDF) ax.set_xlabel(x) ax.set_ylabel(概率密度) ax.legend() ax.set_title(双峰高斯分布模拟直方图与理论PDF)densityTrue会把直方图归一化使得所有 bin 的面积之和为 1这样才能和理论概率密度曲线放在同一个尺度下比较。另一个常见问题是 bin 数量选得太少或太多。选太少双峰可能被“糊”成一个宽峰选太多每个 bin 的随机起伏太大。对于 5 万样本60 到 100 个 bin 是经验安全区间。如果样本量较小可以参考“平均每个 bin 至少 10 个样本”的原则来倒推 bin 数量。4.2 中文标签、坐标刻度与科研绘图风格Matplotlib 的中文支持是老生常谈的问题。默认字体里没有中文字体直接写set_xlabel(样本值)会显示成方框。解决办法是plt.rcParams[font.sans-serif] [SimHei, Noto Sans CJK SC, Microsoft YaHei] plt.rcParams[axes.unicode_minus] False第二行的用处往往被忽略设置了中文字体后坐标轴上的负号有时会显示成乱码因为默认的 Unicode 负号和中文字体不兼容。加上这一行负号就能正常显示。如果要把图放进论文或实验报告我还会进一步微调用figsize(8, 5)控制比例避免默认图形太扁用savefig(pdf_cdf.png, dpi200, bbox_inchestight)保存保证四周不留白将线条颜色改成色盲友好的配色比如深红、深蓝、纯黑图例位置显式指定避免默认位置盖住曲线。4.3 双 y 轴处理 PDF 和 CDF 的组合图如果你想在一张图里同时展示 PDF 和 CDF双 y 轴是最好的方式。PDF 用左侧 y 轴CDF 用右侧 y 轴两者互不干扰fig, ax1 plt.subplots(figsize(8, 5)) ax1.hist(samples, bins80, densityTrue, alpha0.55, label模拟样本直方图) ax1.plot(x, pdf_theory, r-, lw2, label理论PDF) ax1.set_xlabel(x) ax1.set_ylabel(概率密度) ax2 ax1.twinx() ax2.step(sorted_samples, ecdf_y, wherepost, colork, lw1.0, label经验CDF) ax2.plot(x, cdf_theory, b--, lw1.5, label理论CDF) ax2.set_ylabel(累积概率) # 合并两个坐标轴的图例 lines1, labels1 ax1.get_legend_handles_labels() lines2, labels2 ax2.get_legend_handles_labels() ax1.legend(lines1 lines2, labels1 labels2, locbest) plt.tight_layout() plt.savefig(two_gauss_pdf_cdf.png, dpi200, bbox_inchestight)这里要注意twinx()创建的第二个坐标轴只是共用 x 轴y 轴范围是独立的所以即使 PDF 峰值大于 1CDF 也依然能正常显示。5. 参数没选好时图会骗人三个典型陷阱与校准方法5.1 两个峰靠太近时视觉上根本看不出“双峰”如果 μ1 和 μ2 的距离小于 σ1 σ2 的好几倍两个峰叠加后中间不会出现明显谷底整体看起来就是一个扁平的“胖高斯”。这种情况下直接从直方图判断分布形态很容易漏判。我处理这类数据时不会只依赖直方图而是同时观察 CDF 的形态。单高斯分布的 CDF 是一条平滑 S 形曲线而双峰分布尤其是两峰分离较远时CDF 会在中间出现一个相对平缓的“台阶”。所谓台阶就是累积概率在某个区间内增长明显变慢甚至停滞这是双峰结构非常典型的信号。所以画图时我总会把 PDF 和 CDF 一起输出PDF 负责直观展示峰的位置CDF 负责暴露 PDF 中不明显的变化——这是一种很有效的“互查机制”。5.2 权重失衡时小峰容易被淹没当 w1 特别小比如 0.05第一个分量在直方图里可能只会留下一个小小的隆起。这时就算 μ1 和 μ2 离得足够远直方图也可能因为 bin 宽度不合适而看不出这个峰。解决办法有两个。一个是增加样本量让弱势峰积累足够的样本另一个是调整 bin 宽度采用更细的 bin 或者改用核密度估计KDE来补充视觉效果。SciPy 里可以这样快速生成 KDE 曲线kde stats.gaussian_kde(samples, bw_method0.3) pdf_kde kde(x)带宽bw_method越小曲线越细致但噪声也越大。对弱势峰来说0.2 到 0.4 的带宽通常能暴露更多细节。5.3 用最小二乘拟合反推混合参数模拟完了之后有时候还需要用模拟数据反推真实参数。比如你确实有一批真实观测数据猜测它来自双峰高斯分布想验证一下并恢复出两个峰的中心位置和权重。最常用的办法是拟合直方图。from scipy.optimize import curve_fit def two_gauss_pdf(x, w1, mu1, sigma1, mu2, sigma2): return ( w1 * stats.norm.pdf(x, mu1, sigma1) (1 - w1) * stats.norm.pdf(x, mu2, sigma2) ) hist, bin_edges np.histogram(samples, bins80, densityTrue) bin_centers (bin_edges[:-1] bin_edges[1:]) / 2 popt, pcov curve_fit( two_gauss_pdf, bin_centers, hist, p0[0.5, 0.0, 1.0, 3.0, 1.0], bounds([0, -10, 0.01, -10, 0.01], [1, 10, 10, 10, 10]), )拟合完成后把popt里的参数代入理论 PDF和原始直方图叠加对比如果曲线贴合度不错说明双峰高斯的假设是合理的。实操中要特别留意初值p0的设定因为混合模型的似然面经常有多个局部最优解初值离真实值太远容易拟合出负权重或完全错误的峰位置。下表总结了我在调参过程中遇到的现象、可能原因和应对方法现象可能原因常用对策直方图只有一个明显峰两峰间距太小改用 CDF 或 KDE 观察小峰时有时无样本量不足或权重过低增加样本量调细 bin理论 PDF 与直方图整体偏高忘记加 densityTrue直方图归一化CDF 出现“平台阶”两峰分离明显加强双峰模型假设检验拟合参数偏离预期初值太差用分位点估计初值6. 双峰高斯模拟在实际项目里的三种用法6.1 从 CDF 快速反查区间概率模拟并画出 CDF 后我们发现它不仅能直观展示分布形态还能用于概率查询。比如想知道样本落在 0 到 2 之间的概率理论 CDF 直接相减就能得到p_range stats.norm.cdf(2, params[mu1], params[sigma1]) - stats.norm.cdf(0, params[mu1], params[sigma1]) print(p_range)但如果分布的方差或偏态复杂手动用公式不保险最稳妥的办法还是基于大样本模拟做分位数估算。用np.percentile(samples, [2.5, 50, 97.5])可以快速算出 95% 覆盖区间。这个方法在业务里经常用来做基准区间判断比如判断当前指标是否明显偏离正常波动范围。6.2 生成合成数据集做算法验证双峰高斯分布是制造“分类困难样本”的天然工具。当我在测试一个聚类算法时如果样本是两个重叠的高斯簇恰好能检验算法能否正确把混合区域里的样本分开。生成这样的合成数据非常简单前面定义的sample_two_gauss直接就能用而且可以附加标签把来自分量的 0 的样本标为类别 0来自分量的 1 的样本标为类别 1。这种做法在论文实验和模型评估中很实用因为真实标注数据往往稀少而合成数据可以控制簇间距、权重、噪声水平从而系统测试算法的鲁棒性。6.3 后续扩展多维、多峰与核密度对比这篇文章只展示了一维双峰的情况但现实中数据往往更复杂。如果你想继续扩展有三个方向很自然把两个高斯分量替换成多维高斯采样时用rng.multivariate_normalPDF 用scipy.stats.multivariate_normal.pdf把双峰扩展到三峰甚至更多混合权重变成一个概率向量仍然沿用“先按权重选分量再从分量抽样”的思路用核密度估计对比理论 PDF检查模型假设和真实样本的差异。多维情况的蒙特卡洛模拟码会复杂一些但核心理念完全一致先确定模型结构再设计抽样步骤最后用理论和经验分布互相验证。这也是我为什么一直强调要把一维向基础流程彻底搞明白的原因后面的扩展都是在这套流程上做加法。最后分享一个我的习惯每改动一次参数就把参数字典、随机种子、样本量记在代码注释里。双峰分布看起来只是多了三个参数但一旦开始调权重、调均值、调标准差组合爆炸的速度非常快不做好实验记录很容易忘了当前图对应的参数版本。把这些基础工作做扎实整个模拟流程才能真正成为后面分析工作的可靠支撑。
