概率论这门课我当年学的时候最大的感受就是公式推得挺顺一到用代码落地就卡壳。尤其是概型和概率分布这一块教材上给的是定义、性质、期望方差考试考的是计算但真到了要写代码做仿真、做数据分析、做风险评估的时候很多人连这个场景该用哪个分布都判断不出来。这篇内容就是冲着这个痛点来的——把常见的概型古典概型、几何概型、条件概型和几大核心概率分布二项、泊松、几何、均匀、指数、正态用Python从头实现一遍不光给代码还要讲清楚每个分布背后的直觉、适用边界以及我在实际写代码时踩过的那些坑。适合谁看如果你正在学概率统计想把书本知识变成能跑的代码这篇对你有用如果你做数据分析、机器学习经常需要判断数据服从什么分布、怎么采样、怎么验证这篇也能帮你把基础打牢。我不打算写成教科书而是按这个分布解决什么问题—数学长什么样—代码怎么写—实测注意什么的顺序来讲尽量让每个概念都能落到键盘上。1. 先把概型这件事用代码说清楚很多人一上来就急着写分布函数其实概型是更底层的东西——它决定了你怎么数概率。古典概型靠数样本点几何概型靠算面积/长度条件概型靠缩小样本空间。这三者在代码里的实现思路完全不同混着用很容易出错。1.1 古典概型等可能样本点的暴力枚举与组合计数古典概型的核心假设是样本空间有限且每个基本事件等可能。掷骰子、抽扑克、抛硬币都是典型。代码实现有两条路——暴力枚举和组合公式。暴力枚举直观但只适合小规模组合计数效率高但容易在有序/无序放回/不放回上翻车。先看暴力枚举。以掷两枚骰子点数和为7的概率为例import itertools outcomes list(itertools.product(range(1, 7), repeat2)) target [o for o in outcomes if sum(o) 7] print(len(target) / len(outcomes)) # 6/36 0.1667这段代码把36种结果全列出来数出和为7的6种。逻辑没毛病但样本空间一大就爆内存。比如掷10枚骰子样本空间是6的10次方约6000万枚举就不现实了。这时候得上组合计数。用math.comb算组合数from math import comb # 从52张牌里抽5张恰好2张A的概率 total comb(52, 5) favorable comb(4, 2) * comb(48, 3) print(favorable / total) # 约0.0399这里有个我踩过的坑comb(4,2)是从4张A里选2张comb(48,3)是从剩下48张非A里选3张两者相乘才是恰好2张A的有利事件数。很多人会写成comb(4,2)*comb(50,3)忘了抽走2张A后牌堆只剩50张这是典型的样本空间没同步缩小。提示用组合计数时先问自己三个问题——是否有序是否放回样本空间在每一步是否变化这三个问题的答案决定了你用的是排列还是组合、分母要不要跟着变。古典概型的代码实现还有个隐蔽问题浮点精度。当样本空间极大时favorable/total可能因为整数除法丢精度。稳妥做法是用fractions.Fraction保留精确分数最后再转浮点from fractions import Fraction p Fraction(favorable, total) print(float(p))1.2 几何概型用蒙特卡洛把面积比变成可执行代码几何概型的概率等于有利区域测度/总区域测度测度可以是长度、面积、体积。经典例子是在单位正方形内随机取点落在内切圆内的概率是π/4。这类问题解析解往往好求但一旦区域形状不规则就得靠蒙特卡洛仿真。蒙特卡洛的思路很朴素撒点数落在目标区域的比例。代码实现的关键是随机数生成器的选择和样本量控制import numpy as np def estimate_pi(n_samples1_000_000, seed42): rng np.random.default_rng(seed) points rng.random((n_samples, 2)) # [0,1) x [0,1) inside np.sum(np.sum((points - 0.5)**2, axis1) 0.25) return 4 * inside / n_samples print(estimate_pi())这里我特意用了np.random.default_rng(seed)而不是老的np.random.seed。前者是NumPy新版推荐的Generator接口随机数质量更好且支持独立的多路随机流做并行仿真时不会互相干扰。实测下来100万样本估计π的误差通常在0.001量级想要更高精度得加样本量但误差收敛是O(1/√n)样本量翻100倍误差才降10倍性价比不高。几何概型代码里最容易忽略的是均匀性假设。rng.random()生成的是[0,1)上的均匀分布如果你要在一个非矩形区域均匀采样直接对坐标做均匀采样可能不对。比如在圆内均匀取点不能对半径和角度分别均匀采样——那样点会往圆心聚集。正确做法是拒绝采样或对半径做平方根变换# 圆内均匀采样半径取 sqrt(U) * R rng np.random.default_rng(0) n 100000 r np.sqrt(rng.random(n)) * 1.0 theta rng.random(n) * 2 * np.pi x, y r * np.cos(theta), r * np.sin(theta)这个平方根变换的直觉是圆环面积随半径线性增长所以半径的采样密度要正比于r累积分布就是r²反变换采样自然得到sqrt(U)。1.3 条件概型贝叶斯更新在代码里的落地形态条件概型处理的是已知某事件发生求另一事件概率。公式P(A|B)P(AB)/P(B)谁都背得出但代码里真正有价值的是把它做成可复用的贝叶斯更新器。举个医学检测的例子某病患病率1%检测灵敏度99%特异度95%问检测阳性时真患病的概率。def bayes_update(prior, sensitivity, specificity): p_pos_given_disease sensitivity p_pos_given_healthy 1 - specificity p_pos prior * p_pos_given_disease (1 - prior) * p_pos_given_healthy return prior * p_pos_given_disease / p_pos print(bayes_update(0.01, 0.99, 0.95)) # 约0.1667结果只有16.67%很多人第一次看到会惊讶——检测这么准阳性了真患病概率还这么低原因在于先验患病率太低假阳性数量压过了真阳性。这个例子在代码里的价值是提醒你贝叶斯更新的结果对先验极度敏感写代码时先验参数一定要有依据不能拍脑袋。条件概型还有个常见实现场景是马尔可夫链的转移概率。给定当前状态下一状态的概率只依赖当前状态这本质就是条件概型。代码里用转移矩阵表示transition np.array([[0.9, 0.1], [0.3, 0.7]]) # 从状态0出发两步后到状态1的概率 print((transition transition)[0, 1])矩阵乘法的顺序别搞反transition transition才是两步转移transition.T transition是另一回事。我见过有人把转移矩阵转置了还在纳闷结果不对这种错误在调试时特别隐蔽。2. 离散分布从伯努利到泊松的代码谱系离散分布处理的是可数结果的随机变量。伯努利是最小单元二项是伯努利的累加几何是第一次成功前的等待泊松是单位时间内的稀有事件计数。这四个分布之间有清晰的派生关系代码实现时理解这层关系比死记公式有用得多。2.1 伯努利与二项一次和多次的递进关系伯努利分布就是一次成败试验参数p是成功概率。代码实现简单到不需要库def bernoulli(p, rng): return 1 if rng.random() p else 0二项分布是n次独立伯努利试验中成功的次数。NumPy直接有rng.binomial(n, p, size)。但我想强调的是自己实现一遍的价值——用循环累加伯努利能直观看到二项是怎么长出来的def binomial_manual(n, p, rng): return sum(bernoulli(p, rng) for _ in range(n)) rng np.random.default_rng(42) samples [binomial_manual(10, 0.3, rng) for _ in range(10000)] print(np.mean(samples), np.var(samples)) # 约3.0, 约2.1理论均值np3方差np(1-p)2.1实测吻合。这里有个性能坑Python循环跑1万次还行跑100万次就慢了。生产环境直接用rng.binomial底层是C实现快两个数量级。二项分布代码里最需要警惕的是独立性假设。如果试验之间不独立比如抽样不放回二项就不适用得换成超几何分布。我见过做A/B测试的人直接用二项算显著性但用户行为往往有聚集性独立性不成立结论就偏了。2.2 几何分布等待首次成功的两种参数化陷阱几何分布描述首次成功所需的试验次数。这里有个大坑不同教材对几何分布的定义不同。有的定义成首次成功前的失败次数取值0,1,2,...有的定义成首次成功时的试验次数取值1,2,3,...。两者差1均值和方差公式也不同。# 定义A首次成功时的试验次数取值1,2,... def geometric_trials(p, rng): trials 0 while True: trials 1 if rng.random() p: return trials # 定义B首次成功前的失败次数取值0,1,... def geometric_failures(p, rng): return geometric_trials(p, rng) - 1NumPy的rng.geometric(p)用的是定义A取值从1开始。如果你按定义B推导了公式又用NumPy采样结果会系统性偏1。这个坑我在做可靠性分析时踩过——算平均无故障时间公式用的是失败次数版本采样用的是NumPy最后结果差了整整一个单位时间排查了半天才发现是参数化不一致。几何分布还有个无记忆性P(Xmn|Xm)P(Xn)。代码验证很直观p 0.2 # 已经等了5次还没成功再等3次仍未成功的条件概率 # 应该等于从头等3次仍未成功的概率无记忆性是几何分布独有的连续情形里指数分布也有这个性质在建模等待时间时非常有用但也意味着它不适合描述越等越可能成功的场景——那种场景得用别的分布。2.3 泊松分布稀有事件计数的均值方差相等特性泊松分布描述单位时间/空间内稀有事件的发生次数参数λ是平均发生率。它的标志性特征是均值等于方差都等于λ。这个特征在代码里可以用来做快速校验rng np.random.default_rng(7) lam 4.0 samples rng.poisson(lam, size100000) print(samples.mean(), samples.var()) # 都接近4.0泊松和二项的关系值得说清楚当n很大、p很小时二项分布近似泊松λnp。代码验证n, p 10000, 0.0004 binom_samples rng.binomial(n, p, size100000) poisson_samples rng.poisson(n * p, size100000) # 两者分布形态几乎重合这个近似在实际中很有用——当n大到无法逐次模拟时用泊松替代二项能大幅提速。但要注意近似条件一般要求n≥20且p≤0.05n越大p越小近似越好。我见过有人n50、p0.3还用泊松近似误差肉眼可见。泊松分布代码实现时有个容易忽略的点事件必须独立且发生率恒定。如果事件有聚集性比如网站访问量在促销期间暴增泊松就不适用得考虑负二项分布或复合泊松。做流量预测时直接套泊松遇到突发流量就会严重低估尾部风险。3. 连续分布均匀、指数、正态的实现细节连续分布处理的是不可数结果的随机变量概率通过密度函数积分得到。代码实现的核心是采样和密度计算但每个分布都有各自的数值陷阱。3.1 均匀分布看似最简单却最容易用错边界均匀分布U(a,b)的密度在[a,b]上是1/(b-a)外面是0。NumPy的rng.uniform(low, high, size)生成的是[low, high)上的均匀分布——注意是左闭右开。这个细节在离散化时会影响结果rng np.random.default_rng(0) # 生成[0,1)均匀数转成0-9的整数 digits np.floor(rng.random(100000) * 10).astype(int) # 各数字出现频率应接近0.1用np.floor而不是np.round因为rng.random()取不到1乘以10后范围是[0,10)floor后正好是0-9。如果用round0和9的出现概率会只有其他数字的一半因为边界只覆盖了半个区间。这个坑在做随机分组、随机抽样时特别常见。均匀分布还有个用途是逆变换采样——任何分布的采样都可以通过均匀分布累积分布反函数实现。以指数分布为例def exponential_inverse_transform(lam, size, rng): u rng.random(size) return -np.log(1 - u) / lam这就是逆变换采样的标准形态。理解这一点你就能自己实现任意分布的采样器不必依赖现成库。3.2 指数分布与泊松的孪生关系及无记忆性验证指数分布描述泊松过程中相邻事件的时间间隔参数λ是事件发生率。它和泊松是计数与间隔的孪生关系泊松数事件个数指数测事件间隔。rng np.random.default_rng(1) lam 2.0 intervals rng.exponential(1/lam, size100000) print(intervals.mean()) # 约0.5即1/lam注意NumPy的rng.exponential(scale)参数是scale1/λ不是λ本身。这个参数化方向和很多教材相反我第一次用时把λ直接传进去结果均值差了4倍。记忆方法scale是尺度单位和时间一样λ是速率单位是时间的倒数。指数分布的无记忆性可以用代码直观验证# P(X 23 | X 2) 应等于 P(X 3) lam 1.0 # 理论值 p_conditional np.exp(-lam * 3) # 仿真验证 samples rng.exponential(1/lam, size1_000_000) mask samples 2 print(np.mean(samples[mask] 5)) # 接近 exp(-3)无记忆性意味着已经等了多久不影响还要等多久这在建模设备寿命、客户到达间隔时是强假设。如果实际过程有老化或疲劳效应指数分布就不合适得用威布尔分布。3.3 正态分布中心极限定理的代码实证与参数陷阱正态分布是连续分布里用得最多的参数μ和σ分别控制位置和尺度。NumPy的rng.normal(loc, scale, size)里scale就是标准差σ不是方差。这个区分在写代码时经常被混淆rng np.random.default_rng(2) samples rng.normal(loc10, scale2, size100000) print(samples.mean(), samples.std()) # 约10, 约2正态分布最值得用代码验证的是中心极限定理大量独立同分布随机变量之和近似正态。拿均匀分布做实验rng np.random.default_rng(3) n, m 30, 10000 sums rng.random((m, n)).sum(axis1) # sums应近似正态均值n/215方差n/122.5n30时已经相当接近正态了。这个实验的价值在于让你对什么时候可以用正态近似有手感——n小于10时偏态还很明显n到30左右基本可用这也是很多统计方法默认样本量30的由来。正态分布代码里有个隐蔽陷阱用样本标准差估计σ时的自由度。np.std默认ddof0除以n而统计推断里通常要ddof1除以n-1才是无偏估计。做假设检验时用错ddof小样本下结论会偏data rng.normal(0, 1, size10) print(np.std(data, ddof0), np.std(data, ddof1))样本量小时两者差异明显样本量大时趋于一致。这个细节在写统计检验代码时务必注意。4. 分布选择与验证从数据反推分布的实操路径前面讲了各分布的实现但实际工作中更常见的需求是反过来的手上有一堆数据该用什么分布建模这一节讲我的实操路径——先看数据形态再做拟合检验最后用仿真验证。4.1 用直方图和矩特征做初步判断拿到数据第一步不是套分布而是看形态。直方图看偏度和峰度均值方差看离散程度。我通常写一个快速诊断函数import numpy as np from scipy import stats def diagnose(data): print(f均值: {np.mean(data):.3f}) print(f方差: {np.var(data):.3f}) print(f偏度: {stats.skew(data):.3f}) print(f峰度: {stats.kurtosis(data):.3f})判断逻辑均值≈方差且取值为非负整数优先考虑泊松偏度明显为正且取值连续非负考虑指数或对数正态偏度接近0、峰度接近0考虑正态。这个初判不能替代正式检验但能快速排除明显不合适的候选。我踩过的一个坑是只看均值和方差就下结论。有次数据均值方差都是5左右我以为是泊松结果直方图一看是双峰——明显是两个泊松的混合。所以直方图这一步不能省矩特征只能辅助。4.2 KS检验与拟合优度的代码模板初步锁定候选分布后用Kolmogorov-Smirnov检验做拟合优度。SciPy的stats.kstest可以直接用from scipy import stats data rng.poisson(4, size1000) # 检验是否服从泊松分布 stat, p_value stats.kstest(data, poisson, args(4,)) print(fKS统计量: {stat:.4f}, p值: {p_value:.4f})p值大于0.05表示不能拒绝数据服从该分布的原假设。但这里有个重要提醒KS检验对参数估计后的分布会偏保守因为参数是从数据估出来的自由度损失了。严谨做法是用Lilliefors检验或蒙特卡洛校正。实际工作中如果p值在0.05附近徘徊别急着下结论多换几个样本看看稳定性。还有个常见误区是把不能拒绝当成服从。KS检验不拒绝只说明数据与假设分布没有显著矛盾不代表就是那个分布。样本量越大检验越敏感微小偏差都会被检出。所以大样本下p值小不一定意味着分布选错了可能只是偏差在统计上显著但实际影响很小。4.3 采样一致性验证把拟合好的分布跑一遍仿真选定分布并估计参数后我习惯做一步回测用拟合的分布重新采样对比仿真数据和原始数据的关键统计量。这一步能暴露参数估计的偏差def backtest(data, dist_name, params, n_sim10000): sim getattr(rng, dist_name)(*params, sizen_sim) print(f原始均值: {np.mean(data):.3f}, 仿真均值: {np.mean(sim):.3f}) print(f原始方差: {np.var(data):.3f}, 仿真方差: {np.var(sim):.3f}) # 对比分位数 for q in [0.05, 0.5, 0.95]: print(fq{q}: 原始{np.quantile(data, q):.3f}, 仿真{np.quantile(sim, q):.3f})分位数对比比均值方差更能反映尾部行为。如果原始数据在95%分位数上明显高于仿真说明拟合分布低估了尾部风险——这在金融、可靠性领域是致命问题得换重尾分布如t分布、帕累托分布。回测还有个作用是发现参数估计方法的问题。比如用矩估计拟合正态遇到离群点时σ会被高估用最大似然估计则相对稳健。代码里scipy.stats.norm.fit默认用MLE但你可以手动指定mu_mle, sigma_mle stats.norm.fit(data, methodMM) # 矩估计不同估计方法在小样本或含噪数据下结果差异明显回测能帮你判断哪种更合适。5. 几个我反复踩过的实现坑与规避手法这一节单独拎出来讲因为下面这些问题不是某个分布特有的而是跨分布通用的且往往在调试时最难定位。5.1 随机种子管理与可复现性做概率仿真可复现性是底线。NumPy新版推荐用default_rng(seed)创建独立的Generator实例而不是全局的np.random.seed。区别在于全局种子会影响所有后续随机调用多个函数共享时互相干扰独立Generator各管各的互不影响。# 推荐独立Generator rng1 np.random.default_rng(42) rng2 np.random.default_rng(42) # rng1和rng2产生相同序列但互不干扰 # 不推荐全局种子 np.random.seed(42)做并行仿真时给每个进程分配不同的种子偏移避免所有进程产生相同序列from concurrent.futures import ProcessPoolExecutor def simulate(seed): rng np.random.default_rng(seed) return rng.normal(size1000).mean() with ProcessPoolExecutor() as ex: results list(ex.map(simulate, range(8)))如果所有进程都用同一个种子8个结果会完全一样仿真就白做了。这个坑我在第一次写并行蒙特卡洛时踩得结结实实。5.2 浮点精度与极端概率的处理概率计算里浮点精度问题在尾部概率上尤其突出。比如计算正态分布在极端值处的概率直接算会下溢到0from scipy import stats # 直接算 print(stats.norm.pdf(40)) # 下溢为0 # 用对数密度 print(stats.norm.logpdf(40)) # 有限值做似然计算、贝叶斯推断时一律用对数域运算最后再取指数。这个习惯能避免大量数值问题。另一个精度坑是累积概率的相减。P(aXb)F(b)-F(a)当a和b很接近时两个接近1的数相减会丢有效数字。稳妥做法是用stats.norm.cdf(b) - stats.norm.cdf(a)在b-a较小时改用stats.norm.sf(a) - stats.norm.sf(b)sf是生存函数算的是右尾或者直接用区间概率函数。5.3 分布参数化的方向混淆前面提过指数分布的scale和rate其实还有几个分布有类似问题分布NumPy参数常见教材参数易错点指数scale1/λλ速率传错方向差λ²倍正态scaleσσ标准差误传方差伽马scaleθ, shapek有时用rate1/θ方向反了贝塔a, b有时用均值和集中度参数含义不同规避方法写代码时在变量名里带上单位或含义比如rate_lambda、scale_theta别用裸的p、a、b。这个习惯看着啰嗦但能省下大量调试时间。5.4 大样本下的性能与内存做百万级采样时内存和速度都要考虑。rng.normal(size10_000_000)一次生成1000万个浮点数占80MB内存还行。但如果要生成多个分布做对比累积起来就吃紧了。技巧是分块生成、流式处理def streaming_mean(dist_func, total, chunk100000): acc 0.0 for _ in range(total // chunk): acc dist_func(chunk).sum() return acc / total这样内存占用恒定在chunk级别适合超大规模仿真。另外能用向量化就别用循环NumPy的向量化操作比Python循环快一到两个数量级这个差距在大样本下非常明显。概率分布的代码实现说到底就是把数学定义翻译成可执行、可验证、可复现的计算过程。我个人的体会是每学一个分布别急着背公式先用代码把它采样出来、画出直方图、算一遍均值和方差和理论值对一对。这个采样—验证的闭环走多了你对每个分布的直觉会越来越准遇到新场景时判断该用哪个分布也会越来越快。至于那些参数化方向、边界条件、精度陷阱都是在一次次跑代码、对结果的过程中慢慢刻进肌肉记忆的——光看是记不住的得亲手跑一遍。
概率论这门课我当年学的时候最大的感受就是公式推得挺顺一到用代码落地就卡壳。尤其是概型和概率分布这一块教材上给的是定义、性质、期望方差考试考的是计算但真到了要写代码做仿真、做数据分析、做风险评估的时候很多人连这个场景该用哪个分布都判断不出来。这篇内容就是冲着这个痛点来的——把常见的概型古典概型、几何概型、条件概型和几大核心概率分布二项、泊松、几何、均匀、指数、正态用Python从头实现一遍不光给代码还要讲清楚每个分布背后的直觉、适用边界以及我在实际写代码时踩过的那些坑。适合谁看如果你正在学概率统计想把书本知识变成能跑的代码这篇对你有用如果你做数据分析、机器学习经常需要判断数据服从什么分布、怎么采样、怎么验证这篇也能帮你把基础打牢。我不打算写成教科书而是按这个分布解决什么问题—数学长什么样—代码怎么写—实测注意什么的顺序来讲尽量让每个概念都能落到键盘上。1. 先把概型这件事用代码说清楚很多人一上来就急着写分布函数其实概型是更底层的东西——它决定了你怎么数概率。古典概型靠数样本点几何概型靠算面积/长度条件概型靠缩小样本空间。这三者在代码里的实现思路完全不同混着用很容易出错。1.1 古典概型等可能样本点的暴力枚举与组合计数古典概型的核心假设是样本空间有限且每个基本事件等可能。掷骰子、抽扑克、抛硬币都是典型。代码实现有两条路——暴力枚举和组合公式。暴力枚举直观但只适合小规模组合计数效率高但容易在有序/无序放回/不放回上翻车。先看暴力枚举。以掷两枚骰子点数和为7的概率为例import itertools outcomes list(itertools.product(range(1, 7), repeat2)) target [o for o in outcomes if sum(o) 7] print(len(target) / len(outcomes)) # 6/36 0.1667这段代码把36种结果全列出来数出和为7的6种。逻辑没毛病但样本空间一大就爆内存。比如掷10枚骰子样本空间是6的10次方约6000万枚举就不现实了。这时候得上组合计数。用math.comb算组合数from math import comb # 从52张牌里抽5张恰好2张A的概率 total comb(52, 5) favorable comb(4, 2) * comb(48, 3) print(favorable / total) # 约0.0399这里有个我踩过的坑comb(4,2)是从4张A里选2张comb(48,3)是从剩下48张非A里选3张两者相乘才是恰好2张A的有利事件数。很多人会写成comb(4,2)*comb(50,3)忘了抽走2张A后牌堆只剩50张这是典型的样本空间没同步缩小。提示用组合计数时先问自己三个问题——是否有序是否放回样本空间在每一步是否变化这三个问题的答案决定了你用的是排列还是组合、分母要不要跟着变。古典概型的代码实现还有个隐蔽问题浮点精度。当样本空间极大时favorable/total可能因为整数除法丢精度。稳妥做法是用fractions.Fraction保留精确分数最后再转浮点from fractions import Fraction p Fraction(favorable, total) print(float(p))1.2 几何概型用蒙特卡洛把面积比变成可执行代码几何概型的概率等于有利区域测度/总区域测度测度可以是长度、面积、体积。经典例子是在单位正方形内随机取点落在内切圆内的概率是π/4。这类问题解析解往往好求但一旦区域形状不规则就得靠蒙特卡洛仿真。蒙特卡洛的思路很朴素撒点数落在目标区域的比例。代码实现的关键是随机数生成器的选择和样本量控制import numpy as np def estimate_pi(n_samples1_000_000, seed42): rng np.random.default_rng(seed) points rng.random((n_samples, 2)) # [0,1) x [0,1) inside np.sum(np.sum((points - 0.5)**2, axis1) 0.25) return 4 * inside / n_samples print(estimate_pi())这里我特意用了np.random.default_rng(seed)而不是老的np.random.seed。前者是NumPy新版推荐的Generator接口随机数质量更好且支持独立的多路随机流做并行仿真时不会互相干扰。实测下来100万样本估计π的误差通常在0.001量级想要更高精度得加样本量但误差收敛是O(1/√n)样本量翻100倍误差才降10倍性价比不高。几何概型代码里最容易忽略的是均匀性假设。rng.random()生成的是[0,1)上的均匀分布如果你要在一个非矩形区域均匀采样直接对坐标做均匀采样可能不对。比如在圆内均匀取点不能对半径和角度分别均匀采样——那样点会往圆心聚集。正确做法是拒绝采样或对半径做平方根变换# 圆内均匀采样半径取 sqrt(U) * R rng np.random.default_rng(0) n 100000 r np.sqrt(rng.random(n)) * 1.0 theta rng.random(n) * 2 * np.pi x, y r * np.cos(theta), r * np.sin(theta)这个平方根变换的直觉是圆环面积随半径线性增长所以半径的采样密度要正比于r累积分布就是r²反变换采样自然得到sqrt(U)。1.3 条件概型贝叶斯更新在代码里的落地形态条件概型处理的是已知某事件发生求另一事件概率。公式P(A|B)P(AB)/P(B)谁都背得出但代码里真正有价值的是把它做成可复用的贝叶斯更新器。举个医学检测的例子某病患病率1%检测灵敏度99%特异度95%问检测阳性时真患病的概率。def bayes_update(prior, sensitivity, specificity): p_pos_given_disease sensitivity p_pos_given_healthy 1 - specificity p_pos prior * p_pos_given_disease (1 - prior) * p_pos_given_healthy return prior * p_pos_given_disease / p_pos print(bayes_update(0.01, 0.99, 0.95)) # 约0.1667结果只有16.67%很多人第一次看到会惊讶——检测这么准阳性了真患病概率还这么低原因在于先验患病率太低假阳性数量压过了真阳性。这个例子在代码里的价值是提醒你贝叶斯更新的结果对先验极度敏感写代码时先验参数一定要有依据不能拍脑袋。条件概型还有个常见实现场景是马尔可夫链的转移概率。给定当前状态下一状态的概率只依赖当前状态这本质就是条件概型。代码里用转移矩阵表示transition np.array([[0.9, 0.1], [0.3, 0.7]]) # 从状态0出发两步后到状态1的概率 print((transition transition)[0, 1])矩阵乘法的顺序别搞反transition transition才是两步转移transition.T transition是另一回事。我见过有人把转移矩阵转置了还在纳闷结果不对这种错误在调试时特别隐蔽。2. 离散分布从伯努利到泊松的代码谱系离散分布处理的是可数结果的随机变量。伯努利是最小单元二项是伯努利的累加几何是第一次成功前的等待泊松是单位时间内的稀有事件计数。这四个分布之间有清晰的派生关系代码实现时理解这层关系比死记公式有用得多。2.1 伯努利与二项一次和多次的递进关系伯努利分布就是一次成败试验参数p是成功概率。代码实现简单到不需要库def bernoulli(p, rng): return 1 if rng.random() p else 0二项分布是n次独立伯努利试验中成功的次数。NumPy直接有rng.binomial(n, p, size)。但我想强调的是自己实现一遍的价值——用循环累加伯努利能直观看到二项是怎么长出来的def binomial_manual(n, p, rng): return sum(bernoulli(p, rng) for _ in range(n)) rng np.random.default_rng(42) samples [binomial_manual(10, 0.3, rng) for _ in range(10000)] print(np.mean(samples), np.var(samples)) # 约3.0, 约2.1理论均值np3方差np(1-p)2.1实测吻合。这里有个性能坑Python循环跑1万次还行跑100万次就慢了。生产环境直接用rng.binomial底层是C实现快两个数量级。二项分布代码里最需要警惕的是独立性假设。如果试验之间不独立比如抽样不放回二项就不适用得换成超几何分布。我见过做A/B测试的人直接用二项算显著性但用户行为往往有聚集性独立性不成立结论就偏了。2.2 几何分布等待首次成功的两种参数化陷阱几何分布描述首次成功所需的试验次数。这里有个大坑不同教材对几何分布的定义不同。有的定义成首次成功前的失败次数取值0,1,2,...有的定义成首次成功时的试验次数取值1,2,3,...。两者差1均值和方差公式也不同。# 定义A首次成功时的试验次数取值1,2,... def geometric_trials(p, rng): trials 0 while True: trials 1 if rng.random() p: return trials # 定义B首次成功前的失败次数取值0,1,... def geometric_failures(p, rng): return geometric_trials(p, rng) - 1NumPy的rng.geometric(p)用的是定义A取值从1开始。如果你按定义B推导了公式又用NumPy采样结果会系统性偏1。这个坑我在做可靠性分析时踩过——算平均无故障时间公式用的是失败次数版本采样用的是NumPy最后结果差了整整一个单位时间排查了半天才发现是参数化不一致。几何分布还有个无记忆性P(Xmn|Xm)P(Xn)。代码验证很直观p 0.2 # 已经等了5次还没成功再等3次仍未成功的条件概率 # 应该等于从头等3次仍未成功的概率无记忆性是几何分布独有的连续情形里指数分布也有这个性质在建模等待时间时非常有用但也意味着它不适合描述越等越可能成功的场景——那种场景得用别的分布。2.3 泊松分布稀有事件计数的均值方差相等特性泊松分布描述单位时间/空间内稀有事件的发生次数参数λ是平均发生率。它的标志性特征是均值等于方差都等于λ。这个特征在代码里可以用来做快速校验rng np.random.default_rng(7) lam 4.0 samples rng.poisson(lam, size100000) print(samples.mean(), samples.var()) # 都接近4.0泊松和二项的关系值得说清楚当n很大、p很小时二项分布近似泊松λnp。代码验证n, p 10000, 0.0004 binom_samples rng.binomial(n, p, size100000) poisson_samples rng.poisson(n * p, size100000) # 两者分布形态几乎重合这个近似在实际中很有用——当n大到无法逐次模拟时用泊松替代二项能大幅提速。但要注意近似条件一般要求n≥20且p≤0.05n越大p越小近似越好。我见过有人n50、p0.3还用泊松近似误差肉眼可见。泊松分布代码实现时有个容易忽略的点事件必须独立且发生率恒定。如果事件有聚集性比如网站访问量在促销期间暴增泊松就不适用得考虑负二项分布或复合泊松。做流量预测时直接套泊松遇到突发流量就会严重低估尾部风险。3. 连续分布均匀、指数、正态的实现细节连续分布处理的是不可数结果的随机变量概率通过密度函数积分得到。代码实现的核心是采样和密度计算但每个分布都有各自的数值陷阱。3.1 均匀分布看似最简单却最容易用错边界均匀分布U(a,b)的密度在[a,b]上是1/(b-a)外面是0。NumPy的rng.uniform(low, high, size)生成的是[low, high)上的均匀分布——注意是左闭右开。这个细节在离散化时会影响结果rng np.random.default_rng(0) # 生成[0,1)均匀数转成0-9的整数 digits np.floor(rng.random(100000) * 10).astype(int) # 各数字出现频率应接近0.1用np.floor而不是np.round因为rng.random()取不到1乘以10后范围是[0,10)floor后正好是0-9。如果用round0和9的出现概率会只有其他数字的一半因为边界只覆盖了半个区间。这个坑在做随机分组、随机抽样时特别常见。均匀分布还有个用途是逆变换采样——任何分布的采样都可以通过均匀分布累积分布反函数实现。以指数分布为例def exponential_inverse_transform(lam, size, rng): u rng.random(size) return -np.log(1 - u) / lam这就是逆变换采样的标准形态。理解这一点你就能自己实现任意分布的采样器不必依赖现成库。3.2 指数分布与泊松的孪生关系及无记忆性验证指数分布描述泊松过程中相邻事件的时间间隔参数λ是事件发生率。它和泊松是计数与间隔的孪生关系泊松数事件个数指数测事件间隔。rng np.random.default_rng(1) lam 2.0 intervals rng.exponential(1/lam, size100000) print(intervals.mean()) # 约0.5即1/lam注意NumPy的rng.exponential(scale)参数是scale1/λ不是λ本身。这个参数化方向和很多教材相反我第一次用时把λ直接传进去结果均值差了4倍。记忆方法scale是尺度单位和时间一样λ是速率单位是时间的倒数。指数分布的无记忆性可以用代码直观验证# P(X 23 | X 2) 应等于 P(X 3) lam 1.0 # 理论值 p_conditional np.exp(-lam * 3) # 仿真验证 samples rng.exponential(1/lam, size1_000_000) mask samples 2 print(np.mean(samples[mask] 5)) # 接近 exp(-3)无记忆性意味着已经等了多久不影响还要等多久这在建模设备寿命、客户到达间隔时是强假设。如果实际过程有老化或疲劳效应指数分布就不合适得用威布尔分布。3.3 正态分布中心极限定理的代码实证与参数陷阱正态分布是连续分布里用得最多的参数μ和σ分别控制位置和尺度。NumPy的rng.normal(loc, scale, size)里scale就是标准差σ不是方差。这个区分在写代码时经常被混淆rng np.random.default_rng(2) samples rng.normal(loc10, scale2, size100000) print(samples.mean(), samples.std()) # 约10, 约2正态分布最值得用代码验证的是中心极限定理大量独立同分布随机变量之和近似正态。拿均匀分布做实验rng np.random.default_rng(3) n, m 30, 10000 sums rng.random((m, n)).sum(axis1) # sums应近似正态均值n/215方差n/122.5n30时已经相当接近正态了。这个实验的价值在于让你对什么时候可以用正态近似有手感——n小于10时偏态还很明显n到30左右基本可用这也是很多统计方法默认样本量30的由来。正态分布代码里有个隐蔽陷阱用样本标准差估计σ时的自由度。np.std默认ddof0除以n而统计推断里通常要ddof1除以n-1才是无偏估计。做假设检验时用错ddof小样本下结论会偏data rng.normal(0, 1, size10) print(np.std(data, ddof0), np.std(data, ddof1))样本量小时两者差异明显样本量大时趋于一致。这个细节在写统计检验代码时务必注意。4. 分布选择与验证从数据反推分布的实操路径前面讲了各分布的实现但实际工作中更常见的需求是反过来的手上有一堆数据该用什么分布建模这一节讲我的实操路径——先看数据形态再做拟合检验最后用仿真验证。4.1 用直方图和矩特征做初步判断拿到数据第一步不是套分布而是看形态。直方图看偏度和峰度均值方差看离散程度。我通常写一个快速诊断函数import numpy as np from scipy import stats def diagnose(data): print(f均值: {np.mean(data):.3f}) print(f方差: {np.var(data):.3f}) print(f偏度: {stats.skew(data):.3f}) print(f峰度: {stats.kurtosis(data):.3f})判断逻辑均值≈方差且取值为非负整数优先考虑泊松偏度明显为正且取值连续非负考虑指数或对数正态偏度接近0、峰度接近0考虑正态。这个初判不能替代正式检验但能快速排除明显不合适的候选。我踩过的一个坑是只看均值和方差就下结论。有次数据均值方差都是5左右我以为是泊松结果直方图一看是双峰——明显是两个泊松的混合。所以直方图这一步不能省矩特征只能辅助。4.2 KS检验与拟合优度的代码模板初步锁定候选分布后用Kolmogorov-Smirnov检验做拟合优度。SciPy的stats.kstest可以直接用from scipy import stats data rng.poisson(4, size1000) # 检验是否服从泊松分布 stat, p_value stats.kstest(data, poisson, args(4,)) print(fKS统计量: {stat:.4f}, p值: {p_value:.4f})p值大于0.05表示不能拒绝数据服从该分布的原假设。但这里有个重要提醒KS检验对参数估计后的分布会偏保守因为参数是从数据估出来的自由度损失了。严谨做法是用Lilliefors检验或蒙特卡洛校正。实际工作中如果p值在0.05附近徘徊别急着下结论多换几个样本看看稳定性。还有个常见误区是把不能拒绝当成服从。KS检验不拒绝只说明数据与假设分布没有显著矛盾不代表就是那个分布。样本量越大检验越敏感微小偏差都会被检出。所以大样本下p值小不一定意味着分布选错了可能只是偏差在统计上显著但实际影响很小。4.3 采样一致性验证把拟合好的分布跑一遍仿真选定分布并估计参数后我习惯做一步回测用拟合的分布重新采样对比仿真数据和原始数据的关键统计量。这一步能暴露参数估计的偏差def backtest(data, dist_name, params, n_sim10000): sim getattr(rng, dist_name)(*params, sizen_sim) print(f原始均值: {np.mean(data):.3f}, 仿真均值: {np.mean(sim):.3f}) print(f原始方差: {np.var(data):.3f}, 仿真方差: {np.var(sim):.3f}) # 对比分位数 for q in [0.05, 0.5, 0.95]: print(fq{q}: 原始{np.quantile(data, q):.3f}, 仿真{np.quantile(sim, q):.3f})分位数对比比均值方差更能反映尾部行为。如果原始数据在95%分位数上明显高于仿真说明拟合分布低估了尾部风险——这在金融、可靠性领域是致命问题得换重尾分布如t分布、帕累托分布。回测还有个作用是发现参数估计方法的问题。比如用矩估计拟合正态遇到离群点时σ会被高估用最大似然估计则相对稳健。代码里scipy.stats.norm.fit默认用MLE但你可以手动指定mu_mle, sigma_mle stats.norm.fit(data, methodMM) # 矩估计不同估计方法在小样本或含噪数据下结果差异明显回测能帮你判断哪种更合适。5. 几个我反复踩过的实现坑与规避手法这一节单独拎出来讲因为下面这些问题不是某个分布特有的而是跨分布通用的且往往在调试时最难定位。5.1 随机种子管理与可复现性做概率仿真可复现性是底线。NumPy新版推荐用default_rng(seed)创建独立的Generator实例而不是全局的np.random.seed。区别在于全局种子会影响所有后续随机调用多个函数共享时互相干扰独立Generator各管各的互不影响。# 推荐独立Generator rng1 np.random.default_rng(42) rng2 np.random.default_rng(42) # rng1和rng2产生相同序列但互不干扰 # 不推荐全局种子 np.random.seed(42)做并行仿真时给每个进程分配不同的种子偏移避免所有进程产生相同序列from concurrent.futures import ProcessPoolExecutor def simulate(seed): rng np.random.default_rng(seed) return rng.normal(size1000).mean() with ProcessPoolExecutor() as ex: results list(ex.map(simulate, range(8)))如果所有进程都用同一个种子8个结果会完全一样仿真就白做了。这个坑我在第一次写并行蒙特卡洛时踩得结结实实。5.2 浮点精度与极端概率的处理概率计算里浮点精度问题在尾部概率上尤其突出。比如计算正态分布在极端值处的概率直接算会下溢到0from scipy import stats # 直接算 print(stats.norm.pdf(40)) # 下溢为0 # 用对数密度 print(stats.norm.logpdf(40)) # 有限值做似然计算、贝叶斯推断时一律用对数域运算最后再取指数。这个习惯能避免大量数值问题。另一个精度坑是累积概率的相减。P(aXb)F(b)-F(a)当a和b很接近时两个接近1的数相减会丢有效数字。稳妥做法是用stats.norm.cdf(b) - stats.norm.cdf(a)在b-a较小时改用stats.norm.sf(a) - stats.norm.sf(b)sf是生存函数算的是右尾或者直接用区间概率函数。5.3 分布参数化的方向混淆前面提过指数分布的scale和rate其实还有几个分布有类似问题分布NumPy参数常见教材参数易错点指数scale1/λλ速率传错方向差λ²倍正态scaleσσ标准差误传方差伽马scaleθ, shapek有时用rate1/θ方向反了贝塔a, b有时用均值和集中度参数含义不同规避方法写代码时在变量名里带上单位或含义比如rate_lambda、scale_theta别用裸的p、a、b。这个习惯看着啰嗦但能省下大量调试时间。5.4 大样本下的性能与内存做百万级采样时内存和速度都要考虑。rng.normal(size10_000_000)一次生成1000万个浮点数占80MB内存还行。但如果要生成多个分布做对比累积起来就吃紧了。技巧是分块生成、流式处理def streaming_mean(dist_func, total, chunk100000): acc 0.0 for _ in range(total // chunk): acc dist_func(chunk).sum() return acc / total这样内存占用恒定在chunk级别适合超大规模仿真。另外能用向量化就别用循环NumPy的向量化操作比Python循环快一到两个数量级这个差距在大样本下非常明显。概率分布的代码实现说到底就是把数学定义翻译成可执行、可验证、可复现的计算过程。我个人的体会是每学一个分布别急着背公式先用代码把它采样出来、画出直方图、算一遍均值和方差和理论值对一对。这个采样—验证的闭环走多了你对每个分布的直觉会越来越准遇到新场景时判断该用哪个分布也会越来越快。至于那些参数化方向、边界条件、精度陷阱都是在一次次跑代码、对结果的过程中慢慢刻进肌肉记忆的——光看是记不住的得亲手跑一遍。
