Python实战：用NumPy和SciPy快速计算8种概率分布的期望与方差

Python实战用NumPy和SciPy快速计算8种概率分布的期望与方差在数据分析与建模中概率分布的计算往往让初学者感到头疼。传统统计学教材中复杂的积分推导和公式记忆在实际项目中反而成了效率的绊脚石。本文将用Python的NumPy和SciPy两大神器带你用代码直接获取均匀分布、泊松分布等8种常见概率分布的关键参数把数学理论转化为可执行的解决方案。1. 环境配置与基础工具开始前需要确保已安装科学计算三件套。如果使用Anaconda这些库通常已预装。可以通过以下命令检查版本import numpy as np import scipy.stats as stats print(fNumPy版本: {np.__version__}) # 推荐≥1.20 print(fSciPy版本: {stats.__version__}) # 推荐≥1.7关键工具说明NumPy.random生成各种分布的随机样本SciPy.stats提供完整的分布对象包含expect()和var()等现成方法提示Jupyter Notebook中可使用%timeit测试代码执行效率这对大数据量下的分布计算尤为重要2. 连续型分布实战2.1 均匀分布Uniform均匀分布是入门最简单的连续分布描述区间[a,b]内等概率出现的事件。用SciPy计算其期望和方差a, b 2, 5 uniform stats.uniform(loca, scaleb-a) # loc起点, scale区间长度 # 理论值 print(f理论期望: {(ab)/2:.2f}) print(f理论方差: {(b-a)**2/12:.2f}) # 计算值 print(f计算期望: {uniform.expect():.2f}) print(f计算方差: {uniform.var():.2f})实际项目中我们常需要验证数据是否符合均匀分布。这时可以结合直方图与K-S检验samples uniform.rvs(size1000) # 生成1000个样本 hist np.histogram(samples, bins20) stats.kstest(samples, uniform, args(a, b-a)) # p值0.05则接受均匀分布假设2.2 正态分布Normal正态分布参数计算最为直观因其期望和方差直接对应分布参数mu, sigma 0, 1 normal stats.norm(locmu, scalesigma) print(f样本期望: {normal.mean():.2f}) # 直接属性访问 print(f样本方差: {normal.var():.2f})对于非标准正态分布可以通过矩估计法从样本反推参数data np.random.normal(5, 3, 1000) # 生成μ5,σ3的样本 estimated_mu np.mean(data) # 样本均值估计μ estimated_sigma np.std(data, ddof1) # 样本标准差估计σ3. 离散型分布实战3.1 泊松分布Poisson泊松分布常用于描述单位时间内随机事件发生的次数。其特点是期望与方差相等lambda_ 4 # 事件发生率 poisson stats.poisson(mulambda_) print(f计算期望: {poisson.mean()}) print(f计算方差: {poisson.var()}) # 生成符合泊松分布的随机数 samples poisson.rvs(size1000) print(f样本均值: {np.mean(samples):.2f}) # 应接近lambda_3.2 二项分布Binomial二项分布是n次独立伯努利试验的成功次数分布。计算其参数时需注意概率p的范围n, p 10, 0.3 binomial stats.binom(nn, pp) print(f理论期望: {n*p}) print(fSciPy计算方差: {binomial.var()}) # 蒙特卡洛模拟验证 trials 10000 successes np.random.binomial(n, p, trials) print(f模拟方差: {np.var(successes):.2f})4. 特殊分布处理技巧4.1 指数分布的无记忆性指数分布的无记忆性使其在生存分析中广泛应用。验证这一特性lambda_ 0.5 exponential stats.expon(scale1/lambda_) # scale1/λ # 验证P(X10|X5) P(X5) cond_prob np.mean(exponential.rvs(10000) 10) / np.mean(exponential.rvs(10000) 5) print(f条件概率: {cond_prob:.4f} vs 理论值: {np.exp(-lambda_*5):.4f})4.2 超几何分布的抽样验证超几何分布描述不放回抽样中的成功次数。对比其与二项分布的差异N, K, n 100, 30, 10 # 总量成功数抽样数 hypergeom stats.hypergeom(MN, nK, Nn) print(f期望: {hypergeom.mean()}) print(f方差: {hypergeom.var()}) # 与二项分布对比 binom_var n * (K/N) * (1-K/N) print(f二项分布方差(近似): {binom_var:.2f}) print(f有限校正因子: {(N-n)/(N-1):.2f})5. 工程实践中的注意事项5.1 数值稳定性问题计算极端参数下的分布时可能会遇到数值溢出问题。例如泊松分布在大λ值时large_lambda 1e6 poisson_large stats.poisson(mularge_lambda) # 使用对数空间计算避免溢出 log_pmf poisson_large.logpmf(int(large_lambda)) print(fλ1e6时k1e6的对数概率: {log_pmf:.2f})5.2 分布选择建议根据数据特征选择合适分布的快速判断指南数据特征候选分布SciPy实现类连续、对称、钟形正态分布stats.norm连续、非负、右偏指数分布stats.expon离散、计数、稀有事件泊松分布stats.poisson离散、固定试验次数二项分布stats.binom连续、有界区间均匀分布stats.uniform5.3 自定义分布实现当内置分布不满足需求时可以通过子类化rv_continuous或rv_discrete创建自定义分布。以下实现一个三角分布from scipy.stats import rv_continuous class my_triangular(rv_continuous): def _pdf(self, x): return np.where(x 0.5, 2*x, 2*(1-x)) custom_dist my_triangular(a0, b1) print(f自定义分布期望: {custom_dist.expect():.2f})6. 性能优化技巧6.1 向量化计算利用NumPy的向量化操作加速大批量计算# 低效的循环计算 def slow_cdf(x, dist): return np.array([dist.cdf(val) for val in x]) # 高效的向量化计算 x np.linspace(0, 1, 100000) %timeit slow_cdf(x, stats.norm()) # 约120ms %timeit stats.norm().cdf(x) # 约2ms6.2 并行计算对于需要重复生成大量样本的场景可以使用joblib并行化from joblib import Parallel, delayed def simulate(n): return np.mean(stats.norm().rvs(n)) results Parallel(n_jobs4)(delayed(simulate)(1000) for _ in range(10000))7. 可视化验证图形化展示是验证分布特性的重要手段。使用Matplotlib绘制关键分布import matplotlib.pyplot as plt distributions { Normal: stats.norm(0, 1), Exponential: stats.expon(scale1), Poisson: stats.poisson(5) } fig, axes plt.subplots(1, 3, figsize(15, 4)) for ax, (name, dist) in zip(axes, distributions.items()): if hasattr(dist, pdf): x np.linspace(dist.ppf(0.01), dist.ppf(0.99), 100) ax.plot(x, dist.pdf(x)) else: x np.arange(dist.ppf(0.01), dist.ppf(0.99)) ax.bar(x, dist.pmf(x)) ax.set_title(f{name} Distribution) plt.tight_layout()8. 实际案例客户到达时间分析假设某客服中心客户到达服从泊松过程每小时平均到达5人。我们需要计算15分钟内无人到达的概率连续两位客户间隔时间超过30分钟的概率# 问题1泊松分布计算 lambda_15min 5/4 # 将小时率转换为15分钟率 prob_no_arrival stats.poisson(mulambda_15min).pmf(0) # 问题2指数分布计算 lambda_hour 5 prob_interval_30min 1 - stats.expon(scale1/lambda_hour).cdf(0.5) print(f15分钟无人概率: {prob_no_arrival:.2%}) print(f间隔超30分钟概率: {prob_interval_30min:.2%})