概率公式工程落地：从期望方差到贝叶斯与分布采样实战

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(f"n={n:>10}: pi ≈ {pi_est:.6f}, 误差 {abs(pi_est - 3.1415926):.6f}")

我实际跑下来的结果很有代表性：

这个例子想说明两件事。第一，收敛速度确实是“样本量开平方级别”——样本量变成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, key=scores.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块磁盘故障的概率是多少？

这里如果用二项分布，就是 n=1000, p=0.005, k=10；如果使用泊松近似，λ = 1000 * 0.005 = 5，按泊松分布算 k=10。两种结果几乎一致，误差在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(mean=0.0, stddev=1.0): z, _ = box_muller() return mean + stddev * z

这里有个隐藏的坑：u1如果刚好取到0，log(0)会直接报错。Python的random.random()返回 [0,1)，虽然取到0的概率极低，但在跑大规模模拟时，2的几十亿次方样本量下，这种边界情况是真实可能遇到的。稳妥做法是强制跳过0：

def 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): # 使用A&S 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 # 示例：任务 a=3天, m=5天, b=9天 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，两个接近的数相减会灾难性抵消。这两类问题我在前三节已经提到过对应解法，这里再总结成一张表，方便以后对号入座：

这一节想强调的核心理念是：概率公式写起来很优雅，但落到代码里，每一步都要问自己“这个数值运算稳定吗”。写概率相关代码，默认就要用对数、用增量式更新、用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反推法采样泊松分布，预计算了一个概率表到 k=100，然后用均匀随机数在前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里跑这类测试时，我会刻意固定随机种子，保证同一段代码每次得到同样的样本序列，这样测试结果完全可复现。等到准备上线前，再换不同的随机种子做多轮回归，用“多轮结果都稳定”来增强信心。

这一套组合拳打下来，概率公式才能真正变成工程上靠得住的基础设施。哪怕只是一段几十行的采样函数，背后涉及的数值稳定性、随机数质量、边界条件、测试策略，每一项都关系到最终结果是否可信。写代码的人如果把这些细节都拿捏到位，后续用这些工具做性能预测、做模拟实验、做数据建模，才有了站得住的根基。