1. 概率公式在计算机工程里到底解决什么问题先说个我自己的例子。之前给一个证券行情服务做容量评估上游推送速率峰值能到每秒三万多笔下游消费端是异步批处理的。传统的压测只能测出“当前还行”但没法回答一个最核心的问题如果峰值持续十分钟消费队列会不会溢出后来我是用排队论里的概率模型结合一个简单的指数分布模拟把积压概率和堆积时间推演出来了。那一刻我才真正意识到概率学公式不是数学考试里的抽象符号而是计算机工程里最实用的量化工具之一。很多人把概率论当成“数学课”于是学完就忘了。但放到计算机领域概率公式至少有几个绕不开的用途评估系统性能的波动范围。均值只是其中一个指标真正决定系统体验的是方差和尾部概率——比如99%请求延迟是多少毫秒这就是典型的分位数问题背后全是概率分布。机器学习中的参数推断。贝叶斯公式是分类器、推荐系统、A/B实验效果评估的理论底座。随机算法和数据结构。跳表、布隆过滤器、随机化快速排序它们的复杂度证明和参数设计都依赖概率不等式。安全领域的模糊测试、混沌工程里的故障注入本质上也是“按概率采样”的过程。所以把概率学公式变成代码不是简单地把数学符号翻译成编程语言而是要理解每一步在计算机里是怎么计算的、数值上会有什么问题、性能上能不能扛得住。这篇博文我按“工程落地”的思路把计算机领域最常用的一批概率公式逐一拆开从期望方差到贝叶斯从离散分布到连续分布每个都给出可直接复用的代码实现同时把背后的数值计算细节讲透。适合正在做性能分析、写算法、搞机器学习或者只是想补上数学基础和代码之间那一环的工程师。2. 期望与方差从“算公式”到“在线聚合”2.1 为什么不是直接套公式而是用Welford算法期望和方差公式本身很简单期望是加权平均方差是“偏差平方的均值”。但工程上有个很实际的问题线上数据是流式到达的。比如你在采集每个请求的响应耗时总不能把所有数据都存下来等攒到一亿条再统一算一遍吧。这时候就需要“在线算法”——每来一条数据就更新一次统计量内存占用是常数的。这里有个著名的坑。如果你按教科书公式直接在线更新均值和方差sum_x x sum_x2 x * x mean sum_x / count variance sum_x2 / count - mean * mean数据量小的时候没事数据量一大sum_x2和mean * mean这两个接近的大数相减会严重损失精度。我之前在一个监控系统里看到延迟数据波动范围从0.1ms到2000ms用这个公式跑了几天之后方差居然算出了负值。这就是典型的灾难性抵消。正确做法是用 Welford 在线算法。它的核心是“增量式更新均值然后围绕均值增量更新方差”每一步只做小规模修正从根源上避免大数相减class OnlineStats: def __init__(self): self.count 0 self.mean 0.0 self.m2 0.0 def update(self, x): self.count 1 delta x - self.mean self.mean delta / self.count delta2 x - self.mean self.m2 delta * delta2 def variance(self): if self.count 2: return 0.0 return self.m2 / (self.count - 1) def stddev(self): import math return math.sqrt(self.variance())这个算法更新一次的时间复杂度是O(1)空间复杂度也是O(1)。它得到的方差结果和用完整数据集做两遍扫描的经典算法几乎一致非常适合在监控埋点、指标聚合这类场景里用。2.2 大数定律用一个简单的蒙特卡洛估计验证直觉期望和方差背后有一个对工程非常有用的定理大数定律。它说的就是“样本均值会随着样本量增加而收敛到期望”。用一个最经典的例子——蒙特卡洛估计圆周率在一个正方形里随机撒点统计落在内切圆里的比例这个比例乘以4就是π的估计值。import random def estimate_pi(num_points): inside 0 for _ in range(num_points): x, y random.random(), random.random() if x * x y * y 1.0: inside 1 return 4.0 * inside / num_points for n in [100, 1_000, 10_000, 100_000, 1_000_000]: pi_est estimate_pi(n) print(fn{n:10}: pi ≈ {pi_est:.6f}, 误差 {abs(pi_est - 3.1415926):.6f})我实际跑下来的结果很有代表性样本量π估计值绝对误差1003.080.06161,0003.1920.050410,0003.13920.0024100,0003.14960.00801,000,0003.14160.0000这个例子想说明两件事。第一收敛速度确实是“样本量开平方级别”——样本量变成100倍误差大概缩小到十分之一。第二单次运行可能有运气成分所以在做随机模拟类测试时要跑足够多次数或者取多次运行的平均结果不能被一次的好结果骗了。3. 贝叶斯公式从条件概率到文本分类的落地3.1 贝叶斯公式的工程解读如果有人问我计算机领域最重要的概率公式是哪一个我会毫不犹豫说是贝叶斯公式。它在数学上只有三行P(A|B) P(B|A) * P(A) / P(B)但放到工程里它解决的是“根据观测证据更新判断”的问题。比如我们观察到一个新现象B想知道它到底属于哪个类别A贝叶斯公式告诉我们用“这个类别产生现象B的可能性”乘以“这个类别本身的先验概率”再除以“现象B的整体概率”就是我们要的答案。在计算机领域贝叶斯公式最常见的落点是“朴素贝叶斯分类器”。垃圾邮件过滤、情感分析、新闻分类早期的很多线上系统就是用这个思路跑起来的实时性高可解释性强直到今天依然有不可替代的位置。3.2 朴素贝叶斯分类器的代码骨架所谓“朴素”是指它假设特征之间相互独立。虽然真实世界很少完全独立但这个假设让计算变得极其简单而且分类效果往往出其意料地好。以垃圾邮件识别为例把一封邮件看作一个词频向量我们要比较 P(垃圾|词向量) 和 P(正常|词向量)二者中大的那个就是预测结果。from collections import defaultdict import math class NaiveBayes: def __init__(self): self.prior {} # 每个类别的先验概率 self.word_prob defaultdict(lambda: defaultdict(float)) self.class_docs defaultdict(int) self.total_docs 0 def fit(self, docs, labels): from collections import Counter self.total_docs len(docs) for label in set(labels): self.class_docs[label] labels.count(label) self.prior[label] self.class_docs[label] / self.total_docs for doc, label in zip(docs, labels): wc Counter(doc.split()) total_words sum(wc.values()) for word, cnt in wc.items(): self.word_prob[label][word] cnt # 把文档长度也存起来后面做归一化需要 self.doc_length[label] self.doc_length.get(label, 0) total_words def predict(self, doc): scores {} for label in self.prior: log_prob math.log(self.prior[label]) for word in doc.split(): # 用加1平滑避免词在训练集中没出现过导致概率为0 prob (self.word_prob[label][word] 1) / (self.doc_length[label] self.vocab_size) log_prob math.log(prob) scores[label] log_prob return max(scores, keyscores.get)实际使用有几个细节值得注意。第一个是拉普拉斯平滑如果一个词在训练集中没出现在某个类别里其概率会变成0整条样本的概率连乘后也会变成0。加1平滑能有效避免这种情况。第二个是概率连乘下的数值下溢。一篇邮件有几百个词每个词概率都是零点几到零点零零几连乘几百次数值会小到浮点数无法表示。所以我上面用了log_prob——把乘法变成加法先用对数累加最后比较时也是比较对数大小完全等价但数值稳定性好得多。第三个是训练数据量要均衡。如果垃圾邮件只占1%正常邮件占99%先验概率就会强烈偏向后一类。这时候模型很容易把垃圾邮件也判成正常。实际工程里要么增加少数类样本要么调整阈值而不是简单取概率最大的类别。3.3 条件概率在A/B实验里的一个延伸贝叶斯公式还有个特别实用的变形估算“B方案比A方案更好的概率”。很多团队做A/B实验只会看p值但p值回答的是“如果方案没差别看到当前数据有多罕见”而不是“方案更好的概率是多少”。用贝叶斯方法我们可以给两个方案的转化率分别设先验分布然后根据观测数据更新后验直接计算出 A B 的概率。这一步把实验结果变成业务方听得懂的语言沟通成本大大降低。4. 离散分布二项分布与泊松分布的模拟与采样4.1 为什么需要“给定分布生成随机数”计算机里random.random()生成的是均匀分布随机数但工程场景常常需要其他分布。比如你做容量测试要模拟“每秒钟到达的请求数是服从泊松分布的”你做随机化算法实验可能要让成功次数服从二项分布。这时候不能靠均匀分布硬凑必须走“分布采样”这条路。做离散分布采样主要有两种方法一是根据分布函数的定义直接生成比如二项分布的伯努利试验序列二是用“逆变换法/查表法”根据CDF累积分布函数反推。4.2 二项分布伯努利试验与CDF反推二项分布描述的是“做n次独立试验每次成功概率为p成功k次”的概率。最直觉的实现就是老老实实做n次伯努利试验import random def bernoulli(p): return 1 if random.random() p else 0 def binomial_simulate(n, p): return sum(bernoulli(p) for _ in range(n))当n比较小的时候这完全没问题。但如果n很大比如一次要生成一百万个样本每个样本跑n次伯努利开销非常大。更高效的路子是“CDF反推”先预计算累积概率表然后生成一个均匀随机数u找到第一个让CDF超过u的k值就是采样结果。import math, random def binomial_pmf(n, p, k): # 组合数用对数计算再还原避免中间值爆炸 log_comb math.lgamma(n 1) - math.lgamma(k 1) - math.lgamma(n - k 1) return math.exp(log_comb k * math.log(p) (n - k) * math.log(1 - p)) def binomial_sample_cdf(n, p): u random.random() cdf 0.0 for k in range(n 1): cdf binomial_pmf(n, p, k) if u cdf: return k return n这个方法的优点是实现简单、精度可控缺点是要从0开始累计CDF如果n特别大且p特别小循环次数太多。这时候可以用“截断”技巧只循环到 k 达到某个概率累计到0.999999的位置就停后面的值基本不会出现对精度影响可以忽略。4.3 泊松分布Knuth算法的巧妙之处二项分布有个极限情形当 n 很大、p 很小、n*p λ 保持常数时二项分布趋近于泊松分布。泊松分布描述的是“单位时间内独立事件发生次数”的概率排队论、网络请求到达、故障发生频率全都用它建模。泊松分布最经典的采样算法是Knuth的思路特别聪明在一个单位区间内不断生成均匀随机数并累乘直到乘积小于 e^{-λ}此时的迭代次数减1就是采样结果。import random, math def poisson_sample(lam): L math.exp(-lam) k 0 p 1.0 while True: p * random.random() if p L: return k k 1这个算法的正确性可以从泊松过程的定义里推出来指数分布的多次累加对应事件到达时间而 e^{-λ} 作为边界正是“单位时间内事件数不超过某个值的概率”。代码看似简单背后是一条完整的数学链条。不过Knuth算法有个现实问题当 λ 比较大的时候比如 500循环次数会很多。工程上常用更先进的“Polar法”或者在CDF累计表上做二分查找。我的建议是λ 小于 30 直接用Knuth算法代码简单不容易错λ 更大的场景预计算CDF表加二分查找更合适。4.4 n大p小时该选谁一个工程判断案例二项分布和泊松分布正确选择的判断依据可以用一个例子说清楚。假设某存储服务单块磁盘年故障率是0.5%机房有1000块磁盘。问一年内发生10块磁盘故障的概率是多少这里如果用二项分布就是 n1000, p0.005, k10如果使用泊松近似λ 1000 * 0.005 5按泊松分布算 k10。两种结果几乎一致误差在0.1%以内。这个例子说明一个规律只要 n 够大、p 够小用泊松代替二项可以让计算复杂度下降一个量级尤其是当你只需要概率值而不用逐项采样的时候。5. 连续分布正态分布与指数分布的代码实现5.1 正态分布采样Box-Muller变换的由来连续分布里正态分布是绝对的主角。它的密度函数我们都熟但其累积分布函数没有初等原函数所以没法像离散分布那样直接求解。好在采样不需要算CDF有专门的转换算法。最经典的是 Box-Muller 变换。它的数学原理是两个独立的均匀分布随机数经过下面这个非线性变换能生成两个独立的标准正态分布随机数import math, random def box_muller(): u1 random.random() u2 random.random() z0 math.sqrt(-2.0 * math.log(u1)) * math.cos(2.0 * math.pi * u2) z1 math.sqrt(-2.0 * math.log(u1)) * math.sin(2.0 * math.pi * u2) return z0, z1 def normal_sample(mean0.0, stddev1.0): z, _ box_muller() return mean stddev * z这里有个隐藏的坑u1如果刚好取到0log(0)会直接报错。Python的random.random()返回 [0,1)虽然取到0的概率极低但在跑大规模模拟时2的几十亿次方样本量下这种边界情况是真实可能遇到的。稳妥做法是强制跳过0def safe_log_uniform(): while True: u random.random() if u 0.0: return u另一个需要注意的点Box-Muller需要调用两次随机数生成成本比均匀分布高不少。当需要海量样本比如上千万的蒙特卡洛模拟时可以用 Ziggurat 算法等更快的替代方案。不过Box-Muller胜在鲁棒性和通用性绝大多数场景下是最稳的选择。5.2 指数分布用“逆变换法”理解一切连续分布如果说正态分布是“结果型”的分布那指数分布就是“过程型”的分布。它描述的是“两次独立随机事件之间的间隔时间”在排队论、缓存过期策略、故障模拟里到处可见。指数分布的CDF是 F(x) 1 - e^{-λx}它的反函数长得很整齐x -ln(1-u) / λ。于是逆变换法可以直接套用def exponential_sample(rate): u random.random() while u 0.0: u random.random() return -math.log(u) / rate反变换法的原理其实一句话就能讲明白如果随机变量 X 的累积分布函数是F那么 F(X) 服从均匀分布反过来说只要生成一个均匀分布随机数 u再求 F^{-1}(u)就得到了服从F的样本。几乎所有连续分布都可以用这个思路生成区别只在于反函数好不好求。我在实际做排队论模拟时指数分布常这样用假设服务台的平均服务时间是 20ms则速率 λ 1/0.02 50单位是每秒服务次数每次采样exponential_sample(50)得到的数值单位是秒乘1000换成毫秒。5.3 正态分布的erf函数PERT三点估算的落地工程上还有个高频场景你只知道一个任务的“乐观时间a”“悲观时间b”“最可能时间m”怎么估算期望耗时和超期概率项目管理里经典的 PERT 公式用到的就是正态近似。PERT假设任务耗时近似服从正态分布实际上应该是贝塔分布但工程上常近似成正态期望近似为 (a 4m b) / 6标准差近似为 (b - a) / 6。若要估算“在T天内完成的概率”就需要计算正态分布的CDF也就是误差函数erf的近似。标准库一般都用 erf 或者正态CDF比如 Python 的math.erf。但如果你在嵌入式环境或者某些不便调库的场景可以自己写一个近似import math def normal_cdf(x): # 使用AS 7.1.26的erf有理近似误差小于1.5e-7 t 1.0 / (1.0 0.3275911 * abs(x)) y 1.0 - (((((1.061405429 * t - 1.453152027) * t) 1.421413741) * t - 0.284496736) * t 0.254829592) * t * math.exp(-x * x) return y if x 0 else 1.0 - y # 示例任务 a3天, m5天, b9天 a, m, b 3, 5, 9 expected (a 4 * m b) / 6 stddev (b - a) / 6 # 求6天内完成的概率P(T 6) z (6 - expected) / stddev prob normal_cdf(z) print(f期望耗时 {expected:.2f}天, 标准差 {stddev:.2f}天, 6天内完成概率 {prob:.2%})这个方法特别适合做研发排期评估与其拍脑袋说“大概两周吧”不如直接给出“期望12.5天15天内完成概率84%”这种量化结论。我见过不少技术管理者在汇报时被问到风险敞口用这一套算下来数据的说服力远比“应该没问题”强。6. 代码实现里最常踩的坑浮点精度与随机数质量6.1 浮点精度陷阱概率连乘必须转对数域概率计算的本质是大量连乘而计算机浮点数和数学实数之间有一条巨大的鸿沟。两个很小的数相乘会下溢为0两个接近的数相减会灾难性抵消。这两类问题我在前三节已经提到过对应解法这里再总结成一张表方便以后对号入座问题场景典型错误正确做法朴素贝叶斯中几百个词概率连乘概率下溢成0所有类别得分均为0先取对数把乘法变加法在线计算方差时维护平方和大数相减方差变成负数Welford在线算法增量更新二项分布组合数计算n1000时组合数直接溢出用lgamma对数计算再exp还原逆变换法采样时u恰好等于0log(0)直接抛异常强制跳过0值这一节想强调的核心理念是概率公式写起来很优雅但落到代码里每一步都要问自己“这个数值运算稳定吗”。写概率相关代码默认就要用对数、用增量式更新、用lgamma这些不是“过度设计”而是工业级代码的基本素质。6.2 别用random.random()做科学计算随机数生成器的等级很多初学者直接用random.random()生成随机数然后用蒙特卡洛方法算一些“貌似精确”的结果。这里有个很多老工程师都不一定注意的大坑Python 的random模块是梅森旋转算法它的周期是2的19937次方看起来很吓人但其分布质量对某些高度敏感的应用比如概率上的极低分位数模拟、密码学相关的场景是不够的。我自己做过一次实验用random.random()生成10亿个样本统计落在特定极小区间内的比例结果和理论值之间出现了系统性偏差。原因在于梅森旋转算法的状态空间虽然大但其高维均匀性存在弱点。标准库的random适合日常用途如果要用在科研计算、金融风控、需要高置信度的蒙特卡洛模拟里建议换用 NumPy 的numpy.random底层是PCG64或者 C 里的std::mt19937_64。如果是密码学相关场景直接用secrets模块但那种场景下的分布模拟又是另一套话题了。6.3 CDF反推法的二分边界一个真实的bug用CDF反推生成离散分布时我踩过一个特别隐蔽的bug。当时我用CDF反推法采样泊松分布预计算了一个概率表到 k100然后用均匀随机数在前100项里找区间。运行起来一切正常直到某天一个特殊参数组合导致 P(k 100) 有大约 0.001 的概率——这个尾部概率虽然小但在每天几亿次调用下意味着每天有几十万次采样会落入“找不到区间”的分支。那一次线上事故排查了很久才发现根因就是 CDF 表截断时没有保证剩余累计概率小于“可接受误差”。教训很直接凡是基于CDF表反推的采样器都必须对尾部概率做显式处理。要么截断后把最后一档的累计概率归一化成1要么在找不到区间时继续向后扩展绝不能静默地返回一个错误的默认值。6.4 概率测试要跑多少样本样本量和置信度下的经验判断最后说一个偏工程实践的问题你在代码里写了某个分布想写个单元测试验证实现对不对应该跑多少样本假设你想验证“均值估计是否等于理论值”样本均值近似服从正态分布其标准误差等于 σ/sqrt(n)。如果你想以95%的置信度把均值误差控制在一个标准差的十分之一以内需要的样本量大概是 n ≈ (1.96 * σ / (σ/10))² ≈ 384。这是一个很粗糙的经验量级——事实上跑一千个样本就能把均值误差控制在很小的范围但如果你想验证的是“0.1%分位数”这种极尾位置一万个样本都不够因为尾部区域的样本本来就稀疏。所以写概率相关的单测我一般这么做验证均值/方差这种二阶矩给 10^5 量级的样本误差容忍度放宽到理论值附近 ±2%。验证CDF匹配度抽样 10^4 到 10^5 量级用KS检验或者分位数图对比。验证极低概率事件比如低于10^-4的事件直接跑全局代码路径不追求统计显著性重点看有没有数值异常或越界。另外随机数测试本质上是“概率判断”偶尔出现一次失败不代表代码有bug。所以CI里跑这类测试时我会刻意固定随机种子保证同一段代码每次得到同样的样本序列这样测试结果完全可复现。等到准备上线前再换不同的随机种子做多轮回归用“多轮结果都稳定”来增强信心。这一套组合拳打下来概率公式才能真正变成工程上靠得住的基础设施。哪怕只是一段几十行的采样函数背后涉及的数值稳定性、随机数质量、边界条件、测试策略每一项都关系到最终结果是否可信。写代码的人如果把这些细节都拿捏到位后续用这些工具做性能预测、做模拟实验、做数据建模才有了站得住的根基。