未记录事实核验 · 状态：现行内容
记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。
warning 适用边界与维护说明
本文适用于在给定参数和模型假设下学习二项分布计算、模拟、近似与统计推断；业务及医学参数均为教学设定，不构成实际效果证明。各次试验须独立且成功概率相同，有限批次不放回抽样、客户差异或患者异质性需另行判断模型适用性。正态与泊松近似条件是经验判断，尤其尾部概率应结合精确二项计算评估误差；营销案例的95%区间是给定成功率下成功次数的分布区间，而非成功率的置信区间。硬币案例的双尾概率写法限于当前对称零假设及观察值，不宜直接泛化到任意成功概率。医学部分仅示范历史成功率视为已知时的单侧检验，样本量结果限于50至490、步长10的搜索网格，不代表完整临床试验设计或临床有效性证明。本次仅作全文结构审阅，未运行代码或进行外部事实核验。
🎲 二项分布（Binomial Distribution）
🤔 什么是二项分布？
note 📖 **二项分布的定义**
二项分布描述了在n次独立的伯努利试验中，成功次数X的概率分布。每次试验只有两种可能结果：成功（概率p）或失败（概率1-p）。
记法：$X \sim B(n, p)$
概率质量函数： $$P(X = k) = \binom{n}{k} p^k (1-p)^{n-k}$$
其中：$\binom{n}{k} = \frac{n!}{k!(n-k)!}$
生活例子：
- 投掷10次硬币，正面出现的次数
- 100个产品中不合格品的数量
- 向20个客户推销，成功销售的数量
- 100次射击中命中靶心的次数
📊 二项分布的特征
1. 基本性质
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
from scipy.special import comb
def binomial_properties():
"""展示二项分布的基本性质"""
# 参数设定
n = 10 # 试验次数
p = 0.3 # 成功概率
print(f"二项分布 B({n}, {p}):")
print("=" * 25)
# 理论值
mean_theory = n * p
var_theory = n * p * (1 - p)
std_theory = np.sqrt(var_theory)
print(f"理论值:")
print(f" 均值 E(X) = np = {n} × {p} = {mean_theory:.2f}")
print(f" 方差 Var(X) = np(1-p) = {n} × {p} × {1-p} = {var_theory:.2f}")
print(f" 标准差 SD(X) = √Var(X) = {std_theory:.2f}")
# 使用scipy验证
rv = stats.binom(n, p)
mean_scipy = rv.mean()
var_scipy = rv.var()
std_scipy = rv.std()
print(f"\nSciPy验证:")
print(f" 均值: {mean_scipy:.2f}")
print(f" 方差: {var_scipy:.2f}")
print(f" 标准差: {std_scipy:.2f}")
# 概率计算示例
print(f"\n概率计算示例:")
for k in [0, 1, 3, 5, 10]:
prob_manual = comb(n, k) * (p**k) * ((1-p)**(n-k))
prob_scipy = rv.pmf(k)
print(f" P(X = {k}) = {prob_manual:.4f} (手算) = {prob_scipy:.4f} (scipy)")
return n, p, rv
n, p, rv = binomial_properties()
2. 形状特征
def plot_binomial_shapes():
"""绘制不同参数下的二项分布"""
fig, axes = plt.subplots(2, 2, figsize=(12, 10))
axes = axes.flatten()
# 不同参数组合
params = [
(10, 0.1, "n=10, p=0.1 (右偏)"),
(10, 0.5, "n=10, p=0.5 (对称)"),
(10, 0.8, "n=10, p=0.8 (左偏)"),
(30, 0.3, "n=30, p=0.3 (近似正态)")
]
for i, (n, p, title) in enumerate(params):
rv = stats.binom(n, p)
k_values = np.arange(0, n+1)
probabilities = rv.pmf(k_values)
# 绘制概率质量函数
axes[i].bar(k_values, probabilities, alpha=0.7,
color=f'C{i}', edgecolor='black', linewidth=0.5)
# 标注均值和标准差
mean_val = rv.mean()
std_val = rv.std()
axes[i].axvline(mean_val, color='red', linestyle='--',
label=f'μ={mean_val:.1f}')
axes[i].axvline(mean_val - std_val, color='orange', linestyle=':', alpha=0.7)
axes[i].axvline(mean_val + std_val, color='orange', linestyle=':', alpha=0.7)
axes[i].set_title(title)
axes[i].set_xlabel('成功次数 k')
axes[i].set_ylabel('概率 P(X=k)')
axes[i].legend()
axes[i].grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 形状分析
print("二项分布的形状特征:")
print("=" * 25)
print("• p = 0.5时，分布对称")
print("• p < 0.5时，分布右偏")
print("• p > 0.5时，分布左偏")
print("• n增大时，分布趋近正态分布")
plot_binomial_shapes()
💻 Python实现
1. 基本计算
def binomial_calculations():
"""二项分布的各种计算"""
n, p = 20, 0.3
rv = stats.binom(n, p)
print(f"二项分布 B({n}, {p}) 的计算:")
print("=" * 30)
# 1. 概率质量函数 (PMF)
print("1. 概率质量函数:")
for k in [0, 5, 6, 7, 10, 15, 20]:
prob = rv.pmf(k)
print(f" P(X = {k:2d}) = {prob:.4f}")
# 2. 累积分布函数 (CDF)
print("\n2. 累积分布函数:")
for k in [5, 6, 7, 10]:
cdf_val = rv.cdf(k)
print(f" P(X ≤ {k:2d}) = {cdf_val:.4f}")
# 3. 生存函数 (Survival Function)
print("\n3. 生存函数:")
for k in [5, 6, 7, 10]:
sf_val = rv.sf(k)
print(f" P(X > {k:2d}) = {sf_val:.4f}")
# 4. 分位数函数 (PPF)
print("\n4. 分位数:")
for percentile in [0.25, 0.5, 0.75, 0.9]:
quantile = rv.ppf(percentile)
print(f" 第{percentile*100:4.0f}百分位数: {quantile:.0f}")
# 5. 区间概率
print("\n5. 区间概率:")
intervals = [(3, 7), (5, 10), (0, 5)]
for a, b in intervals:
prob_interval = rv.cdf(b) - rv.cdf(a-1)
print(f" P({a} ≤ X ≤ {b}) = {prob_interval:.4f}")
return rv
rv = binomial_calculations()
2. 随机数生成和模拟
def binomial_simulation():
"""二项分布的模拟验证"""
n, p = 15, 0.4
n_simulations = 10000
print(f"二项分布模拟验证:")
print(f"参数: n={n}, p={p}")
print(f"模拟次数: {n_simulations:,}")
print("=" * 30)
# 方法1: 使用numpy生成随机数
np.random.seed(42)
samples_numpy = np.random.binomial(n, p, n_simulations)
# 方法2: 使用scipy生成随机数
np.random.seed(42)
rv = stats.binom(n, p)
samples_scipy = rv.rvs(size=n_simulations)
# 方法3: 手动模拟（通过多次伯努利试验）
np.random.seed(42)
samples_manual = []
for _ in range(n_simulations):
# 进行n次伯努利试验
successes = np.sum(np.random.random(n) < p)
samples_manual.append(successes)
samples_manual = np.array(samples_manual)
# 比较结果
methods = ['NumPy', 'SciPy', '手动模拟']
samples_list = [samples_numpy, samples_scipy, samples_manual]
# 理论值
theory_mean = n * p
theory_var = n * p * (1 - p)
theory_std = np.sqrt(theory_var)
print(f"理论值: 均值={theory_mean:.2f}, 标准差={theory_std:.2f}")
print("\n模拟结果比较:")
print("方法\t\t均值\t标准差\t最小值\t最大值")
print("-" * 50)
for method, samples in zip(methods, samples_list):
mean_sim = np.mean(samples)
std_sim = np.std(samples, ddof=1)
min_sim = np.min(samples)
max_sim = np.max(samples)
print(f"{method}\t\t{mean_sim:.2f}\t{std_sim:.2f}\t{min_sim}\t{max_sim}")
# 可视化比较
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(15, 6))
# 理论PMF vs 模拟频率
k_values = np.arange(0, n+1)
theory_pmf = rv.pmf(k_values)
# 计算模拟频率 (使用NumPy结果)
unique, counts = np.unique(samples_numpy, return_counts=True)
sim_freq = counts / n_simulations
ax1.bar(k_values - 0.2, theory_pmf, width=0.4, alpha=0.7,
label='理论PMF', color='blue')
# 为模拟频率创建完整的数组
sim_pmf = np.zeros(len(k_values))
for val, freq in zip(unique, sim_freq):
if val < len(sim_pmf):
sim_pmf[val] = freq
ax1.bar(k_values + 0.2, sim_pmf, width=0.4, alpha=0.7,
label='模拟频率', color='red')
ax1.set_xlabel('成功次数 k')
ax1.set_ylabel('概率/频率')
ax1.set_title('理论vs模拟分布')
ax1.legend()
ax1.grid(True, alpha=0.3)
# Q-Q图检验
from scipy import stats as scipy_stats
scipy_stats.probplot(samples_numpy, dist=stats.binom(n, p), plot=ax2)
ax2.set_title('Q-Q图: 模拟数据 vs 理论分布')
ax2.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
return samples_numpy
samples = binomial_simulation()
🎯 实际应用案例
案例1：质量控制
def quality_control_case():
"""产品质量控制案例"""
print("案例1: 产品质量控制")
print("=" * 25)
# 背景：某工厂生产的产品不合格率为5%
# 随机抽检50件产品，分析不合格品数量
n_products = 50 # 抽检数量
defect_rate = 0.05 # 不合格率5%
rv = stats.binom(n_products, defect_rate)
print(f"问题设定:")
print(f" 抽检产品数: {n_products}件")
print(f" 不合格率: {defect_rate*100}%")
print(f" 分布: X ~ B({n_products}, {defect_rate})")
# 期望值
expected_defects = rv.mean()
std_defects = rv.std()
print(f"\n统计特征:")
print(f" 期望不合格品数: {expected_defects:.2f}件")
print(f" 标准差: {std_defects:.2f}件")
# 关键概率计算
print(f"\n关键概率:")
# 1. 没有不合格品的概率
prob_zero = rv.pmf(0)
print(f" P(X = 0) = {prob_zero:.4f} ({prob_zero*100:.2f}%)")
# 2. 最多2件不合格的概率
prob_at_most_2 = rv.cdf(2)
print(f" P(X ≤ 2) = {prob_at_most_2:.4f} ({prob_at_most_2*100:.2f}%)")
# 3. 超过5件不合格的概率（异常情况）
prob_more_than_5 = rv.sf(5)
print(f" P(X > 5) = {prob_more_than_5:.4f} ({prob_more_than_5*100:.2f}%)")
# 4. 控制限的设定（3σ原则）
control_limit = expected_defects + 3 * std_defects
prob_exceed_limit = rv.sf(int(control_limit))
print(f"\n质量控制分析:")
print(f" 控制上限 (μ + 3σ): {control_limit:.1f}件")
print(f" 超过控制限概率: {prob_exceed_limit:.6f} ({prob_exceed_limit*10000:.2f}万分之一)")
# 决策建议
print(f"\n决策建议:")
if prob_more_than_5 < 0.01:
print(f" 如果发现超过5件不合格品(概率<1%)，应调查生产过程")
# 可视化
k_values = np.arange(0, 15)
probabilities = rv.pmf(k_values)
plt.figure(figsize=(10, 6))
plt.bar(k_values, probabilities, alpha=0.7, color='steelblue', edgecolor='black')
plt.axvline(expected_defects, color='red', linestyle='--', linewidth=2,
label=f'期望值 = {expected_defects:.2f}')
plt.axvline(control_limit, color='orange', linestyle='--', linewidth=2,
label=f'控制限 = {control_limit:.1f}')
plt.xlabel('不合格品数量')
plt.ylabel('概率')
plt.title(f'不合格品数量分布 B({n_products}, {defect_rate})')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
quality_control_case()
案例2：营销活动效果评估
def marketing_campaign_case():
"""营销活动效果评估案例"""
print("\n\n案例2: 营销活动效果评估")
print("=" * 30)
# 背景：电话营销活动，历史成功率20%
# 今天拨打了80个电话，分析成功销售数量
n_calls = 80 # 拨打电话数
success_rate = 0.20 # 成功率20%
rv = stats.binom(n_calls, success_rate)
print(f"营销活动设定:")
print(f" 拨打电话数: {n_calls}个")
print(f" 历史成功率: {success_rate*100}%")
print(f" 分布: X ~ B({n_calls}, {success_rate})")
# 统计特征
expected_success = rv.mean()
std_success = rv.std()
print(f"\n预期结果:")
print(f" 期望成功数: {expected_success:.1f}个")
print(f" 标准差: {std_success:.2f}个")
# 业务关键问题
print(f"\n业务分析:")
# 1. 至少完成15个销售的概率
target_sales = 15
prob_at_least_target = rv.sf(target_sales - 1)
print(f" P(X ≥ {target_sales}) = {prob_at_least_target:.4f} ({prob_at_least_target*100:.2f}%)")
# 2. 90%把握下的最低成功数
confidence_90 = rv.ppf(0.1) # 10%分位数
print(f" 90%把握下最低成功数: {confidence_90:.0f}个")
# 3. 95%区间估计
lower_95 = rv.ppf(0.025)
upper_95 = rv.ppf(0.975)
print(f" 95%区间: [{lower_95:.0f}, {upper_95:.0f}]个")
# 4. 异常高效果的判断
exceptional_threshold = expected_success + 2 * std_success
prob_exceptional = rv.sf(int(exceptional_threshold))
print(f" 超常发挥阈值 (μ + 2σ): {exceptional_threshold:.1f}个")
print(f" 超常发挥概率: {prob_exceptional:.4f} ({prob_exceptional*100:.2f}%)")
# ROI分析
cost_per_call = 5 # 每通电话成本5元
revenue_per_sale = 200 # 每次成功销售收入200元
total_cost = n_calls * cost_per_call
expected_revenue = expected_success * revenue_per_sale
expected_profit = expected_revenue - total_cost
expected_roi = (expected_profit / total_cost) * 100
print(f"\nROI分析:")
print(f" 总成本: {total_cost}元")
print(f" 期望收入: {expected_revenue:.0f}元")
print(f" 期望利润: {expected_profit:.0f}元")
print(f" 期望ROI: {expected_roi:.1f}%")
# 盈亏平衡分析
breakeven_sales = total_cost / revenue_per_sale
prob_breakeven = rv.sf(int(np.floor(breakeven_sales))) # 严格盈利：成功销售数大于盈亏平衡点
print(f" 盈亏平衡点: {breakeven_sales:.1f}个销售")
print(f" 盈利概率: {prob_breakeven:.4f} ({prob_breakeven*100:.2f}%)")
# 可视化结果
k_values = np.arange(0, 35)
probabilities = rv.pmf(k_values)
plt.figure(figsize=(12, 6))
plt.bar(k_values, probabilities, alpha=0.7, color='lightgreen', edgecolor='black')
# 标注关键点
plt.axvline(expected_success, color='blue', linestyle='-', linewidth=2,
label=f'期望值 = {expected_success:.1f}')
plt.axvline(target_sales, color='red', linestyle='--', linewidth=2,
label=f'目标 = {target_sales}')
plt.axvline(breakeven_sales, color='orange', linestyle=':', linewidth=2,
label=f'盈亏平衡 = {breakeven_sales:.1f}')
# 填充区域
k_success = k_values[k_values >= target_sales]
prob_success = probabilities[k_values >= target_sales]
plt.fill_between(k_success, 0, prob_success, alpha=0.3, color='red',
label=f'达到目标区域')
plt.xlabel('成功销售数量')
plt.ylabel('概率')
plt.title(f'营销活动成功数量分布 B({n_calls}, {success_rate})')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
marketing_campaign_case()
案例3：医学试验设计
def medical_trial_case():
"""医学试验设计案例"""
print("\n\n案例3: 医学试验设计")
print("=" * 20)
# 背景：新药治疗效果试验
# 已知标准治疗成功率60%，新药预期成功率75%
n_patients = 100 # 试验患者数
standard_rate = 0.60 # 标准治疗成功率
new_drug_rate = 0.75 # 新药预期成功率
rv_standard = stats.binom(n_patients, standard_rate)
rv_new_drug = stats.binom(n_patients, new_drug_rate)
print(f"试验设计:")
print(f" 患者数量: {n_patients}人")
print(f" 标准治疗成功率: {standard_rate*100}%")
print(f" 新药预期成功率: {new_drug_rate*100}%")
# 统计特征比较
print(f"\n预期结果比较:")
print(f" 标准治疗: 期望{rv_standard.mean():.1f}人成功, 标准差{rv_standard.std():.2f}")
print(f" 新药治疗: 期望{rv_new_drug.mean():.1f}人成功, 标准差{rv_new_drug.std():.2f}")
# 试验成功标准设定
# 假设患者结局独立、成功概率相同且历史基准0.60视为已知；成功人数超过基准分布的95%分位数时，拒绝H₀: p=0.60，支持H₁: p>0.60，但不据此单独认定临床有效性
threshold = rv_standard.ppf(0.95)
prob_new_drug_success = rv_new_drug.sf(int(threshold))
print(f"\n试验成功标准:")
print(f" 成功阈值: >{threshold:.0f}人")
print(f" 新药达到标准的概率: {prob_new_drug_success:.4f} ({prob_new_drug_success*100:.2f}%)")
# 检验力分析
# 在α=0.05下，新药真实有效时被检出的概率
alpha = 0.05
critical_value = rv_standard.ppf(1 - alpha)
power = rv_new_drug.sf(int(critical_value))
print(f"\n检验力分析:")
print(f" 显著性水平 α = {alpha}")
print(f" 临界值: {critical_value:.0f}人")
print(f" 检验力 (1-β): {power:.4f} ({power*100:.2f}%)")
# 样本量计算（为达到80%检验力所需的样本量）
def calculate_sample_size_for_power(target_power=0.8):
"""计算达到目标检验力所需的样本量"""
for n in range(50, 500, 10):
rv_std_n = stats.binom(n, standard_rate)
rv_new_n = stats.binom(n, new_drug_rate)
crit_val = rv_std_n.ppf(1 - alpha)
power_n = rv_new_n.sf(int(crit_val))
if power_n >= target_power:
return n, power_n
return None, None
optimal_n, optimal_power = calculate_sample_size_for_power()
print(f"\n样本量计算:")
print(f" 在50至490、步长10的搜索网格中首次达到80%检验力的样本量: {optimal_n}人")
print(f" 对应检验力: {optimal_power:.4f} ({optimal_power*100:.2f}%)")
# 风险分析
# 第一类错误：H₀: p=standard_rate成立时，成功人数仍落入拒绝域
type1_error = rv_standard.sf(int(critical_value))
# 第二类错误：新药实际有效，但未能检出差异
type2_error = 1 - power
print(f"\n风险分析:")
print(f" 实际第一类错误率（名义α={alpha}）: {type1_error:.3f} ({type1_error*100:.1f}%)")
print(f" 第二类错误 (β): {type2_error:.3f} ({type2_error*100:.1f}%)")
# 可视化两个分布
k_values = np.arange(40, 90)
prob_standard = rv_standard.pmf(k_values)
prob_new_drug = rv_new_drug.pmf(k_values)
plt.figure(figsize=(12, 6))
plt.bar(k_values - 0.2, prob_standard, width=0.4, alpha=0.7,
label=f'标准治疗 B({n_patients}, {standard_rate})', color='blue')
plt.bar(k_values + 0.2, prob_new_drug, width=0.4, alpha=0.7,
label=f'新药治疗 B({n_patients}, {new_drug_rate})', color='red')
plt.axvline(critical_value, color='black', linestyle='--', linewidth=2,
label=f'临界值 = {critical_value:.0f}')
# 填充检验力区域
k_power = k_values[k_values > critical_value]
prob_power = rv_new_drug.pmf(k_power)
plt.fill_between(k_power + 0.2, 0, prob_power, alpha=0.3, color='green',
label=f'检验力区域')
plt.xlabel('成功治疗人数')
plt.ylabel('概率')
plt.title('医学试验：标准治疗 vs 新药治疗')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
medical_trial_case()
🔄 二项分布的近似
1. 正态近似
def normal_approximation():
"""二项分布的正态近似"""
print("\n\n二项分布的正态近似:")
print("=" * 25)
# 大样本情况
n = 100
p = 0.3
rv_binom = stats.binom(n, p)
# 正态近似参数
mu = n * p
sigma = np.sqrt(n * p * (1 - p))
rv_normal = stats.norm(mu, sigma)
print(f"原分布: B({n}, {p})")
print(f"正态近似: N({mu:.1f}, {sigma:.2f}²)")
# 近似条件检查
print(f"\n近似条件检查:")
print(f" np = {n*p:.1f} ≥ 5? {'是' if n*p >= 5 else '否'}")
print(f" n(1-p) = {n*(1-p):.1f} ≥ 5? {'是' if n*(1-p) >= 5 else '否'}")
# 比较概率计算
test_values = [25, 30, 35, 40]
print(f"\n概率比较:")
print("k\t二项分布\t正态近似\t连续性校正\t误差")
print("-" * 60)
for k in test_values:
# 精确二项概率
prob_exact = rv_binom.pmf(k)
# 正态近似（不校正）
prob_normal = rv_normal.pdf(k)
# 连续性校正
prob_corrected = rv_normal.cdf(k + 0.5) - rv_normal.cdf(k - 0.5)
error = abs(prob_exact - prob_corrected)
print(f"{k}\t{prob_exact:.6f}\t{prob_normal:.6f}\t{prob_corrected:.6f}\t{error:.6f}")
# 累积概率比较
print(f"\n累积概率比较 P(X ≤ k):")
print("k\t二项分布\t正态近似\t连续性校正\t误差")
print("-" * 60)
for k in test_values:
# 精确累积概率
cdf_exact = rv_binom.cdf(k)
# 正态近似（不校正）
cdf_normal = rv_normal.cdf(k)
# 连续性校正
cdf_corrected = rv_normal.cdf(k + 0.5)
error = abs(cdf_exact - cdf_corrected)
print(f"{k}\t{cdf_exact:.6f}\t{cdf_normal:.6f}\t{cdf_corrected:.6f}\t{error:.6f}")
# 可视化比较
k_values = np.arange(15, 50)
prob_binom = rv_binom.pmf(k_values)
# 正态密度（连续）
x_continuous = np.linspace(15, 50, 200)
density_normal = rv_normal.pdf(x_continuous)
# 正态近似（离散化）
prob_normal_discrete = rv_normal.pdf(k_values)
plt.figure(figsize=(12, 6))
# 二项分布
plt.bar(k_values, prob_binom, alpha=0.7, color='blue',
label=f'二项分布 B({n}, {p})', edgecolor='black', linewidth=0.5)
# 正态近似（连续）
plt.plot(x_continuous, density_normal, 'r-', linewidth=2,
label=f'正态近似 N({mu:.1f}, {sigma:.2f}²)')
# 正态近似（离散点）
plt.scatter(k_values, prob_normal_discrete, color='red', s=20, alpha=0.7)
plt.xlabel('k')
plt.ylabel('概率/密度')
plt.title('二项分布的正态近似')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
print(f"\n💡 使用建议:")
print(f"- 当np ≥ 5且n(1-p) ≥ 5时，正态近似效果较好")
print(f"- 计算累积概率时使用连续性校正：P(X ≤ k) ≈ Φ((k+0.5-μ)/σ)")
print(f"- p接近0.5时近似效果最好")
normal_approximation()
2. 泊松近似
def poisson_approximation():
"""二项分布的泊松近似"""
print("\n\n二项分布的泊松近似:")
print("=" * 25)
# 大n小p的情况
n = 1000
p = 0.003
lambda_param = n * p
rv_binom = stats.binom(n, p)
rv_poisson = stats.poisson(lambda_param)
print(f"原分布: B({n}, {p})")
print(f"泊松近似: Poisson({lambda_param})")
# 近似条件检查
print(f"\n近似条件检查:")
print(f" n = {n} (大)")
print(f" p = {p} (小)")
print(f" λ = np = {lambda_param} (适中)")
# 比较概率
test_values = [0, 1, 2, 3, 4, 5, 6, 7, 8]
print(f"\n概率比较 P(X = k):")
print("k\t二项分布\t泊松近似\t误差")
print("-" * 40)
for k in test_values:
prob_binom = rv_binom.pmf(k)
prob_poisson = rv_poisson.pmf(k)
error = abs(prob_binom - prob_poisson)
print(f"{k}\t{prob_binom:.6f}\t{prob_poisson:.6f}\t{error:.6f}")
# 累积概率比较
print(f"\n累积概率比较 P(X ≤ k):")
print("k\t二项分布\t泊松近似\t误差")
print("-" * 40)
for k in test_values:
cdf_binom = rv_binom.cdf(k)
cdf_poisson = rv_poisson.cdf(k)
error = abs(cdf_binom - cdf_poisson)
print(f"{k}\t{cdf_binom:.6f}\t{cdf_poisson:.6f}\t{error:.6f}")
# 可视化比较
k_range = np.arange(0, 15)
prob_binom = rv_binom.pmf(k_range)
prob_poisson = rv_poisson.pmf(k_range)
plt.figure(figsize=(12, 6))
width = 0.35
plt.bar(k_range - width/2, prob_binom, width, alpha=0.7,
label=f'二项分布 B({n}, {p})', color='blue')
plt.bar(k_range + width/2, prob_poisson, width, alpha=0.7,
label=f'泊松近似 Poisson({lambda_param})', color='red')
plt.xlabel('k')
plt.ylabel('概率')
plt.title('二项分布的泊松近似 (大n小p情况)')
plt.legend()
plt.grid(True, alpha=0.3)
plt.show()
print(f"\n💡 使用条件:")
print(f"- n很大 (通常n ≥ 100)")
print(f"- p很小 (通常p ≤ 0.01)")
print(f"- np适中 (通常np < 10)")
print(f"- 经验法则: n ≥ 100, p ≤ 0.01, np ≤ 10")
poisson_approximation()
🛠️ 实用技巧和最佳实践
1. 参数估计
def parameter_estimation():
"""二项分布参数估计"""
print("\n\n二项分布参数估计:")
print("=" * 20)
# 模拟真实情况：已知试验次数n，估计成功概率p
true_n = 50
true_p = 0.35
# 生成观察数据
np.random.seed(42)
observed_successes = np.random.binomial(true_n, true_p, 100) # 100次实验
print(f"真实参数: n = {true_n}, p = {true_p}")
print(f"观察数据: 100次实验，每次{true_n}次试验")
# 最大似然估计
sample_mean = np.mean(observed_successes)
p_hat = sample_mean / true_n # MLE估计
print(f"\n参数估计:")
print(f" 样本均值: {sample_mean:.2f}")
print(f" p的MLE估计: {p_hat:.4f}")
print(f" 估计误差: {abs(p_hat - true_p):.4f}")
# 置信区间（正态近似）
total_trials = true_n * len(observed_successes)
se_p = np.sqrt(p_hat * (1 - p_hat) / total_trials) # 全部独立试验下的标准误
confidence_level = 0.95
z_alpha = stats.norm.ppf(1 - (1 - confidence_level)/2)
ci_lower = p_hat - z_alpha * se_p
ci_upper = p_hat + z_alpha * se_p
print(f"\n95%置信区间:")
print(f" [{ci_lower:.4f}, {ci_upper:.4f}]")
print(f" 真实值在区间内: {'是' if ci_lower <= true_p <= ci_upper else '否'}")
# Wilson置信区间（更准确）
def wilson_confidence_interval(x, n, confidence=0.95):
"""Wilson置信区间"""
z = stats.norm.ppf(1 - (1 - confidence)/2)
p_tilde = (x + z**2/2) / (n + z**2)
margin = z / (1 + z**2 / n) * np.sqrt((x / n) * (1 - x / n) / n + z**2 / (4 * n**2))
return p_tilde - margin, p_tilde + margin
wilson_lower, wilson_upper = wilson_confidence_interval(np.sum(observed_successes), total_trials)
print(f"\nWilson 95%置信区间:")
print(f" [{wilson_lower:.4f}, {wilson_upper:.4f}]")
print(f" 真实值在区间内: {'是' if wilson_lower <= true_p <= wilson_upper else '否'}")
parameter_estimation()
2. 假设检验应用
def binomial_hypothesis_test():
"""二项分布假设检验"""
print("\n\n二项分布假设检验:")
print("=" * 20)
# 案例：检验硬币是否公平
n_tosses = 100 # 投掷次数
observed_heads = 58 # 观察到的正面次数
null_p = 0.5 # 原假设：硬币公平
print(f"案例: 硬币公平性检验")
print(f" 投掷次数: {n_tosses}")
print(f" 观察正面: {observed_heads}次")
print(f" H₀: p = {null_p} (硬币公平)")
print(f" H₁: p ≠ {null_p} (硬币不公平)")
# 方法1: 精确二项检验
p_value_exact = 2 * min(stats.binom.cdf(observed_heads, n_tosses, null_p),
stats.binom.sf(observed_heads - 1, n_tosses, null_p))
print(f"\n精确二项检验:")
print(f" p值: {p_value_exact:.4f}")
# 方法2: 正态近似检验（连续性校正）
expected = n_tosses * null_p
variance = n_tosses * null_p * (1 - null_p)
std_error = np.sqrt(variance)
# 连续性校正
if observed_heads > expected:
z_stat = (observed_heads - 0.5 - expected) / std_error
else:
z_stat = (observed_heads + 0.5 - expected) / std_error
p_value_normal = 2 * (1 - stats.norm.cdf(abs(z_stat)))
print(f"\n正态近似检验 (连续性校正):")
print(f" 期望值: {expected}")
print(f" 标准误: {std_error:.2f}")
print(f" z统计量: {z_stat:.3f}")
print(f" p值: {p_value_normal:.4f}")
# 决策
alpha = 0.05
print(f"\n检验结论 (α = {alpha}):")
if p_value_exact < alpha:
print(f" 拒绝原假设 (p = {p_value_exact:.4f} < {alpha})")
print(f" 结论: 在独立同分布投掷假设下，有统计证据支持硬币正面概率不等于0.5")
else:
print(f" 不拒绝原假设 (p = {p_value_exact:.4f} ≥ {alpha})")
print(f" 结论: 没有足够证据表明硬币不公平")
# 置信区间
from statsmodels.stats.proportion import proportion_confint
ci_lower, ci_upper = proportion_confint(observed_heads, n_tosses, alpha=0.05)
print(f"\n观察比例的95%置信区间:")
print(f" 观察比例: {observed_heads/n_tosses:.3f}")
print(f" 置信区间: [{ci_lower:.3f}, {ci_upper:.3f}]")
print(f" 0.5在区间内: {'是' if ci_lower <= 0.5 <= ci_upper else '否'}")
binomial_hypothesis_test()
🔗 相关知识
#统计学 #二项分布 #概率分布 #离散分布 #伯努利试验 #基础