德国彩票Monte Carlo模拟：理论概率计算与代码优化咨询

问题背景 #

我正在编写Monte Carlo模拟程序，研究德国彩票中最冷门数字13在3570次抽奖中被抽取次数≤361的概率，已完成Python模拟代码。现需要：

1. 对应场景的理论概率计算公式
2. 针对random.sample的代码性能优化方案

原模拟代码 #

import random as rnd
import numpy as np

#Variablen setzen 
anzahl_sim=10000
lottozahlen_alle = list(range(1,50))
sim_wie_realität = 0
liste_niedrgiste_häufigkeit=[]

#Schleife für Simulation
for i in range(anzahl_sim):
    häufigkeit_zahlen_aus_lottoziehung = {zahl: 0 for zahl in lottozahlen_alle}
    for j in range(3570):           #sci.py.stats.hypergeom für matrixbasierte Simulatione
        x = rnd.sample(lottozahlen_alle, 6)     #rnd.sample funktioniert wie "Urnenziehung"
        for zahl in x:
                häufigkeit_zahlen_aus_lottoziehung[zahl] += 1
    summe_häufigkeit_bei_3570_ziehungen = sum(häufigkeit_zahlen_aus_lottoziehung.values())
    niedrigste_zahl = min(häufigkeit_zahlen_aus_lottoziehung, key=häufigkeit_zahlen_aus_lottoziehung.get)
    niedrigste_häufigkeit = häufigkeit_zahlen_aus_lottoziehung[niedrigste_zahl]
    if niedrigste_häufigkeit <= 361:      #Prüfen, ob Simulation von 3570 Lottoziehungen Ergebnis wie Realität aufweist
        sim_wie_realität += 1
    liste_niedrgiste_häufigkeit.append(niedrigste_häufigkeit)   #Sammlung von seltensten Kugeln bei 3570 Ziehungen

array_niedrgiste_häufigkeit=np.array(liste_niedrgiste_häufigkeit)

p_sim_wie_realität = sim_wie_realität/anzahl_sim
avg_experiment = array_niedrgiste_häufigkeit.sum(axis=0)/anzahl_sim
std_experiment = array_niedrgiste_häufigkeit.std(axis=0)
var_experiment = array_niedrgiste_häufigkeit.var(axis=0)
entfernung_in_std = abs((361-avg_experiment)/std_experiment)

print("Mittelwert von Anzahl niedrigste Ziehung: ", avg_experiment)
print("Standardabweichung: ", std_experiment)
print("Varianz: ", var_experiment)
print("Entfernung in Standardabweichung realtität von Mittelwert Simulation: ", entfernung_in_std)
print("P(eine Zahl ist seltener als 361 mal oder 361 mal gezogen worden) für ", anzahl_sim, "Simmulationen: ", p_sim_wie_realität)

一、理论概率计算公式 #

注意：原代码统计的是所有49个数字中至少有一个数字出现次数≤361的概率，而问题描述提到的是特定数字13出现次数≤361的概率，两者需区分：

1. 特定数字（如13）出现次数≤361的概率 #

德国彩票每次从1-49中抽取6个不重复数字，单次抽奖中数字13被抽到的概率为：
$$p = \frac{\binom{48}{5}}{\binom{49}{6}} = \frac{6}{49}$$

3570次抽奖中，数字13被抽到的次数$X$服从二项分布$Binomial(n=3570, p=\frac{6}{49})$，其累积概率为：
$$P(X \leq 361) = \sum_{k=0}^{361} \binom{3570}{k} \left(\frac{6}{49}\right)^k \left(\frac{43}{49}\right)^{3570-k}$$

由于$n$较大，可采用正态近似简化计算：

• 均值：$\mu = n \cdot p = 3570 \times \frac{6}{49} \approx 437.14$
• 方差：$\sigma^2 = n \cdot p \cdot (1-p) = 3570 \times \frac{6}{49} \times \frac{43}{49} \approx 383.61$
• 标准差：$\sigma \approx 19.59$

加入连续性修正后，计算Z值：
$$Z = \frac{361 + 0.5 - \mu}{\sigma} \approx \frac{361.5 - 437.14}{19.59} \approx -3.86$$

查标准正态分布表可得$P(Z \leq -3.86)$，即近似概率。

2. 所有数字中至少一个出现次数≤361的概率 #

由于49个数字的出现次数总和固定为$3570 \times 6 = 21420$，变量间存在相依性，理论计算复杂度极高。可通过补集思想近似：
$$P(\text{至少一个数字次数≤361}) = 1 - P(\text{所有数字次数>361})$$

但因变量非独立，精确计算需涉及多元超几何分布的联合概率，实际中Monte Carlo模拟是更可行的方案。

二、代码性能优化方案 #

原代码嵌套循环调用random.sample，效率极低（10000次模拟需执行3570万次random.sample），可通过numpy向量化操作大幅提升性能：

优化思路 #

1. 用numpy.random.choice代替random.sample，一次性生成所有模拟的抽奖结果，避免循环开销
2. 用np.bincount批量统计数字出现次数，替代字典计数的循环操作
3. 移除冗余计算（如原代码中summe_häufigkeit_bei_3570_ziehungen无实际作用）

优化后代码 #

import numpy as np

# 参数设置
anzahl_sim = 10000
total_draws = 3570
numbers = np.arange(1, 50)
draws_per_trial = 6

# 一次性生成所有模拟的抽奖结果：形状(模拟次数, 每次模拟抽奖次数, 单次抽奖数字数)
all_draws = np.random.choice(numbers, size=(anzahl_sim, total_draws, draws_per_trial), replace=False)

# 统计每个模拟中各数字的出现次数
# 展平每个模拟的抽奖结果，再用bincount统计计数
counts = np.array([np.bincount(draws.flatten(), minlength=50) for draws in all_draws])
# 提取数字1-49的计数，计算每个模拟的最低出现次数
min_counts = counts[:, 1:50].min(axis=1)

# 计算统计量
sim_wie_realität = (min_counts <= 361).sum()
p_sim_wie_realität = sim_wie_realität / anzahl_sim
avg_experiment = min_counts.mean()
std_experiment = min_counts.std()
var_experiment = min_counts.var()
entfernung_in_std = abs((361 - avg_experiment) / std_experiment)

# 输出结果
print(f"Mittelwert von Anzahl niedrigste Ziehung: {avg_experiment:.4f}")
print(f"Standardabweichung: {std_experiment:.4f}")
print(f"Varianz: {var_experiment:.4f}")
print(f"Entfernung in Standardabweichung realität von Mittelwert Simulation: {entfernung_in_std:.4f}")
print(f"P(eine Zahl ist seltener als 361 mal oder 361 mal gezogen worden) für {anzahl_sim} Simulationen: {p_sim_wie_realität:.4f}")

进一步优化（内存友好） #

若一次性生成所有数据内存不足，可分批次处理：

import numpy as np

anzahl_sim = 10000
batch_size = 1000
total_draws = 3570
numbers = np.arange(1, 50)
draws_per_trial = 6

min_counts = []
for _ in range(anzahl_sim // batch_size):
    # 生成批次内的抽奖结果
    batch_draws = np.random.choice(numbers, size=(batch_size, total_draws, draws_per_trial), replace=False)
    # 统计批次内的最低次数
    batch_counts = np.array([np.bincount(draws.flatten(), minlength=50) for draws in batch_draws])
    batch_min = batch_counts[:, 1:50].min(axis=1)
    min_counts.extend(batch_min)

min_counts = np.array(min_counts)
# 后续统计量计算同优化后代码

内容的提问来源于stack exchange，提问作者Clems