未记录事实核验 · 状态：维护中
记录仅说明核验时间与维护状态，不保证适用于所有环境；使用前请留意正文中的版本、适用条件与风险提示。
warning 适用边界与维护说明
本文适合学习泊松计数概率、参数估计、近似比较及业务模拟，不宜直接作为排班、容量或质量管理的已验证方案；本次仅作全文结构审阅，未执行代码或开展外部事实核验。拟合检验仍删除低期望频数类别并遗漏观测最大值之外的理论尾部，且对预先固定的λ错误扣除一个估计参数的自由度；默认λ=3、10小时样本没有类别达到期望频数5，因而不会进入检验输出分支，不能据此获得有效拟合结论。网站检测采用全天均值阈值，未处理模拟数据的昼夜变化，所谓“异常点”只能理解为混合阈值下的超限点。座席计算仅是固定小时处理能力下的工作量估算，另设服务率的单座席M/M/1模型不能验证该排班；质量比较和控制图要求批次检验规模或缺陷机会量可比，评级、预警规则及改进幅度均为示例设定。二项到泊松近似的数值条件仅作经验参考，不保证指定概率区域的精度；截断过滤后的KL项部分和不是完整KL散度。
📈 泊松分布（Poisson Distribution）
🤔 什么是泊松分布？
note 📖 **泊松分布的定义**
泊松分布是一种描述非负整数计数的概率分布，常用于固定时间或空间内的事件次数建模，但并不限于稀有事件。若二项分布的n趋于无穷、p趋于0且np趋于固定的λ，则其收敛于Poisson(λ)。
数学表达式： $$P(X = k) = \frac{\lambda^k e^{-\lambda}}{k!}, \quad k = 0, 1, 2, \ldots$$
记号：$X \sim \text{Poisson}(\lambda)$
其中 $\lambda > 0$ 是参数，表示所讨论观察区间内事件发生的期望次数；若齐次泊松过程的事件率为 $r$，观察时长为 $t$，则该区间计数的参数为 $\lambda=rt$。
生活例子：
- 每小时到达银行的顾客数量
- 每天收到的垃圾邮件数量
- 每年某地区发生的地震次数
- 网站每分钟的访问次数
- 电话交换机每小时接到的电话数
📐 泊松分布的性质
1. 基本特征
import numpy as np
import matplotlib.pyplot as plt
from scipy import stats
import seaborn as sns
def plot_poisson_properties():
"""展示泊松分布的基本性质"""
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 1. 不同λ值的泊松分布
lambda_values = [1, 3, 5, 10]
colors = ['red', 'blue', 'green', 'orange']
x = np.arange(0, 20)
for lam, color in zip(lambda_values, colors):
pmf = stats.poisson.pmf(x, lam)
ax1.plot(x, pmf, 'o-', color=color, linewidth=2, markersize=6,
label=f'λ = {lam}')
ax1.set_title('不同λ值的泊松分布')
ax1.set_xlabel('k (事件发生次数)')
ax1.set_ylabel('P(X = k)')
ax1.legend()
ax1.grid(True, alpha=0.3)
# 2. 泊松分布的对称性变化
lambda_sym = [1, 5, 15]
colors_sym = ['red', 'blue', 'green']
for lam, color in zip(lambda_sym, colors_sym):
x_range = np.arange(0, max(25, lam + 4*np.sqrt(lam)))
pmf = stats.poisson.pmf(x_range, lam)
ax2.plot(x_range, pmf, 'o-', color=color, linewidth=2, markersize=4,
label=f'λ = {lam}')
ax2.set_title('λ增大时分布形状的变化')
ax2.set_xlabel('k')
ax2.set_ylabel('P(X = k)')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 3. 累积分布函数
lam = 5
x_cdf = np.arange(0, 15)
cdf = stats.poisson.cdf(x_cdf, lam)
pmf = stats.poisson.pmf(x_cdf, lam)
ax3.bar(x_cdf, pmf, alpha=0.6, color='skyblue', label='PMF')
ax3.step(x_cdf, cdf, where='post', color='red', linewidth=2, label='CDF')
ax3.set_title(f'泊松分布的PMF和CDF (λ = {lam})')
ax3.set_xlabel('k')
ax3.set_ylabel('概率')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 4. 泊松过程的可加性
lam1, lam2 = 3, 2
lam_sum = lam1 + lam2
x_add = np.arange(0, 15)
pmf1 = stats.poisson.pmf(x_add, lam1)
pmf2 = stats.poisson.pmf(x_add, lam2)
pmf_sum = stats.poisson.pmf(x_add, lam_sum)
ax4.plot(x_add, pmf1, 'o-', color='red', label=f'Poisson({lam1})')
ax4.plot(x_add, pmf2, 's-', color='blue', label=f'Poisson({lam2})')
ax4.plot(x_add, pmf_sum, '^-', color='green', linewidth=2,
label=f'Poisson({lam_sum}) = Sum')
ax4.set_title('泊松分布的可加性')
ax4.set_xlabel('k')
ax4.set_ylabel('P(X = k)')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
plot_poisson_properties()
# 泊松分布的重要性质
print("泊松分布的重要性质:")
print("=" * 40)
print("1. 离散性: 只能取非负整数值")
print("2. 单参数: 仅由参数λ确定")
print("3. 均值 = 方差 = λ")
print("4. 可加性: 若X₁与X₂相互独立，且X₁ ~ Poisson(λ₁)、X₂ ~ Poisson(λ₂)，则X₁+X₂ ~ Poisson(λ₁+λ₂)")
print("5. 渐近正态性: λ很大时近似正态分布")
# 验证均值和方差
lambda_test = 8
theoretical_mean = lambda_test
theoretical_var = lambda_test
print(f"\n验证 λ = {lambda_test} 时:")
print(f"理论均值: {theoretical_mean}")
print(f"理论方差: {theoretical_var}")
# 通过模拟验证
sample_data = np.random.poisson(lambda_test, 10000)
sample_mean = np.mean(sample_data)
sample_var = np.var(sample_data, ddof=1)
print(f"模拟均值: {sample_mean:.3f}")
print(f"模拟方差: {sample_var:.3f}")
2. 泊松过程
def demonstrate_poisson_process():
"""演示泊松过程"""
np.random.seed(42)
# 参数设置
lambda_rate = 3 # 每小时平均事件发生次数
time_period = 10 # 观察10小时
print("泊松过程演示:")
print("=" * 40)
print(f"参数λ = {lambda_rate} (每小时平均{lambda_rate}次事件)")
print(f"观察时间: {time_period}小时")
# 1. 生成事件发生时间
# 泊松过程中，事件间隔时间服从指数分布
inter_arrival_times = np.random.exponential(1/lambda_rate, 100)
arrival_times = np.cumsum(inter_arrival_times)
# 只取观察期内的事件
arrival_times = arrival_times[arrival_times <= time_period]
n_events = len(arrival_times)
# 2. 计算每个时间单位内的事件数
hourly_counts = []
for hour in range(time_period):
count = np.sum((arrival_times >= hour) & (arrival_times < hour + 1))
hourly_counts.append(count)
# 3. 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 事件时间线
ax1.scatter(arrival_times, np.ones(len(arrival_times)), alpha=0.7, s=50)
ax1.set_xlim(0, time_period)
ax1.set_ylim(0.5, 1.5)
ax1.set_xlabel('时间 (小时)')
ax1.set_title(f'泊松过程事件时间线 (总共{n_events}个事件)')
ax1.grid(True, alpha=0.3)
# 每小时事件计数
hours = np.arange(time_period)
ax2.bar(hours, hourly_counts, alpha=0.7, color='skyblue', edgecolor='black')
ax2.axhline(y=lambda_rate, color='red', linestyle='--',
label=f'理论均值 λ = {lambda_rate}')
ax2.set_xlabel('时间 (小时)')
ax2.set_ylabel('事件数量')
ax2.set_title('每小时事件计数')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 事件计数分布
unique_counts, count_freq = np.unique(hourly_counts, return_counts=True)
count_prob = count_freq / len(hourly_counts)
# 理论泊松分布
x_theory = np.arange(0, max(hourly_counts) + 3)
theory_prob = stats.poisson.pmf(x_theory, lambda_rate)
ax3.bar(unique_counts, count_prob, alpha=0.7, color='lightgreen',
width=0.4, label='观测分布')
ax3.plot(x_theory, theory_prob, 'ro-', linewidth=2, markersize=6,
label=f'理论Poisson({lambda_rate})')
ax3.set_xlabel('每小时事件数')
ax3.set_ylabel('概率')
ax3.set_title('事件计数分布对比')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 累积事件数
cumulative_counts = np.arange(1, len(arrival_times) + 1)
ax4.step(arrival_times, cumulative_counts, where='post', linewidth=2)
# 理论期望累积
time_points = np.linspace(0, time_period, 100)
expected_cumulative = lambda_rate * time_points
ax4.plot(time_points, expected_cumulative, 'r--', linewidth=2,
label=f'理论期望 λt = {lambda_rate}t')
ax4.set_xlabel('时间 (小时)')
ax4.set_ylabel('累积事件数')
ax4.set_title('累积事件数 vs 时间')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 统计分析
print(f"\n统计分析结果:")
print(f"总事件数: {n_events}")
print(f"平均每小时事件数: {np.mean(hourly_counts):.2f}")
print(f"方差: {np.var(hourly_counts, ddof=1):.2f}")
print(f"理论均值 = 方差 = {lambda_rate}")
# 拟合优度检验
observed_freq = np.bincount(hourly_counts, minlength=max(hourly_counts)+1)
expected_freq = stats.poisson.pmf(np.arange(len(observed_freq)), lambda_rate) * time_period
# 只考虑期望频数 >= 5 的类别
valid_categories = expected_freq >= 5
if np.sum(valid_categories) >= 2:
chi2_stat = np.sum((observed_freq[valid_categories] - expected_freq[valid_categories])**2 /
expected_freq[valid_categories])
df = np.sum(valid_categories) - 1 - 1 # -1 for estimated parameter
p_value = 1 - stats.chi2.cdf(chi2_stat, df) if df > 0 else None
print(f"\n拟合优度检验:")
print(f"卡方统计量: {chi2_stat:.4f}")
if p_value is not None:
print(f"p值: {p_value:.4f}")
if p_value > 0.05:
print("结论: 若检验分类及前提有效，在0.05显著性水平下未拒绝指定泊松分布假设；这不证明数据服从泊松分布")
else:
print("结论: 若检验分类及前提有效，在0.05显著性水平下拒绝指定泊松分布假设")
demonstrate_poisson_process()
💻 Python实现
1. 基本计算
from scipy import stats
import numpy as np
def poisson_calculations():
"""泊松分布的各种计算"""
lam = 4 # λ = 4
print(f"泊松分布 Poisson({lam}) 的计算:")
print("=" * 40)
# 1. 概率质量函数值
k_values = [0, 1, 2, 3, 4, 5, 8, 10]
print("概率质量函数值:")
for k in k_values:
pmf_value = stats.poisson.pmf(k, lam)
print(f"P(X = {k}) = {pmf_value:.4f}")
# 2. 累积分布函数值
print(f"\n累积分布函数值:")
for k in [2, 4, 6, 8]:
cdf_value = stats.poisson.cdf(k, lam)
print(f"P(X ≤ {k}) = {cdf_value:.4f}")
# 3. 分位数
print(f"\n分位数:")
percentiles = [0.1, 0.25, 0.5, 0.75, 0.9, 0.95]
for p in percentiles:
quantile = stats.poisson.ppf(p, lam)
print(f"{p:.0%}分位数: {quantile}")
# 4. 区间概率
print(f"\n区间概率:")
intervals = [(2, 6), (1, 5), (3, 8)]
for a, b in intervals:
prob = stats.poisson.cdf(b, lam) - stats.poisson.cdf(a-1, lam)
print(f"P({a} ≤ X ≤ {b}) = {prob:.4f}")
# 5. 特殊概率
print(f"\n特殊概率:")
# 至少发生一次
prob_at_least_one = 1 - stats.poisson.pmf(0, lam)
print(f"P(X ≥ 1) = {prob_at_least_one:.4f}")
# 超过平均值
prob_above_mean = 1 - stats.poisson.cdf(lam, lam)
print(f"P(X > λ) = P(X > {lam}) = {prob_above_mean:.4f}")
# 在均值±标准差范围内
std_dev = np.sqrt(lam)
lower = max(0, int(lam - std_dev))
upper = int(lam + std_dev)
prob_within_1std = stats.poisson.cdf(upper, lam) - stats.poisson.cdf(lower-1, lam)
print(f"P({lower} ≤ X ≤ {upper}) = {prob_within_1std:.4f}")
poisson_calculations()
2. 参数估计
def poisson_parameter_estimation():
"""泊松分布的参数估计"""
# 真实参数
true_lambda = 6.5
n_samples = 1000
np.random.seed(42)
sample_data = np.random.poisson(true_lambda, n_samples)
print("泊松分布参数估计:")
print("=" * 40)
print(f"真实参数 λ = {true_lambda}")
print(f"样本量 n = {n_samples}")
# 1. 矩估计法 (方法1)
sample_mean = np.mean(sample_data)
lambda_moment = sample_mean # 泊松分布的矩估计就是样本均值
print(f"\n矩估计法:")
print(f"样本均值: {sample_mean:.4f}")
print(f"λ的矩估计: {lambda_moment:.4f}")
# 2. 最大似然估计 (方法2)
# 对于泊松分布，MLE也是样本均值
lambda_mle = sample_mean
print(f"\n最大似然估计:")
print(f"λ的MLE: {lambda_mle:.4f}")
# 3. 置信区间
# 使用正态近似 (当λ较大时)
se = np.sqrt(lambda_mle / n_samples) # 标准误差
z_975 = stats.norm.ppf(0.975) # 97.5%分位数
ci_lower = lambda_mle - z_975 * se
ci_upper = lambda_mle + z_975 * se
print(f"\n95%置信区间:")
print(f"[{ci_lower:.4f}, {ci_upper:.4f}]")
print(f"真实值是否在置信区间内: {ci_lower <= true_lambda <= ci_upper}")
# 4. 精确置信区间 (基于卡方分布)
# 对于泊松分布的精确置信区间
sum_x = np.sum(sample_data)
# 下界: χ²_{2Σx, α/2} / (2n)
# 上界: χ²_{2(Σx+1), 1-α/2} / (2n)
chi2_lower = 0.0 if sum_x == 0 else stats.chi2.ppf(0.025, 2 * sum_x)
chi2_upper = stats.chi2.ppf(0.975, 2 * (sum_x + 1))
exact_ci_lower = chi2_lower / (2 * n_samples)
exact_ci_upper = chi2_upper / (2 * n_samples)
print(f"\n精确95%置信区间:")
print(f"[{exact_ci_lower:.4f}, {exact_ci_upper:.4f}]")
# 5. 不同样本量下的估计精度
sample_sizes = [50, 100, 500, 1000, 5000]
print(f"\n不同样本量下的估计精度:")
print("样本量\t估计值\t\t标准误差\t95%置信区间宽度")
print("-" * 60)
for n in sample_sizes:
sample = np.random.poisson(true_lambda, n)
estimate = np.mean(sample)
se_n = np.sqrt(estimate / n)
ci_width = 2 * z_975 * se_n
print(f"{n}\t{estimate:.4f}\t\t{se_n:.4f}\t\t{ci_width:.4f}")
poisson_parameter_estimation()
🎯 实际应用场景
1. 客服呼叫中心分析
def call_center_analysis():
"""呼叫中心的泊松分布应用"""
np.random.seed(42)
# 模拟一天24小时的呼叫数据
# 不同时段的呼叫率不同
hourly_rates = [
1, 1, 0.5, 0.5, 0.5, 1, # 0-5点 (夜间)
3, 5, 8, 10, 12, 12, # 6-11点 (上午)
15, 14, 12, 13, 15, 16, # 12-17点 (下午)
12, 10, 8, 6, 4, 2 # 18-23点 (晚上)
]
print("呼叫中心泊松分布分析:")
print("=" * 50)
# 生成一周的数据 (7天)
daily_calls = []
hourly_calls_all = []
for day in range(7):
daily_total = 0
for hour, rate in enumerate(hourly_rates):
n_calls = np.random.poisson(rate)
hourly_calls_all.append(n_calls)
daily_total += n_calls
daily_calls.append(daily_total)
# 分析结果
print(f"一周总呼叫量: {sum(daily_calls)}")
print(f"日均呼叫量: {np.mean(daily_calls):.1f}")
print(f"日呼叫量标准差: {np.std(daily_calls, ddof=1):.1f}")
# 高峰时段分析 (12-17点)
peak_hours = hourly_rates[12:18]
peak_mean = np.mean(peak_hours)
print(f"\n高峰时段分析 (12-17点):")
print(f"平均每小时呼叫量: {peak_mean:.1f}")
# 服务水平分析
# 假设每个客服每小时能处理8个电话
agents_per_hour = 8
print(f"\n人员配置分析:")
for hour in range(24):
expected_calls = hourly_rates[hour]
min_agents = int(np.ceil(expected_calls / agents_per_hour))
# 假设每人每小时固定处理8通，以小时呼叫量95%分位数估算座席数；不保证等待时间或接通服务水平
# P(X ≤ k) ≥ 0.95
calls_95 = stats.poisson.ppf(0.95, expected_calls)
agents_95 = int(np.ceil(calls_95 / agents_per_hour))
if hour in [9, 10, 14, 15]: # 只显示几个关键时段
print(f"{hour:2d}:00 - 期望{expected_calls:4.1f}通, "
f"按期望工作量估算{min_agents}人, 按小时呼叫量95%分位数估算{agents_95}人（非服务水平保证）")
# 等待时间分析
# 假设服务时间服从指数分布，平均3分钟
service_rate = 20 # 每小时服务20个客户 (3分钟/个)
print(f"\n等待时间分析 (高峰期):")
peak_arrival_rate = 15 # 高峰期每小时15个电话
# 单座席M/M/1稳态模型：泊松到达、独立指数服务时间，且到达率小于服务率；此处服务率20通/小时与前述排班假设不同
utilization = peak_arrival_rate / service_rate
print(f"系统利用率: {utilization:.1%}")
if utilization < 1:
avg_wait_time = utilization / (service_rate * (1 - utilization)) * 60 # 转换为分钟
print(f"预期等待时间: {avg_wait_time:.1f}分钟")
else:
print("系统过载，需要增加座席")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 24小时呼叫率
hours = np.arange(24)
ax1.bar(hours, hourly_rates, alpha=0.7, color='skyblue')
ax1.set_xlabel('小时')
ax1.set_ylabel('平均呼叫量')
ax1.set_title('24小时呼叫率分布')
ax1.grid(True, alpha=0.3)
# 某个高峰时段的实际分布
peak_hour_calls = [np.random.poisson(15) for _ in range(100)] # 模拟100天的14点数据
unique_calls, call_counts = np.unique(peak_hour_calls, return_counts=True)
call_probs = call_counts / len(peak_hour_calls)
# 理论分布
x_theory = np.arange(0, max(peak_hour_calls) + 3)
theory_probs = stats.poisson.pmf(x_theory, 15)
ax2.bar(unique_calls, call_probs, alpha=0.7, width=0.4,
color='lightgreen', label='实际观测')
ax2.plot(x_theory, theory_probs, 'ro-', linewidth=2,
label='理论Poisson(15)')
ax2.set_xlabel('每小时呼叫数')
ax2.set_ylabel('概率')
ax2.set_title('高峰时段呼叫分布 (14:00)')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 日呼叫量分布
days = np.arange(1, 8)
ax3.bar(days, daily_calls, alpha=0.7, color='orange')
ax3.axhline(y=np.mean(daily_calls), color='red', linestyle='--',
label=f'均值 = {np.mean(daily_calls):.1f}')
ax3.set_xlabel('星期')
ax3.set_ylabel('日呼叫量')
ax3.set_title('一周日呼叫量')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 座席需求分析
expected_agents = [int(np.ceil(rate / 8)) for rate in hourly_rates]
agents_95 = [int(np.ceil(stats.poisson.ppf(0.95, rate) / 8)) for rate in hourly_rates]
ax4.plot(hours, expected_agents, 'b-o', linewidth=2, label='最低需求')
ax4.plot(hours, agents_95, 'r-s', linewidth=2, label='小时呼叫量95%分位数对应座席估算')
ax4.fill_between(hours, expected_agents, agents_95, alpha=0.3, color='yellow')
ax4.set_xlabel('小时')
ax4.set_ylabel('所需座席数')
ax4.set_title('座席需求分析')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 业务建议
print(f"\n业务优化建议:")
print("1. 高峰时段 (12-17点) 需要配置更多座席")
print("2. 夜间时段 (0-5点) 可以减少座席或采用自动应答")
print("3. 建议实施弹性排班制，根据历史数据动态调整人员")
print("4. 考虑设置回电服务以降低高峰期等待时间")
call_center_analysis()
2. 网站流量分析
def website_traffic_analysis():
"""网站流量的泊松分布分析"""
np.random.seed(42)
print("网站流量泊松分布分析:")
print("=" * 50)
# 模拟不同页面的访问率 (每分钟访问次数)
pages = {
'首页': 25,
'产品页': 15,
'关于我们': 5,
'联系我们': 3,
'博客': 8
}
# 生成24小时的数据
time_minutes = 24 * 60 # 1440分钟
traffic_data = {}
for page, rate in pages.items():
# 考虑时间因素 - 白天流量高，夜间流量低
hourly_traffic = []
for hour in range(24):
hour_minutes = 60
# 流量波动因子 (模拟白天夜间差异)
if 6 <= hour <= 22: # 白天
factor = 1.0 + 0.5 * np.sin((hour - 6) * np.pi / 16)
else: # 夜间
factor = 0.3
adjusted_rate = rate * factor
# 生成该小时的访问数据
minute_counts = np.random.poisson(adjusted_rate, hour_minutes)
hourly_traffic.extend(minute_counts)
traffic_data[page] = hourly_traffic
# 分析各页面流量
print("页面流量统计:")
print("-" * 40)
total_visits = {}
for page, data in traffic_data.items():
total = sum(data)
mean_per_minute = np.mean(data)
variance = np.var(data, ddof=1)
total_visits[page] = total
print(f"{page}:")
print(f" 总访问量: {total}")
print(f" 每分钟均值: {mean_per_minute:.2f}")
print(f" 方差: {variance:.2f}")
print(f" 均值/方差比: {mean_per_minute/variance:.3f} (仅为离散程度诊断，接近1不能证明泊松拟合；此处全天数据还混合了不同时段的访问率)")
print()
# 服务器负载分析
print("服务器负载分析:")
print("-" * 30)
# 计算每分钟总访问量
total_traffic_per_minute = [sum(traffic_data[page][i] for page in pages)
for i in range(time_minutes)]
max_load = max(total_traffic_per_minute)
avg_load = np.mean(total_traffic_per_minute)
print(f"平均每分钟访问量: {avg_load:.1f}")
print(f"峰值每分钟访问量: {max_load}")
# 容量规划
# 99%的时间内能够处理的访问量
capacity_99 = np.percentile(total_traffic_per_minute, 99)
print(f"本次模拟分钟访问量的经验99%分位数: {capacity_99:.0f}（不保证未来容量覆盖率）")
# 如果服务器容量是平均负载的3倍
server_capacity = avg_load * 3
overload_prob = np.mean(np.array(total_traffic_per_minute) > server_capacity)
print(f"服务器容量设为{server_capacity:.0f}时，过载概率: {overload_prob:.1%}")
# 异常检测
print(f"\n异常流量检测:")
print("-" * 20)
# 基于泊松分布的异常检测
# 如果某分钟的访问量超过95%分位数，视为异常
for page, data in traffic_data.items():
mean_rate = np.mean(data)
# 95%上分位数
threshold_95 = stats.poisson.ppf(0.95, mean_rate)
# 找出异常点
anomalies = [i for i, count in enumerate(data) if count > threshold_95]
if len(anomalies) > 0:
print(f"{page}: 检测到{len(anomalies)}个异常点")
print(f" 阈值: {threshold_95:.0f} (95%分位数)")
print(f" 异常率: {len(anomalies)/len(data):.1%}")
else:
print(f"{page}: 无异常流量")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 各页面总访问量
pages_list = list(total_visits.keys())
visits_list = list(total_visits.values())
ax1.bar(pages_list, visits_list, alpha=0.7, color='skyblue')
ax1.set_ylabel('总访问量')
ax1.set_title('24小时各页面总访问量')
ax1.tick_params(axis='x', rotation=45)
ax1.grid(True, alpha=0.3)
# 首页的24小时流量模式
homepage_hourly = []
for hour in range(24):
hour_start = hour * 60
hour_end = (hour + 1) * 60
hour_total = sum(traffic_data['首页'][hour_start:hour_end])
homepage_hourly.append(hour_total)
hours = np.arange(24)
ax2.plot(hours, homepage_hourly, 'o-', linewidth=2, markersize=6)
ax2.set_xlabel('小时')
ax2.set_ylabel('访问量')
ax2.set_title('首页24小时流量模式')
ax2.grid(True, alpha=0.3)
# 首页某小时的访问分布 (选择第14小时)
hour_14_data = traffic_data['首页'][14*60:15*60]
unique_counts, count_freq = np.unique(hour_14_data, return_counts=True)
count_prob = count_freq / len(hour_14_data)
# 理论泊松分布
mean_rate_14 = np.mean(hour_14_data)
x_theory = np.arange(0, max(hour_14_data) + 3)
theory_prob = stats.poisson.pmf(x_theory, mean_rate_14)
ax3.bar(unique_counts, count_prob, alpha=0.7, width=0.4,
color='lightgreen', label='实际分布')
ax3.plot(x_theory, theory_prob, 'ro-', linewidth=2,
label=f'理论Poisson({mean_rate_14:.1f})')
ax3.set_xlabel('每分钟访问量')
ax3.set_ylabel('概率')
ax3.set_title('首页14:00-15:00访问分布')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 总流量时间序列 (显示前4小时)
time_points = np.arange(0, 240) # 前4小时 = 240分钟
traffic_series = total_traffic_per_minute[:240]
ax4.plot(time_points, traffic_series, alpha=0.7, linewidth=1)
ax4.axhline(y=avg_load, color='red', linestyle='--',
label=f'日均值 = {avg_load:.1f}')
ax4.set_xlabel('时间 (分钟)')
ax4.set_ylabel('总访问量')
ax4.set_title('网站总流量时间序列 (前4小时)')
ax4.legend()
ax4.grid(True, alpha=0.3)
plt.tight_layout()
plt.show()
# 优化建议
print(f"\n网站优化建议:")
print("1. 高峰时段 (上午10点-下午3点) 需要确保服务器性能")
print("2. 可以根据历史流量数据进行自动扩缩容")
print("3. 建议设置实时监控，当流量超过95%分位数时发出警报")
print("4. 对于低流量页面，可以考虑CDN缓存以提高响应速度")
website_traffic_analysis()
3. 质量控制中的缺陷分析
def quality_defect_analysis():
"""制造业缺陷分析的泊松分布应用"""
np.random.seed(42)
print("制造业缺陷泊松分布分析:")
print("=" * 50)
# 模拟不同产品线的缺陷数据
# λ表示每批次平均缺陷数
product_lines = {
'A产品线': 2.1, # 平均每批次2.1个缺陷
'B产品线': 1.5, # 平均每批次1.5个缺陷
'C产品线': 0.8, # 平均每批次0.8个缺陷
'D产品线': 3.2 # 平均每批次3.2个缺陷
}
n_batches = 100 # 每条产品线观察100个批次
defect_data = {}
for product, lambda_rate in product_lines.items():
defects = np.random.poisson(lambda_rate, n_batches)
defect_data[product] = defects
# 基础统计分析
print("各产品线缺陷统计:")
print("-" * 40)
for product, defects in defect_data.items():
mean_defects = np.mean(defects)
var_defects = np.var(defects, ddof=1)
max_defects = np.max(defects)
zero_defect_rate = np.mean(defects == 0)
print(f"{product}:")
print(f" 平均缺陷数: {mean_defects:.2f}")
print(f" 缺陷方差: {var_defects:.2f}")
print(f" 最大缺陷数: {max_defects}")
print(f" 零缺陷率: {zero_defect_rate:.1%}")
print()
# 质量水平评估
print("质量水平评估:")
print("-" * 20)
# 设定质量标准
acceptable_defect_rate = 2.0 # 可接受的平均缺陷率
for product, lambda_rate in product_lines.items():
if lambda_rate <= acceptable_defect_rate:
quality_level = "优秀"
elif lambda_rate <= 3.0:
quality_level = "良好"
elif lambda_rate <= 5.0:
quality_level = "一般"
else:
quality_level = "需要改进"
print(f"{product}: {quality_level} (λ = {lambda_rate})")
# 控制图分析
print(f"\n控制图分析:")
print("-" * 15)
for product, defects in defect_data.items():
lambda_est = np.mean(defects)
# 泊松分布的控制限
# UCL = λ + 3√λ, LCL = max(0, λ - 3√λ)
ucl = lambda_est + 3 * np.sqrt(lambda_est)
lcl = max(0, lambda_est - 3 * np.sqrt(lambda_est))
# 检查失控点
out_of_control = np.sum((defects > ucl) | (defects < lcl))
print(f"{product}:")
print(f" 中心线 CL = {lambda_est:.2f}")
print(f" 上控制限 UCL = {ucl:.2f}")
print(f" 下控制限 LCL = {lcl:.2f}")
print(f" 失控点数量: {out_of_control}")
print()
# 过程能力分析
print("过程能力分析:")
print("-" * 15)
target_defect_rate = 1.0 # 目标缺陷率
for product, lambda_rate in product_lines.items():
# 计算超过目标的概率
prob_exceed_target = 1 - stats.poisson.cdf(target_defect_rate, lambda_rate)
# 计算零缺陷概率
prob_zero_defects = stats.poisson.pmf(0, lambda_rate)
print(f"{product}:")
print(f" 超过目标({target_defect_rate})的概率: {prob_exceed_target:.1%}")
print(f" 零缺陷概率: {prob_zero_defects:.1%}")
print()
# 改进效果预测
print("改进效果预测:")
print("-" * 15)
improvement_factor = 0.7 # 假设改进后缺陷率降低30%
print(f"如果各产品线缺陷率都降低{(1-improvement_factor)*100:.0f}%:")
for product, lambda_rate in product_lines.items():
new_lambda = lambda_rate * improvement_factor
new_zero_prob = stats.poisson.pmf(0, new_lambda)
old_zero_prob = stats.poisson.pmf(0, lambda_rate)
print(f"{product}: λ {lambda_rate:.1f} → {new_lambda:.1f}, "
f"零缺陷率 {old_zero_prob:.1%} → {new_zero_prob:.1%}")
# 可视化
fig, ((ax1, ax2), (ax3, ax4)) = plt.subplots(2, 2, figsize=(15, 12))
# 各产品线缺陷分布对比
products = list(product_lines.keys())
lambdas = list(product_lines.values())
ax1.bar(products, lambdas, alpha=0.7, color='lightcoral')
ax1.axhline(y=acceptable_defect_rate, color='green', linestyle='--',
label=f'可接受水平 = {acceptable_defect_rate}')
ax1.set_ylabel('平均缺陷数 (λ)')
ax1.set_title('各产品线缺陷水平')
ax1.legend()
ax1.tick_params(axis='x', rotation=45)
ax1.grid(True, alpha=0.3)
# A产品线控制图
batch_numbers = np.arange(1, n_batches + 1)
a_defects = defect_data['A产品线']
a_lambda = np.mean(a_defects)
a_ucl = a_lambda + 3 * np.sqrt(a_lambda)
a_lcl = max(0, a_lambda - 3 * np.sqrt(a_lambda))
ax2.plot(batch_numbers, a_defects, 'bo-', markersize=4, alpha=0.7)
ax2.axhline(y=a_lambda, color='green', linewidth=2, label=f'CL = {a_lambda:.2f}')
ax2.axhline(y=a_ucl, color='red', linestyle='--', label=f'UCL = {a_ucl:.2f}')
ax2.axhline(y=a_lcl, color='red', linestyle='--', label=f'LCL = {a_lcl:.2f}')
ax2.set_xlabel('批次号')
ax2.set_ylabel('缺陷数')
ax2.set_title('A产品线控制图')
ax2.legend()
ax2.grid(True, alpha=0.3)
# 缺陷分布直方图 (A产品线)
unique_defects, defect_counts = np.unique(a_defects, return_counts=True)
defect_probs = defect_counts / len(a_defects)
# 理论泊松分布
x_theory = np.arange(0, max(a_defects) + 2)
theory_probs = stats.poisson.pmf(x_theory, a_lambda)
ax3.bar(unique_defects, defect_probs, alpha=0.7, width=0.4,
color='skyblue', label='实际分布')
ax3.plot(x_theory, theory_probs, 'ro-', linewidth=2,
label=f'理论Poisson({a_lambda:.1f})')
ax3.set_xlabel('缺陷数')
ax3.set_ylabel('概率')
ax3.set_title('A产品线缺陷分布')
ax3.legend()
ax3.grid(True, alpha=0.3)
# 零缺陷率对比
zero_defect_rates = []
for product, lambda_rate in product_lines.items():
zero_rate = stats.poisson.pmf(0, lambda_rate)
zero_defect_rates.append(zero_rate)
ax4.bar(products, zero_defect_rates, alpha=0.7, color='lightgreen')
ax4.set_ylabel('零缺陷率')
ax4.set_title('各产品线零缺陷率对比')
ax4.tick_params(axis='x', rotation=45)
ax4.grid(True, alpha=0.3)
# 在柱子上显示数值
for i, rate in enumerate(zero_defect_rates):
ax4.text(i, rate + 0.01, f'{rate:.1%}', ha='center', va='bottom')
plt.tight_layout()
plt.show()
# 管理建议
print(f"\n质量管理建议:")
print("1. D产品线缺陷率过高，需要立即进行根因分析")
print("2. 建立基于泊松分布的实时质量监控系统")
print("3. 设定预警机制，当连续3个批次超过UCL时发出警报")
print("4. 基于稳定基线回顾控制限；更新控制限本身不能确保过程受控，也不应掩盖异常")
print("5. 示例中C产品线的设定缺陷均值较低；真实场景需核对检验规模并调查原因，不能仅据此认定存在可推广的最佳实践")
quality_defect_analysis()
🧮 泊松分布的近似
1. 正态近似
def poisson_normal_approximation():
"""演示泊松分布的正态近似"""
lambda_values = [1, 5, 10, 25, 50]
print("泊松分布的正态近似:")
print("=" * 40)
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
for i, lam in enumerate(lambda_values):
if i >= len(axes):
break
# 泊松分布
x_max = int(lam + 4 * np.sqrt(lam))
x = np.arange(0, x_max + 1)
poisson_pmf = stats.poisson.pmf(x, lam)
# 正态近似 (连续性修正)
x_normal = np.linspace(0, x_max, 1000)
normal_pdf = stats.norm.pdf(x_normal, lam, np.sqrt(lam))
# 绘图
axes[i].bar(x, poisson_pmf, alpha=0.7, width=0.8,
color='skyblue', label=f'Poisson({lam})')
axes[i].plot(x_normal, normal_pdf, 'r-', linewidth=2,
label=f'N({lam}, {lam})')
axes[i].set_title(f'λ = {lam}')
axes[i].set_xlabel('k')
axes[i].set_ylabel('概率/密度')
axes[i].legend()
axes[i].grid(True, alpha=0.3)
# 计算近似误差
# 选择几个点计算误差
test_points = [int(lam), int(lam + np.sqrt(lam)), int(lam - np.sqrt(lam))]
test_points = [max(0, p) for p in test_points if p <= x_max]
print(f"\nλ = {lam} 时的近似误差:")
for point in test_points:
poisson_prob = stats.poisson.pmf(point, lam)
# 连续性修正: P(X = k) ≈ P(k-0.5 < Y < k+0.5)
normal_prob = (stats.norm.cdf(point + 0.5, lam, np.sqrt(lam)) -
stats.norm.cdf(point - 0.5, lam, np.sqrt(lam)))
error = abs(poisson_prob - normal_prob)
print(f" P(X = {point}): 泊松={poisson_prob:.4f}, "
f"正态={normal_prob:.4f}, 误差={error:.4f}")
# 删除多余的子图
if len(lambda_values) < len(axes):
for j in range(len(lambda_values), len(axes)):
fig.delaxes(axes[j])
plt.tight_layout()
plt.show()
# 近似准则
print(f"\n正态近似的使用准则:")
print("1. λ增大时，标准化泊松分布渐近于标准正态分布")
print("2. λ=5、10、30可作为比较近似误差的示例，不是通用精度阈值")
print("3. 是否采用正态近似需结合概率区域及绝对或相对误差要求，尤其应检查尾部概率")
print("4. 使用连续性修正可以提高小λ时的近似精度")
poisson_normal_approximation()
2. 与二项分布的关系
def binomial_to_poisson_limit():
"""演示二项分布收敛到泊松分布"""
# 固定 λ = np = 5
lambda_fixed = 5
# 不同的 n 值，相应调整 p = λ/n
n_values = [10, 50, 100, 500, 1000]
print("二项分布收敛到泊松分布:")
print("=" * 40)
print(f"固定 λ = np = {lambda_fixed}")
fig, axes = plt.subplots(2, 3, figsize=(15, 10))
axes = axes.flatten()
for i, n in enumerate(n_values):
if i >= len(axes):
break
p = lambda_fixed / n
# 计算合理的x范围
x_max = min(n, int(lambda_fixed + 4 * np.sqrt(lambda_fixed)))
x = np.arange(0, x_max + 1)
# 二项分布
binomial_pmf = stats.binom.pmf(x, n, p)
# 泊松分布
poisson_pmf = stats.poisson.pmf(x, lambda_fixed)
# 绘图
axes[i].plot(x, binomial_pmf, 'bo-', markersize=4, alpha=0.7,
label=f'Binomial({n}, {p:.4f})')
axes[i].plot(x, poisson_pmf, 'r^-', markersize=4,
label=f'Poisson({lambda_fixed})')
axes[i].set_title(f'n = {n}, p = {p:.4f}')
axes[i].set_xlabel('k')
axes[i].set_ylabel('P(X = k)')
axes[i].legend()
axes[i].grid(True, alpha=0.3)
# 计算KL散度作为距离度量
# 只在两个分布都有意义的范围内计算
valid_range = (binomial_pmf > 1e-10) & (poisson_pmf > 1e-10)
if np.sum(valid_range) > 0:
kl_div = np.sum(binomial_pmf[valid_range] *
np.log(binomial_pmf[valid_range] / poisson_pmf[valid_range]))
else:
kl_div = np.inf
print(f"n = {n:4d}, p = {p:.6f}: 截断并过滤后的KL项部分和 = {kl_div:.6f}（非完整KL散度）")
# 删除多余的子图
if len(n_values) < len(axes):
for j in range(len(n_values), len(axes)):
fig.delaxes(axes[j])
plt.tight_layout()
plt.show()
# 收敛条件
print(f"\n泊松近似的使用条件:")
print("1. n 很大 (通常 n ≥ 20)")
print("2. p 很小 (通常 p ≤ 0.1)")
print("3. np 为中等大小 (通常 1 ≤ np ≤ 10)")
print("4. 满足以上条件时，Binomial(n,p) ≈ Poisson(np)")
binomial_to_poisson_limit()
🔗 相关知识
#统计学 #泊松分布 #概率分布 #稀有事件 #离散分布 #基础