从“杯子放球”到“射击命中”用Python模拟帮你彻底搞懂离散随机变量概率论中的离散随机变量概念常常让初学者感到抽象难懂。传统的数学推导虽然严谨但缺乏直观性。本文将带你用Python代码亲手模拟几个经典概率问题通过可视化手段让这些概念变得触手可及。1. 环境准备与基础概念在开始之前我们需要配置Python环境并安装必要的库。推荐使用Jupyter Notebook进行交互式编程它能实时显示代码运行结果和可视化图表。# 安装必要库 !pip install numpy matplotlib离散随机变量的核心特征是它只能取有限或可数个值。比如掷骰子的结果1到6抛硬币的结果正面或反面射击命中次数0到n关键区别与连续随机变量不同离散随机变量的概率分布可以用概率质量函数(PMF)完整描述。2. 经典问题模拟杯子放球实验这个经典问题描述为将n个球随机放入N个杯子中每个杯子被选中的概率相同。我们关心的是特定杯子中球的数量。2.1 问题建模与代码实现import numpy as np import matplotlib.pyplot as plt def cup_balls_experiment(n_balls10, n_cups5, n_trials10000): 模拟杯子放球实验 :param n_balls: 球的数量 :param n_cups: 杯子的数量 :param n_trials: 实验次数 :return: 特定杯子中球数的概率分布 results np.zeros(n_trials) for i in range(n_trials): # 随机分配球到杯子 allocations np.random.randint(0, n_cups, sizen_balls) # 统计第一个杯子中的球数 results[i] np.sum(allocations 0) # 计算概率分布 unique, counts np.unique(results, return_countsTrue) probabilities counts / n_trials return unique, probabilities2.2 结果可视化与分析运行上述代码并绘制结果values, probs cup_balls_experiment(n_balls10, n_cups5) plt.bar(values, probs) plt.xlabel(Number of balls in the first cup) plt.ylabel(Probability) plt.title(Probability Distribution of Balls in a Cup) plt.show()理论分析这实际上是一个二项分布问题每个球独立地以概率1/N落入特定杯子。当N5时理论概率应为$$ P(k) C_{10}^k \left(\frac{1}{5}\right)^k \left(\frac{4}{5}\right)^{10-k} $$3. 伯努利试验与射击命中问题射击命中问题是典型的伯努利试验每次射击独立命中概率为p未命中概率为1-p。n次射击中命中次数k服从二项分布。3.1 模拟实现def binomial_simulation(n10, p0.3, trials10000): 模拟二项分布 :param n: 试验次数 :param p: 成功概率 :param trials: 模拟次数 :return: 成功次数的概率分布 results np.random.binomial(n, p, trials) unique, counts np.unique(results, return_countsTrue) probabilities counts / trials return unique, probabilities3.2 不同参数下的对比我们可以比较不同命中概率p对分布的影响p值分布形状特征期望值方差0.1右偏集中在0附近10.90.5对称分布52.50.9左偏集中在n附近90.9# 比较不同p值的分布 p_values [0.1, 0.5, 0.9] plt.figure(figsize(12,4)) for i, p in enumerate(p_values): plt.subplot(1,3,i1) values, probs binomial_simulation(pp) plt.bar(values, probs) plt.title(fp{p}) plt.xlabel(Number of successes) plt.ylabel(Probability) plt.tight_layout() plt.show()4. 泊松分布稀有事件建模泊松分布常用于描述单位时间内稀有事件发生的次数。例如网站每分钟的访问量放射性物质单位时间的衰变次数4.1 泊松过程模拟def poisson_process(rate3, time1, trials10000): 模拟泊松过程 :param rate: 事件发生率 :param time: 观察时间 :param trials: 模拟次数 :return: 事件次数的概率分布 results np.random.poisson(rate*time, trials) unique, counts np.unique(results, return_countsTrue) probabilities counts / trials return unique, probabilities4.2 泊松与二项分布的关系当n很大而p很小时二项分布近似于泊松分布。我们可以通过模拟验证这一点# 比较二项分布和泊松分布 n, p 1000, 0.003 # λ np 3 binomial_values, binomial_probs binomial_simulation(nn, pp) poisson_values, poisson_probs poisson_process(raten*p) plt.bar(binomial_values-0.1, binomial_probs, width0.2, labelBinomial) plt.bar(poisson_values0.1, poisson_probs, width0.2, labelPoisson) plt.legend() plt.xlabel(Number of events) plt.ylabel(Probability) plt.title(Comparison of Binomial and Poisson Distributions) plt.show()5. 几何分布等待首次成功几何分布描述了在独立伯努利试验中首次成功所需的试验次数。例如射击直到首次命中打电话直到首次接通5.1 模拟实现def geometric_simulation(p0.2, trials10000): 模拟几何分布 :param p: 每次试验成功概率 :param trials: 模拟次数 :return: 首次成功所需试验次数的概率分布 results [] for _ in range(trials): count 1 while np.random.random() p: count 1 results.append(count) unique, counts np.unique(results, return_countsTrue) probabilities counts / trials return unique, probabilities5.2 无记忆性验证几何分布的一个重要特性是无记忆性已经失败了k次后首次成功还需要等待的次数的分布与原始分布相同我们可以用代码验证这一特性# 验证无记忆性 p 0.1 values, probs geometric_simulation(pp) # 条件概率已知前5次都失败首次成功在第6次及以后的分布 conditional_results [] for _ in range(10000): count 1 # 确保前5次都失败 while True: trials [] for _ in range(5): trials.append(np.random.random() p) if all(trials): break # 从第6次开始记录首次成功 count 6 while np.random.random() p: count 1 conditional_results.append(count-5) # 额外的等待次数 # 比较原始分布和条件分布 cond_unique, cond_counts np.unique(conditional_results, return_countsTrue) cond_probs cond_counts / sum(cond_counts) plt.bar(values, probs, alpha0.5, labelOriginal) plt.bar(cond_unique, cond_probs, alpha0.5, labelConditional) plt.legend() plt.xlabel(Additional waiting time) plt.ylabel(Probability) plt.title(Memoryless Property Verification) plt.show()6. 负二项分布多次成功的等待时间负二项分布是几何分布的推广描述了获得r次成功所需的试验次数。6.1 模拟实现def negative_binomial_simulation(r3, p0.2, trials10000): 模拟负二项分布 :param r: 需要的成功次数 :param p: 每次试验成功概率 :param trials: 模拟次数 :return: 达到r次成功所需试验次数的概率分布 results [] for _ in range(trials): successes 0 count 0 while successes r: count 1 if np.random.random() p: successes 1 results.append(count) unique, counts np.unique(results, return_countsTrue) probabilities counts / trials return unique, probabilities6.2 与几何分布的关系当r1时负二项分布退化为几何分布。我们可以验证这一点# 比较r1的负二项分布和几何分布 nb_values, nb_probs negative_binomial_simulation(r1, p0.2) geo_values, geo_probs geometric_simulation(p0.2) plt.bar(nb_values-0.1, nb_probs, width0.2, labelNegative Binomial (r1)) plt.bar(geo_values0.1, geo_probs, width0.2, labelGeometric) plt.legend() plt.xlabel(Number of trials until first success) plt.ylabel(Probability) plt.title(Negative Binomial vs Geometric (r1)) plt.show()7. 超几何分布不放回抽样超几何分布描述了在不放回抽样中特定类别的物品被抽中的次数。例如从包含5个红球和15个黑球的袋子中抽取10个球红球的数量质量控制中的缺陷品抽样检查7.1 模拟实现def hypergeometric_simulation(N20, K5, n10, trials10000): 模拟超几何分布 :param N: 总体大小 :param K: 成功物品数量 :param n: 抽样数量 :param trials: 模拟次数 :return: 样本中成功物品数量的概率分布 population np.array([1]*K [0]*(N-K)) results [] for _ in range(trials): sample np.random.choice(population, sizen, replaceFalse) results.append(np.sum(sample)) unique, counts np.unique(results, return_countsTrue) probabilities counts / trials return unique, probabilities7.2 与二项分布的比较当N很大而n相对较小时超几何分布近似于二项分布。我们可以观察这一现象# 比较超几何分布和二项分布 N, K, n 1000, 200, 50 # p K/N 0.2 hyp_values, hyp_probs hypergeometric_simulation(NN, KK, nn) bin_values, bin_probs binomial_simulation(nn, pK/N) plt.bar(hyp_values-0.1, hyp_probs, width0.2, labelHypergeometric) plt.bar(bin_values0.1, bin_probs, width0.2, labelBinomial) plt.legend() plt.xlabel(Number of successes in sample) plt.ylabel(Probability) plt.title(Hypergeometric vs Binomial (Large N)) plt.show()8. 实际应用案例质量控制模拟让我们将这些分布应用于一个实际的质量控制场景。假设一个工厂生产灯泡历史数据显示有5%的缺陷率。质检部门每天随机抽取50个灯泡进行测试。8.1 问题建模我们可以用二项分布建模每天的缺陷品数量def quality_control_simulation(defect_rate0.05, sample_size50, days365): 模拟一年的质量控制过程 :param defect_rate: 缺陷率 :param sample_size: 每日抽样数量 :param days: 模拟天数 :return: 每日缺陷品数量 daily_defects np.random.binomial(sample_size, defect_rate, days) return daily_defects8.2 异常检测我们可以计算均值和标准差并设置3σ控制限defects quality_control_simulation() mean np.mean(defects) std np.std(defects) plt.plot(defects, b-, labelDaily defects) plt.axhline(mean, colorg, linestyle--, labelMean) plt.axhline(mean 3*std, colorr, linestyle:, labelUpper control limit) plt.axhline(mean - 3*std, colorr, linestyle:, labelLower control limit) plt.xlabel(Day) plt.ylabel(Number of defective items) plt.title(Quality Control Chart) plt.legend() plt.show()实际应用建议当数据点超出控制限时可能表明生产过程出现了异常波动需要调查原因。