利用scipy distributions模拟指定双硬币抛掷实验及结果分布表征

先理清楚实验逻辑 #

我们有两枚硬币：

• 公平硬币（记为Coin A）：正面（H）概率0.5，反面（T）概率0.5
• 偏态硬币（记为Coin B）：反面（T）概率0.3，因此正面（H）概率为0.7

每次"双硬币抛掷"同时抛A和B各一次，会得到4种组合：HH、HT、TH、TT。我们需要重复这个双抛实验10次为一组，生成100组这样的独立样本，再对四类组合的分布做分析。

用Scipy实现模拟 #

Scipy的scipy.stats.bernoulli可以高效模拟单次伯努利试验（硬币抛掷），结合numpy批量处理样本：

import numpy as np
from scipy.stats import bernoulli
import matplotlib.pyplot as plt

# 定义硬币参数
p_coinA_H = 0.5  # 公平硬币正面概率
p_coinB_H = 0.7  # 偏态硬币正面概率（反面概率0.3）

# 生成100个样本，每个样本包含10次双硬币抛掷
# 数组形状：(样本数, 实验次数, 2)，最后一维对应两枚硬币的结果（1=正面，0=反面）
samples = np.zeros((100, 10, 2), dtype=int)
# 模拟Coin A的所有结果
samples[:, :, 0] = bernoulli.rvs(p=p_coinA_H, size=(100, 10))
# 模拟Coin B的所有结果
samples[:, :, 1] = bernoulli.rvs(p=p_coinB_H, size=(100, 10))

# 将每次双抛的数字结果转换为组合标签
def map_to_combination(trial):
    a, b = trial
    if a == 1 and b == 1:
        return "HH"
    elif a == 1 and b == 0:
        return "HT"
    elif a == 0 and b == 1:
        return "TH"
    else:
        return "TT"

# 为每个样本（10次实验）生成组合序列
sample_combinations = [list(map(map_to_combination, sample)) for sample in samples]

四类组合的分布表征 #

可以从模拟结果的频率分布和理论分布验证两个维度来表征：

1. 样本内组合次数的频率分布

统计每个10次实验样本中，四类组合各自出现的次数，再用直方图展示100个样本的次数分布：

# 统计每个样本中四类组合的出现次数
counts_per_sample = []
for combos in sample_combinations:
    counts = {
        "HH": combos.count("HH"),
        "HT": combos.count("HT"),
        "TH": combos.count("TH"),
        "TT": combos.count("TT")
    }
    counts_per_sample.append(counts)

# 提取每类组合的次数数组
hh_counts = np.array([c["HH"] for c in counts_per_sample])
ht_counts = np.array([c["HT"] for c in counts_per_sample])
th_counts = np.array([c["TH"] for c in counts_per_sample])
tt_counts = np.array([c["TT"] for c in counts_per_sample])

# 绘制分布直方图
fig, axes = plt.subplots(2, 2, figsize=(12, 8))
plot_params = {"bins": np.arange(0, 11)-0.5, "edgecolor": "black"}

axes[0,0].hist(hh_counts, **plot_params)
axes[0,0].set_title("HH出现次数分布")
axes[0,0].set_xlabel("次数")
axes[0,0].set_ylabel("样本数")

axes[0,1].hist(ht_counts, **plot_params)
axes[0,1].set_title("HT出现次数分布")
axes[0,1].set_xlabel("次数")

axes[1,0].hist(th_counts, **plot_params)
axes[1,0].set_title("TH出现次数分布")
axes[1,0].set_xlabel("次数")
axes[1,0].set_ylabel("样本数")

axes[1,1].hist(tt_counts, **plot_params)
axes[1,1].set_title("TT出现次数分布")
axes[1,1].set_xlabel("次数")

plt.tight_layout()
plt.show()

2. 理论分布验证

从概率学角度，每次双抛实验中四类组合的概率为：

• ( P(HH) = 0.5 \times 0.7 = 0.35 )
• ( P(HT) = 0.5 \times 0.3 = 0.15 )
• ( P(TH) = 0.5 \times 0.7 = 0.35 )
• ( P(TT) = 0.5 \times 0.3 = 0.15 )

每个样本是10次独立实验，因此每类组合的出现次数服从二项分布：

• HH次数 ~ ( \text{Binomial}(n=10, p=0.35) )
• HT次数 ~ ( \text{Binomial}(n=10, p=0.15) )
• TH次数与HH一致，TT次数与HT一致

可以用Scipy的binom生成理论分布，和模拟结果对比：

from scipy.stats import binom

n_trials = 10
x = np.arange(0, 11)

# 生成HH的理论概率质量函数
pmf_hh = binom.pmf(x, n_trials, 0.35)
# 生成HT的理论概率质量函数
pmf_ht = binom.pmf(x, n_trials, 0.15)

# 在直方图上叠加理论曲线（以HH为例）
axes[0,0].plot(x, pmf_hh * 100, "r-", linewidth=2, label="理论分布")
axes[0,0].legend()

关键说明 #

• Scipy的伯努利分布工具非常适合这类独立重复试验的模拟，批量生成样本效率很高
• 模拟结果的分布应该和理论二项分布高度吻合，样本量越大（比如增加到1000个样本），拟合效果越好
• 四类组合中，HH和TH的理论概率相同，HT和TT的理论概率相同，因此它们的分布形态也会一致

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