做数据分析和统计模拟的兄弟应该都遇到过直方图画出两个峰的尴尬单高斯拟合怎么看都不对硬塞一个均值又两头不讨好。这种数据背后往往是两类群体混在一起对应的模型就是双峰高斯分布。要摸清它的脾气比如算某个区间的概率、找90%分位点手推公式太折磨人我一般直接上蒙特卡洛模拟——抽一两万个样本把PDF和CDF顺着画出来要什么有什么。这篇文章就把整套流程摊开讲一遍从分布参数设定、随机采样实现到两种绘图方式和踩坑经验适合刚接触统计模拟或者想做分布可视化但还没下过手的同学参考。之所以聊这个话题是因为最近帮朋友收拾一个数据分析脚本问题就出在用户使用时长分布上。数据画出来两个峰清清楚楚他想用单峰正态去近似结果拟合优度惨不忍睹。我帮他改成双峰高斯混合模型然后加了一段蒙特卡洛模拟半小时就把PDF和CDF出图了业务方一眼就看懂了。这个经验值得沉淀下来。1. 双峰高斯分布先搞懂模型再动手模拟1.1 混合模型公式与参数含义双峰高斯分布本质上不是一个新的分布而是两个正态分布按权重做混合标准写法是f(x) p * N(x; μ1, σ1) (1 - p) * N(x; μ2, σ2)其中p是第一个成分的混合权重取值范围在0到1之间(1-p)自然就是第二个成分的权重。这里的N就是普通正态分布的概率密度函数各自带着自己的均值和标准差。关键点在于这是加权平均而不是简单相加。把两个密度函数直接相加曲线下面积变成2积分不再是1就不符合概率密度函数的定义了。加权平均才能保证归一化这一点新手特别容易忽略。我这次项目里取的参数是μ1-2、σ11、μ23、σ21.3、p0.4。两个峰的中心相距5个标准差左右图型上可以清晰看到两个明显分离的峰。p0.4意味着40%的样本来自第一个成分60%来自第二个成分当样本量足够大时直方图的两个峰高大致也按这个比例分布。如果p接近0.5两个峰等高看起来更像对称的双峰p偏到0.8以上第一个峰就会明显高出一截。1.2 双峰出现在真实数据里的典型场景一个特别典型的场景是APP用户每日使用时长。轻度用户一天只打开几次累计使用时间集中在10分钟以内重度用户可能刷剧、刷短视频几个小时挂在上面。这两拨人各自内部的波动通常接近正态分布但把所有用户的数据放在一起直方图就会呈现两个峰。如果只看平均值大概率落在某个没多少人待着的中间位置对产品运营来说这个数字没有指导意义。要回答有多少人属于轻度用户、重度用户的典型使用时长是多少必须用混合模型分开建模。另一个例子是道路车速观测。白天通勤高峰以私家车为主车速分布有一批集中在60-80km/h深夜以货车为主车速集中在40-55km/h。混合后的车速数据就是典型的双峰结构。做交通仿真或者红绿灯配时方案时用单峰正态去拟合会严重低估昼夜两类车流的差异预测结果自然失真。这两个例子说明双峰高斯分布不是一个数学玩具它是对观测数据由两类群体混合而成这一事实的简洁建模方式。理解了模型结构后面做蒙特卡洛模拟才有根基。2. 蒙特卡洛模拟思路为什么采样比解析推导更实用2.1 解析方法难在哪如果只是画出PDF曲线公式摆在那边直接代值计算就行根本不需要蒙特卡洛。但真实工作中经常要回答另外两类问题一是求分位数比如90%的用户时长小于多少二是求区间概率比如车速落在40到60之间的概率有多大。混合分布的CDF确实可以写出来就是两个正态CDF的加权平均但它要求逆、求分位数时没有闭式解只能用数值迭代或者查表。区间概率虽然可以直接积分算但每换一个条件、每换一个区间就要重新推导一遍效率很低。更麻烦的是如果你把模型从双峰高斯换成三个成分、四个成分甚至换成t分布和正态分布的混合解析推导的工作量会迅速失控。蒙特卡洛模拟的思路很简单既然解析算不出来那就直接从这个分布里抽一大批样本然后用样本的经验分布去逼近真实分布。大数定律保证样本量足够大时经验CDF逐点收敛到理论CDF样本分位数也会收敛到理论分位数。布丰投针实验就是几百年前的雏形——朝格子里随机扔针来估算圆周率本质上就是用随机试验代替解析推导。现在计算机一秒钟生成百万级的随机数这种方法成本极其便宜。2.2 条件采样的实现逻辑从混合分布里抽样本标准做法是条件采样分两步走。第一步决定当前样本属于哪个成分生成一个在(0,1)区间均匀分布的随机数u如果u小于p就归到成分1否则归到成分2。第二步根据上一步选中的成分从对应的正态分布里抽一个值。数学依据是全概率分解混合分布的密度等于权重乘以各成分密度之和所以先按权重随机选成分、再按该成分分布采样最终样本的分布就自动满足混合分布。这样一段逻辑如果一开始没有想清楚很容易写出一个for循环逐个判断效率很低。用numpy的布尔掩码一次就能把样本分成两组然后分别用正态随机数发生器填充。样本量到十万级时效率差距还不明显到了百万级别for循环就会慢得让人怀疑人生。2.3 样本量怎么选才够稳样本量直接决定模拟结果的质量这也是很多人第一次跑蒙特卡洛最容易忽略的点。样本太少直方图都是毛刺双峰只能隐约看到两个鼓包KDE估计也不稳定样本太多计算时间上去了收益却边际递减。误差随样本量n按根号反比缩小从100到10000误差缩小为原来的十分之一但从10000到100000只是再缩小三倍。所以我自己画轮廓图时一般取10000到50000性价比最高。如果还需要做稳定的分位数估计有一个更实用的检验方法换几个不同的随机种子跑同一套参数看关键分位数的波动范围。比如95%分位数在多次运行之间波动小于0.1说明样本量够用如果每次都差出半个单位就该加样本量了。这比空谈理论更直接也更适合项目实战。3. 完整实现代码采样、PDF与CDF绘图3.1 环境准备与依赖安装实现整个模拟只需要三个库numpy负责随机数生成和数组操作scipy提供正态分布的概率密度函数和累积分布函数matplotlib负责绘图。Python 3.8以上的环境基本都能跑如果机器上还没有这些库装一下也很简单pip install numpy scipy matplotlib有Anaconda环境的同学直接用conda install也一样。这三个库是数据科学领域最常见的基础组合版本选择上没有特殊要求用最新的稳定版就行。3.2 采样函数实现我习惯把采样逻辑封装成一个函数方便换参数反复调用。参数放在函数签名里每次模拟只要改几行参数就能跑出新结果。这里用default_rng而不是老的np.random.seed加np.random.rand的组合主要是因为default_rng是numpy当前推荐的新接口每次生成独立的随机数流不会污染全局随机状态。固定种子之后两次跑出来的样本完全一致这对调试和复现非常重要。import numpy as np from scipy.stats import norm, gaussian_kde import matplotlib.pyplot as plt def sample_mixture(mu1, sigma1, mu2, sigma2, p, n20000, seed42): rng np.random.default_rng(seed) u rng.random(n) mask_first u p samples np.empty(n) samples[mask_first] rng.normal(mu1, sigma1, sizemask_first.sum()) samples[~mask_first] rng.normal(mu2, sigma2, size(~mask_first).sum()) return samples mu1, sigma1 -2.0, 1.0 mu2, sigma2 3.0, 1.3 p 0.4 n 20000 samples sample_mixture(mu1, sigma1, mu2, sigma2, p, nn, seed42)mask_first.sum()在n足够大时约等于pn但由于随机性会有微小波动这不是错误而是随机抽样的正常表现。如果你强制让第一个成分的样本数精确等于pn反而破坏了随机性得到的样本分布会有微小偏差。3.3 PDF绘图直方图、理论曲线与KDEPDF的全称是概率密度函数。理论PDF可以直接按加权公式计算在横轴上生成一批密集的x值分别计算两个正态分布在每个x处的密度再按权重加起来。同时把采样的样本画成直方图再补一条核密度估计KDE曲线。KDE的好处在于只从样本出发就能平滑地还原密度形态不需要事先知道真实分布这在探索性数据分析里非常实用。x np.linspace(-8, 9, 1000) pdf_theory (p * norm.pdf(x, mu1, sigma1) (1 - p) * norm.pdf(x, mu2, sigma2)) kde gaussian_kde(samples, bw_method0.3) pdf_kde kde(x) plt.figure(figsize(7, 5)) plt.hist(samples, bins80, densityTrue, alpha0.35, colorskyblue, label样本直方图) plt.plot(x, pdf_theory, r-, lw2, label理论PDF) plt.plot(x, pdf_kde, g--, lw2, labelKDE估计) plt.xlabel(x) plt.ylabel(概率密度) plt.legend() plt.show()直方图里有个细节必须注意densityTrue一定要写否则纵轴是频数曲线是概率密度两个量纲根本对不上。另外bins的数量也会影响直方图的表现80个bins在这个数据量下足够看到清晰的峰型太少会把双峰糊成一个宽包太多又会出现间隙和毛刺。从KDE曲线可以直观看到带宽的作用。bw_method0.3时曲线基本贴合理论PDF两个峰都清楚。带宽如果调到0.8双峰被磨平只剩一个鼓包调到0.1曲线会剧烈波动出现很多假的抖动细节。带宽的选择没有绝对标准和两个成分的标准差有关后面会专门展开讲。3.4 CDF绘图经验分布与理论分布对照CDF是累积分布函数表示随机变量取值小于等于某个x的概率。理论CDF和理论PDF一样用加权公式直接计算。经验CDF则是从样本出发把样本排序后第i个点的纵坐标就是i/n。这条线用样本直接画出不需要任何分布假设。cdf_theory (p * norm.cdf(x, mu1, sigma1) (1 - p) * norm.cdf(x, mu2, sigma2)) sorted_samples np.sort(samples) ecdf_y np.arange(1, len(sorted_samples) 1) / len(sorted_samples) plt.figure(figsize(7, 5)) plt.plot(sorted_samples, ecdf_y, b-, lw1.5, label经验CDF) plt.plot(x, cdf_theory, r--, lw2, label理论CDF) plt.xlabel(x) plt.ylabel(累积概率) plt.legend() plt.grid(alpha0.3) plt.show()两条线基本重合这是蒙特卡洛有效性的直观证据。如果你发现两条线明显错开问题大概率出在采样参数或者权重上不可能是理论的锅。得到CDF之后分位数和区间概率就都能直接读了。90%分位数就是经验CDF纵坐标等于0.9时对应的横坐标值在数组里直接取索引q90 sorted_samples[int(n * 0.90) - 1] print(f90% 分位数: {q90:.3f}) p_interval ((samples -1) (samples 2)).mean() print(fP(-1 X 2) 模拟值: {p_interval:.3f})上面这行((samples -1) (samples 2)).mean()非常简洁numpy的布尔数组求均值就是True的比例因而一步就得到区间概率的模拟值。这段代码就是蒙特卡洛的核心魅力换一个问题只改一行判断条件不需要重新推公式。4. 绘图细节处理布局、中文显示与文件导出4.1 图面排版与样式调整如果要把PDF和CDF放在同一张图里对比展示用subplots画左右两个子图更合适。子图并排的好处是两张图共享同一套横轴逻辑视觉上便于对分布形态和累积概率进行关联理解。我会把图幅设成12x5左右子图各占一半横轴范围保持一致的linspace。PDF子图重点展示形态CDF子图重点展示累积概率和分位数两者成对看效果最好。样式上的几个细节也值得注意x轴标签写清楚变量名和单位PDF的纵轴是概率密度CDF的纵轴是累积概率这两个词不能混。图例放在不遮挡曲线的位置通常右上角或者右下角。网格线透明度设到0.3既能看到刻度参考又不会喧宾夺主。线条宽度PDF曲线设2CDF经验线设1.5层次上让理论曲线更突出。4.2 中文显示与PDF矢量导出matplotlib默认字体不支持中文一出现中文标签就是满屏方块。这个问题在Linux服务器上尤其常见。我用的配置是这样的plt.rcParams[font.sans-serif] [SimHei, Microsoft YaHei, Noto Sans CJK SC] plt.rcParams[axes.unicode_minus] False第一行把黑体、微软雅黑、Noto字体都放进候选列表系统里装了哪个就用哪个Windows一般命中SimHeiLinux一般命中Noto Sans CJK SC。第二行专门解决负号显示成方块的问题因为设置了中文字体后默认的ASCII负号可能被替换成不支持的字符必须显式关闭unicode_minus。出图保存也有讲究。日常预览用PNGdpi设200就够了体积小、加载快。如果要放到论文或者交付文档里应该保存PDF矢量图放大不模糊。一行代码搞定fig.savefig(mixture_pdf_cdf.pdf, bbox_inchestight) fig.savefig(mixture_pdf_cdf.png, dpi200, bbox_inchestight)bbox_inchestight会自动裁掉多余留白只保留绘图区内容这让图片在文档里插入时不需要再手动裁剪。5. 常见问题与避坑实录5.1 分布形状对不上理论曲线怎么办先说一个我踩过的坑权重参数p用反。p0.4代表第一个成分为0.4、第二个成分为0.6。如果你在采样时把p当成第二个成分的权重抽样分布就和理论PDF完全错位。检查方法很简单画完图看第一个峰和第二个峰的相对高度再用样本均值做个粗略验证。理论均值是pμ1 (1-p)μ2代入这次参数是0.4(-2)0.631.0样本均值在1.0附近说明权重方向没问题。如果两个峰糊成了一个宽峰先检查μ1和μ2的间距以及σ1、σ2的取值。两个成分的中心如果只差1个标准差混合出来的峰型大概率只有一个宽包这是模型本身的特性不是画图错误。大致的经验是两个中心间距大于2倍标准差之和时双峰才肉眼可见。想要更早看到峰谷可以做峰度检验或bimodality系数比肉眼更客观。5.2 KDE带宽到底调多少KDE的带宽bw_method是这条曲线最关键的参数它控制平滑程度。带宽太小KDE会跟着每个样本抖动出现许多虚假的小波动带宽太大真实的双峰结构会被磨平。scipy的gaussian_kde默认会根据数据自动选带宽但自动选出来的值在这个场景下往往偏大双峰会显得不够锐利。我实测下来的经验是先用默认值看整体形态再逐步调小到0.3左右对比。这次项目里两个成分的标准差分别是1和1.3用bw_method0.3效果好两个峰和理论曲线贴合度很高。如果两个成分的标准差都是0.5左右的尖峰带宽0.2会更合适。没有万能参数多试几个值选一个在平滑和保真之间平衡的就好。5.3 直方图和PDF的y轴对不上这是新手最容易疑惑的点为什么PDF曲线高度和直方图柱子高度差距这么大是不是画错了大概率不是是直方图的density参数没设置。直方图默认纵轴是频数也就是每个bin里的样本计数而PDF的纵轴是概率密度表示单位区间上的概率分布强度。两者的量纲完全不同叠在一起当然对不上。设置了densityTrue之后直方图的纵轴就变成密度和PDF曲线同量纲。这时要注意概率密度在一个点上的取值可以大于1这不违反概率规则因为单点的概率是0只有在区间上积分才有概率意义。我在带新人的时候经常说一句PDF的高度不是概率面积才是概率。5.4 复现问题怎么解决模拟结果第二天复现不出来是很多团队的痛。这里有个三件套固定随机种子、锁定依赖版本、记录关键参数。default_rng的固定种子在同一个numpy版本下结果是稳定的但如果numpy跨大版本升级随机数算法有过变动结果可能会不一样。我的习惯是在项目仓库里写一个requirements.txt把numpy、scipy、matplotlib的主版本都钉住。另外如果你在Notebook里跑代码中途执行过其他随机数操作再用同一个种子跑采样结果也可能不同。default_rng比全局随机状态好在每次实例化都是独立的随机数流只要你每次从同一个种子实例化就不会被外部随机操作污染。建议明确一条规矩要做复现实验就从创建rng开始重跑整个cell不要依赖断点续跑。5.5 工具横向对比和选型建议最后聊一下工具选择。matplotlib是数据探索阶段最顺手的代码写完立刻出图适合快速迭代seaborn对KDE和分布图封装得更简洁一行kdeplot就能出密度曲线R语言里的ggplot2做统计图形在排版上更讲究很多论文级配图用它。Origin则是交互式操作的代表适合不喜欢写代码、想靠鼠标完成精修的人但批量处理多个参数组合时还是脚本工具更省事。至于qt绘图、canvas、echarts这些是桌面软件和Web前端的绘图方案面向的是把控件嵌到应用里和网页上和大屏上做交互展示这样的场景跟统计模拟不是同一个赛道。如果你只是在分析一个分布形态花半小时搭一套Web可视化得不偿失。先打通模拟逻辑再考虑炫酷的展示这是我在项目里的优先级排序。做这类模拟我个人的体会是别急着一次跑出完美图先把采样函数、理论公式、画图代码拆成独立模块每个模块单独验证。采样出的样本均值要和理论均值对得上理论PDF曲线下的面积积分要接近1经验CDF的终点要落在1附近。这些小验证加起来只需要几行代码却能筛掉绝大部分低级错误。我早期有一次就把p写反了图怎么改都不对最后还是用均值这个最简单的检查暴露了问题。把这些检查养成习惯蒙特卡洛模拟基本就是一条顺畅的流水线。