Python概率统计可视化——概率分布、假设检验与分子运动模型
Python概率统计可视化——概率分布、假设检验与分子运动模型
前言
概率统计作为描述不确定性和随机现象的数学工具，广泛应用于物理学、生物学、经济学等领域。然而，抽象的概率分布和统计推断过程往往难以直观理解。可视化技术通过将概率密度、假设检验逻辑和分子运动规律转化为图形，为数据分析和理论验证提供了直观的洞察途径。本节结合Python的可视化库，实现从基础概率分布到复杂统计模型的可视化，揭示数据背后的统计规律与物理本质。
一、基本概念与理论基础
1. 常见概率分布
概率分布描述随机变量的取值规律，分为离散型和连续型：
-
正态分布（Normal Distribution）：
又称高斯分布，形如钟形曲线，描述大量独立随机变量的统计规律，参数为均值 μ \mu μ 和标准差 σ \sigma σ，记作 N ( μ , σ 2 ) N(\mu, \sigma^2) N(μ,σ2)。
特点：对称性、集中性，68-95-99.7法则描述数据在1σ、2σ、3σ区间的概率。 -
泊松分布（Poisson Distribution）：
描述单位时间内稀有事件的发生次数，参数为 λ \lambda λ（平均发生率），概率质量函数为：
P ( k ) = λ k e − λ k ! P(k) = \frac{\lambda^k e^{-\lambda}}{k!} P(k)=k!λke−λ
应用：交通事故计数、细胞突变频率等。 -
二项分布（Binomial Distribution）：
描述 n n n 次独立试验中成功次数的分布，参数为 n n n（试验次数）和 p p p（单次成功概率），记作 B ( n , p ) B(n, p) B(n,p)。 -
指数分布（Exponential Distribution）：
描述事件发生的时间间隔，参数为 β \beta β（平均间隔时间），概率密度函数为：
f ( x ) = 1 β e − x / β ( x ≥ 0 ) f(x) = \frac{1}{\beta}e^{-x/\beta} \quad (x \geq 0) f(x)=β1e−x/β(x≥0)
特性：无记忆性，常用于可靠性分析。 -
t分布（t-Distribution）：
用于小样本均值的推断，参数为自由度 d f df df，形状类似正态分布但尾部更厚，随 d f df df 增大趋近于标准正态分布。
2. 假设检验基础
假设检验通过样本数据推断总体参数的假设是否成立：
-
零假设（H0）与备择假设（H1）：
零假设通常为“无效应”或“无差异”，如“新药与安慰剂效果相同”；备择假设为待验证的假设，如“新药效果优于安慰剂”。 -
P值与显著性水平（ α \alpha α）：
- P值：在零假设成立时，观察到极端样本结果的概率；
- α \alpha α：预设的拒绝零假设的阈值（通常为0.05），若 P < α P < \alpha P<α 则拒绝H0，可能犯I型错误（弃真错误）。
-
两类错误：
- I型错误：错误拒绝H0，概率为 α \alpha α；
- II型错误：错误接受H0，概率为 β \beta β，检验功效为 1 − β 1-\beta 1−β（正确拒绝H0的概率）。
3. 分子运动统计模型
麦克斯韦-玻尔兹曼分布描述理想气体分子的速率分布，揭示微观运动与宏观热力学的联系：
- 物理意义：
分子速率 v v v 的概率密度函数为：
f ( v ) = 4 π ( m 2 π k T ) 3 / 2 v 2 e − m v 2 2 k T f(v) = 4\pi \left(\frac{m}{2\pi kT}\right)^{3/2} v^2 e^{-\frac{mv^2}{2kT}} f(v)=4π(2πkTm)3/2v2e−2kTmv2
其中 m m m 为分子质量， T T T 为温度， k k k 为玻尔兹曼常数。 - 特征速率：
- 最概然速率（ v p v_p vp）：出现概率最高的速率， v p = 2 k T m v_p = \sqrt{\frac{2kT}{m}} vp=m2kT；
- 平均速率（ v ˉ \bar{v} vˉ）： v ˉ = 8 k T π m \bar{v} = \sqrt{\frac{8kT}{\pi m}} vˉ=πm8kT；
- 方均根速率（ v r m s v_{rms} vrms）： v r m s = 3 k T m v_{rms} = \sqrt{\frac{3kT}{m}} vrms=m3kT。
二、概率分布可视化工具
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson, binom, expon, t
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider, Button
import matplotlib.gridspec as gridspecdef plot_probability_distributions(dist_type='normal', params=None, x_range=None, title=None, output_path="probability_distribution.png"):"""可视化常见概率分布参数:dist_type: 分布类型 ('normal', 'poisson', 'binomial', 'exponential', 't')params: 分布参数 (字典形式)x_range: x轴范围title: 图像标题output_path: 输出文件路径"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 设置默认参数if params is None:params = {}# 根据分布类型设置默认参数和范围if dist_type == 'normal':mu = params.get('mu', 0)sigma = params.get('sigma', 1)if x_range is None:x_range = (mu - 4*sigma, mu + 4*sigma)x = np.linspace(x_range[0], x_range[1], 500)y = norm.pdf(x, mu, sigma)dist_label = f"正态分布: μ={mu}, σ={sigma}"fill_label = "68.27% (1σ)" # 1σ区间填充elif dist_type == 'poisson':mu = params.get('mu', 5)if x_range is None:x_range = (0, max(10, 2*mu))x = np.arange(np.floor(x_range[0]), np.ceil(x_range[1])+1)y = poisson.pmf(x, mu)dist_label = f"泊松分布: λ={mu}"fill_label = Noneelif dist_type == 'binomial':n = params.get('n', 10)p = params.get('p', 0.5)if x_range is None:x_range = (0, n)x = np.arange(0, n+1)y = binom.pmf(x, n, p)dist_label = f"二项分布: n={n}, p={p:.2f}"fill_label = Noneelif dist_type == 'exponential':beta = params.get('beta', 1)if x_range is None:x_range = (0, 4*beta)x = np.linspace(x_range[0], x_range[1], 500)y = expon.pdf(x, scale=beta)dist_label = f"指数分布: β={beta}"fill_label = Noneelif dist_type == 't':df = params.get('df', 5)if x_range is None:x_range = (-4, 4)x = np.linspace(x_range[0], x_range[1], 500)y = t.pdf(x, df)# 对比正态分布y_norm = norm.pdf(x)plt.plot(x, y_norm, 'r--', alpha=0.7, label="标准正态分布")dist_label = f"t分布: df={df}"fill_label = None# 绘制分布曲线plt.plot(x, y, 'b-', linewidth=2.5, label=dist_label)# 添加特殊区域填充if fill_label:if dist_type == 'normal':# 填充1σ区间x_fill = np.linspace(mu-sigma, mu+sigma, 100)y_fill = norm.pdf(x_fill, mu, sigma)plt.fill_between(x_fill, y_fill, color='skyblue', alpha=0.5, label=fill_label)# 标记关键点plt.axvline(mu, color='k', linestyle='--', alpha=0.7)plt.text(mu, norm.pdf(mu, mu, sigma)*1.05, 'μ', ha='center', fontsize=12)# 添加2σ和3σ标记for i in [2, 3]:x_i = mu + i*sigmaplt.axvline(x_i, color='gray', linestyle=':', alpha=0.5)plt.text(x_i, norm.pdf(x_i, mu, sigma)*0.9, f'{i}σ', ha='center', fontsize=10)# 设置标题if title is None:title = f"{dist_label}概率分布"plt.title(title, fontsize=16, pad=15)plt.xlabel("x" if dist_type != 'binomial' else "成功次数", fontsize=12)plt.ylabel("概率密度" if dist_type != 'binomial' else "概率", fontsize=12)# 添加图例plt.legend(loc='best', fontsize=11)# 添加网格plt.grid(True, linestyle='--', alpha=0.4)# 保存图像plt.savefig(output_path, bbox_inches='tight')plt.close()print(f"概率分布图已保存至: {output_path}")
三、假设检验P值动态演示
def interactive_hypothesis_test(output_path="hypothesis_test.gif"):"""创建交互式假设检验P值演示动画"""# 设置参数mu0 = 100 # 零假设均值sigma = 15 # 标准差alpha = 0.05 # 显著性水平n = 30 # 样本量# 创建图形fig, ax = plt.subplots(figsize=(12, 7), dpi=130)plt.subplots_adjust(bottom=0.25)# 初始样本均值sample_mean = 105# 计算抽样分布 (正态分布)se = sigma / np.sqrt(n) # 标准误x = np.linspace(mu0 - 4*se, mu0 + 4*se, 500)y = norm.pdf(x, mu0, se)# 绘制抽样分布line, = plt.plot(x, y, 'b-', linewidth=2.5, label="零假设抽样分布")plt.axvline(mu0, color='k', linestyle='--', alpha=0.7, label=f"零假设均值: μ={mu0}")# 绘制样本均值线sample_line = plt.axvline(sample_mean, color='r', linestyle='-', linewidth=2, label=f"样本均值: {sample_mean:.2f}")# 计算并绘制P值区域p_value = 2 * (1 - norm.cdf(abs(sample_mean - mu0)/se)) # 双侧检验if sample_mean > mu0:x_p = np.linspace(sample_mean, mu0 + 4*se, 100)else:x_p = np.linspace(mu0 - 4*se, sample_mean, 100)y_p = norm.pdf(x_p, mu0, se)p_fill = plt.fill_between(x_p, y_p, color='red', alpha=0.3, label=f"P值: {p_value:.4f}")# 绘制临界值z_crit = norm.ppf(1 - alpha/2)crit_lower = mu0 - z_crit * secrit_upper = mu0 + z_crit * seplt.axvline(crit_lower, color='g', linestyle=':', linewidth=2, label=f"临界值: ±{z_crit:.2f}σ")plt.axvline(crit_upper, color='g', linestyle=':', linewidth=2)# 设置标题和标签plt.title("假设检验P值演示 (双侧检验)", fontsize=16, pad=15)plt.xlabel("样本均值", fontsize=12)plt.ylabel("概率密度", fontsize=12)plt.legend(loc='upper left', fontsize=10)plt.grid(True, linestyle='--', alpha=0.4)# 添加样本均值滑块ax_slider = plt.axes([0.25, 0.1, 0.55, 0.03])slider = Slider(ax_slider, '样本均值', mu0-3*se, mu0+3*se, valinit=sample_mean)# 添加重置按钮ax_button = plt.axes([0.8, 0.05, 0.1, 0.04])button = Button(ax_button, '重置')# 更新函数def update(val):new_mean = slider.val# 更新样本均值线sample_line.set_xdata([new_mean, new_mean])# 更新P值区域p_value = 2 * (1 - norm.cdf(abs(new_mean - mu0)/se))# 清除旧填充global p_fillp_fill.remove()# 绘制新P值区域if new_mean > mu0:x_p = np.linspace(new_mean, mu0 + 4*se, 100)else:x_p = np.linspace(mu0 - 4*se, new_mean, 100)y_p = norm.pdf(x_p, mu0, se)p_fill = plt.fill_between(x_p, y_p, color='red', alpha=0.3, label=f"P值: {p_value:.4f}")# 更新图例handles, labels = ax.get_legend_handles_labels()# 保留除P值外的所有标签keep_labels = [label for label in labels if not label.startswith("P值")]keep_handles = [handles[i] for i, label in enumerate(labels) if not label.startswith("P值")]# 添加新P值标签keep_handles.append(p_fill)keep_labels.append(f"P值: {p_value:.4f}")# 更新图例ax.legend(keep_handles, keep_labels, loc='upper left', fontsize=10)# 判断是否拒绝零假设decision = "拒绝H₀" if p_value < alpha else "不拒绝H₀"plt.title(f"假设检验P值演示 (双侧检验) - {decision}", fontsize=16, pad=15)fig.canvas.draw_idle()# 重置函数def reset(event):slider.set_val(mu0)# 注册事件slider.on_changed(update)button.on_clicked(reset)# 初始更新update(sample_mean)# 保存为GIFdef save_animation():# 创建动画帧means = np.linspace(mu0-3*se, mu0+3*se, 50)means = np.concatenate([means, means[::-1]])def anim_update(frame):slider.set_val(means[frame])return sample_line, p_fillanim = FuncAnimation(fig, anim_update, frames=len(means), interval=100)anim.save(output_path, writer='pillow', fps=10)print(f"假设检验动画已保存至: {output_path}")# 保存动画save_animation()plt.close()def plot_type_errors(output_path="type_errors.png"):"""可视化I型和II型错误"""# 参数设置mu0 = 100 # 零假设均值mu1 = 110 # 备择假设均值sigma = 15 # 总体标准差n = 30 # 样本量alpha = 0.05 # 显著性水平# 计算标准误se = sigma / np.sqrt(n)# 创建图形fig, ax = plt.subplots(figsize=(12, 8), dpi=130)# 零假设分布x0 = np.linspace(mu0 - 4*se, mu0 + 4*se, 500)y0 = norm.pdf(x0, mu0, se)plt.plot(x0, y0, 'b-', linewidth=2.5, label=f"零假设分布 (μ={mu0})")# 备择假设分布x1 = np.linspace(mu1 - 4*se, mu1 + 4*se, 500)y1 = norm.pdf(x1, mu1, se)plt.plot(x1, y1, 'r-', linewidth=2.5, label=f"备择假设分布 (μ={mu1})")# 计算临界值 (双侧检验)z_crit = norm.ppf(1 - alpha/2)crit_lower = mu0 - z_crit * secrit_upper = mu0 + z_crit * se# 绘制临界线plt.axvline(crit_lower, color='g', linestyle=':', linewidth=2, label=f"临界值: {crit_lower:.1f}, {crit_upper:.1f}")plt.axvline(crit_upper, color='g', linestyle=':', linewidth=2)# 填充I型错误区域 (α)x_alpha1 = np.linspace(mu0 - 4*se, crit_lower, 100)y_alpha1 = norm.pdf(x_alpha1, mu0, se)plt.fill_between(x_alpha1, y_alpha1, color='blue', alpha=0.3, label=f"I型错误 (α={alpha})")x_alpha2 = np.linspace(crit_upper, mu0 + 4*se, 100)y_alpha2 = norm.pdf(x_alpha2, mu0, se)plt.fill_between(x_alpha2, y_alpha2, color='blue', alpha=0.3)# 填充II型错误区域 (β)x_beta = np.linspace(crit_lower, crit_upper, 100)y_beta = norm.pdf(x_beta, mu1, se)plt.fill_between(x_beta, y_beta, color='red', alpha=0.3, label="II型错误 (β)")# 计算β值beta = norm.cdf(crit_upper, mu1, se) - norm.cdf(crit_lower, mu1, se)power = 1 - beta# 添加统计信息stats_text = (f"统计量:\n"f"α = {alpha} (显著性水平)\n"f"β = {beta:.3f} (II型错误概率)\n"f"检验功效 = {power:.3f}")plt.text(0.05, 0.95, stats_text, transform=ax.transAxes, fontsize=12,bbox=dict(facecolor='white', alpha=0.8))# 设置标题和标签plt.title("假设检验中的I型和II型错误", fontsize=16, pad=15)plt.xlabel("样本均值", fontsize=12)plt.ylabel("概率密度", fontsize=12)plt.legend(loc='upper left', fontsize=10)plt.grid(True, linestyle='--', alpha=0.4)# 保存图像plt.savefig(output_path, bbox_inches='tight')plt.close()print(f"错误类型图已保存至: {output_path}")
四、跨学科案例：分子运动速率分布
def maxwell_boltzmann_distribution(T=300, m=4.65e-26, output_path="maxwell_boltzmann.png"):"""可视化麦克斯韦-玻尔兹曼速率分布 (气体分子速率分布)参数:T: 温度 (K)m: 分子质量 (kg)output_path: 输出文件路径"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 物理常数k = 1.38e-23 # 玻尔兹曼常数 (J/K)# 计算最概然速率v_p = np.sqrt(2 * k * T / m)# 设置速率范围v_max = 4 * v_pv = np.linspace(0, v_max, 500)# 计算分布函数f = 4 * np.pi * (m/(2 * np.pi * k * T))**(3/2) * v**2 * np.exp(-m * v**2 / (2 * k * T))# 绘制分布曲线plt.plot(v, f, 'b-', linewidth=2.5, label=f"T={T}K")# 标记最概然速率plt.axvline(v_p, color='r', linestyle='--', alpha=0.7)plt.text(v_p, f.max()*0.9, f'最概然速率 $v_p$={v_p:.1f} m/s', ha='center', fontsize=11)# 计算并标记平均速率和方均根速率v_avg = np.sqrt(8 * k * T / (np.pi * m))v_rms = np.sqrt(3 * k * T / m)plt.axvline(v_avg, color='g', linestyle='--', alpha=0.7)plt.text(v_avg, f.max()*0.8, f'平均速率 $\\bar{{v}}$={v_avg:.1f} m/s', ha='center', fontsize=11)plt.axvline(v_rms, color='m', linestyle='--', alpha=0.7)plt.text(v_rms, f.max()*0.7, f'方均根速率 $v_{{rms}}$={v_rms:.1f} m/s', ha='center', fontsize=11)# 设置标题和标签plt.title(f"麦克斯韦-玻尔兹曼速率分布 (分子质量={m:.2e} kg)", fontsize=16, pad=15)plt.xlabel("速率 (m/s)", fontsize=12)plt.ylabel("概率密度", fontsize=12)# 添加网格plt.grid(True, linestyle='--', alpha=0.4)# 添加图例plt.legend(loc='best', fontsize=11)# 添加物理公式formula_text = (r"$f(v) = 4\pi \left(\frac{m}{2\pi k T}\right)^{3/2} v^2 \exp\left(-\frac{m v^2}{2kT}\right)$")plt.text(0.05, 0.95, formula_text, transform=ax.transAxes, fontsize=14,bbox=dict(facecolor='white', alpha=0.8))# 保存图像plt.savefig(output_path, bbox_inches='tight')plt.close()print(f"分子速率分布图已保存至: {output_path}")def plot_temperature_effect(m=4.65e-26, output_path="temperature_effect.png"):"""可视化温度对分子速率分布的影响"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 物理常数k = 1.38e-23 # 玻尔兹曼常数 (J/K)# 不同温度temperatures = [100, 300, 1000] # Kcolors = ['b', 'g', 'r']# 计算最概然速率确定范围v_p_max = np.sqrt(2 * k * max(temperatures) / m)v = np.linspace(0, 1.5 * v_p_max, 500)for T, color in zip(temperatures, colors):# 计算分布函数f = 4 * np.pi * (m/(2 * np.pi * k * T))**(3/2) * v**2 * np.exp(-m * v**2 / (2 * k * T))# 绘制分布曲线plt.plot(v, f, color=color, linewidth=2.5, label=f"T={T}K")# 标记最概然速率v_p = np.sqrt(2 * k * T / m)plt.axvline(v_p, color=color, linestyle='--', alpha=0.5)# 设置标题和标签plt.title("温度对分子速率分布的影响", fontsize=16, pad=15)plt.xlabel("速率 (m/s)", fontsize=12)plt.ylabel("概率密度", fontsize=12)# 添加网格plt.grid(True, linestyle='--', alpha=0.4)# 添加图例plt.legend(loc='best', fontsize=11)# 添加物理解释plt.text(0.6, 0.8, "温度升高:\n- 分布变宽\n- 峰值右移\n- 峰值降低", transform=ax.transAxes, fontsize=12,bbox=dict(facecolor='white', alpha=0.8))# 保存图像plt.savefig(output_path, bbox_inches='tight')plt.close()print(f"温度效应图已保存至: {output_path}")def plot_molecular_speeds(output_path="molecular_speeds.png"):"""可视化不同气体的分子速率分布"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 物理常数k = 1.38e-23 # 玻尔兹曼常数 (J/K)T = 298 # 室温 (298K)# 不同气体 (分子质量, kg)gases = {"氢气 (H₂)": 3.34e-27,"氦气 (He)": 6.64e-27,"氮气 (N₂)": 4.65e-26,"氧气 (O₂)": 5.31e-26,"二氧化碳 (CO₂)": 7.31e-26}colors = plt.cm.viridis(np.linspace(0, 1, len(gases)))# 计算最重气体的最概然速率确定范围heaviest_mass = max(gases.values())v_p_max = np.sqrt(2 * k * T / heaviest_mass)v = np.linspace(0, 1.5 * v_p_max, 500)for (name, mass), color in zip(gases.items(), colors):# 计算分布函数f = 4 * np.pi * (mass/(2 * np.pi * k * T))**(3/2) * v**2 * np.exp(-mass * v**2 / (2 * k * T))# 绘制分布曲线plt.plot(v, f, color=color, linewidth=2.5, label=name)# 设置标题和标签plt.title(f"不同气体的分子速率分布 (T={T}K)", fontsize=16, pad=15)plt.xlabel("速率 (m/s)", fontsize=12)plt.ylabel("概率密度", fontsize=12)# 添加网格plt.grid(True, linestyle='--', alpha=0.4)# 添加图例plt.legend(loc='best', fontsize=11)# 添加物理解释plt.text(0.6, 0.8, "分子质量增加:\n- 分布变窄\n- 峰值左移\n- 峰值升高", transform=ax.transAxes, fontsize=12,bbox=dict(facecolor='white', alpha=0.8))# 保存图像plt.savefig(output_path, bbox_inches='tight')plt.close()print(f"气体分子速率比较图已保存至: {output_path}")
五、使用示例
if __name__ == "__main__":# 1. 概率分布可视化# 正态分布plot_probability_distributions('normal', {'mu': 100, 'sigma': 15}, title="智商(IQ)分布: μ=100, σ=15",output_path="normal_distribution.png")# 泊松分布plot_probability_distributions('poisson', {'mu': 3.5}, title="泊松分布: λ=3.5 (单位时间事件发生次数)",output_path="poisson_distribution.png")# 2. 假设检验演示interactive_hypothesis_test(output_path="hypothesis_test.gif")plot_type_errors(output_path="type_errors.png")# 3. 分子运动速率分布# 氮气分子在室温下的速率分布m_n2 = 4.65e-26 # 氮气分子质量 (kg)maxwell_boltzmann_distribution(T=298, m=m_n2, output_path="n2_velocity_distribution.png")# 温度效应plot_temperature_effect(m=m_n2, output_path="temperature_effect.png")# 不同气体比较plot_molecular_speeds(output_path="gas_velocity_comparison.png")
六、生成图像说明
-
概率分布图 (
normal_distribution.png,poisson_distribution.png):- 展示不同概率分布的形状特征
- 正态分布包含1σ、2σ、3σ区间标记
- 离散分布以柱状图形式显示
- 包含分布参数和关键统计量
-
假设检验动画 (
hypothesis_test.gif):- 动态展示样本均值变化对P值的影响
- 蓝色曲线表示零假设下的抽样分布
- 红色区域表示P值区域
- 绿色虚线表示临界值
- 实时显示检验决策结果
-
错误类型图 (
type_errors.png):- 蓝色曲线表示零假设分布
- 红色曲线表示备择假设分布
- 蓝色区域表示I型错误(α)
- 红色区域表示II型错误(β)
- 包含统计量计算和检验功效
-
分子速率分布图 (
n2_velocity_distribution.png):- 展示麦克斯韦-玻尔兹曼分布曲线
- 标记最概然速率、平均速率和方均根速率
- 包含分布函数的数学公式
- 显示温度和分子质量参数
-
温度效应图 (
temperature_effect.png):- 比较不同温度下的速率分布
- 温度升高导致分布变宽、峰值右移
- 直观展示温度对分子动能的影响
-
气体分子比较图 (
molecular_speeds.png):- 比较不同气体的分子速率分布
- 分子质量增加导致分布变窄、峰值左移
- 使用不同颜色区分气体类型
- 包含分子运动论的物理解释
七、教学要点
-
概率分布基础:
- 正态分布的性质与68-95-99.7法则
- 泊松分布与稀有事件建模
- 离散分布与连续分布的区别
- 分布参数的实际意义
-
假设检验原理:
- 零假设与备择假设的概念
- P值的定义与解释
- I型错误与II型错误的权衡
- 检验功效与样本量的关系
-
统计物理学应用:
- 麦克斯韦-玻尔兹曼分布的推导
- 温度与分子动能的统计关系
- 最概然速率、平均速率和方均根速率的计算
- 分子质量对速率分布的影响
-
科学可视化技巧:
- 概率分布的直观表示
- 动态交互式演示的实现
- 复杂概念的图解方法
- 跨学科知识的整合展示
- 物理公式的LaTeX渲染
本节提供了强大的概率统计可视化工具，从基础分布到高级假设检验概念。分子运动速率分布案例展示了统计学在物理化学中的实际应用，揭示了宏观现象背后的微观统计规律。
总结
通过概率分布可视化，我们将抽象的统计规律转化为可观测的图形；假设检验的动态演示揭示了统计推断的逻辑本质；分子运动模型则展现了统计学在物理学中的跨学科应用。这些工具不仅帮助理解概率统计的核心概念，更为数据分析、科学建模提供了直观的方法论。可视化作为连接数学理论与实际应用的桥梁，将复杂的随机现象转化为人类认知可理解的形式，是现代科学研究的重要辅助手段。
完整代码
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm, poisson, binom, expon, t
from matplotlib.animation import FuncAnimation
from matplotlib.widgets import Slider, Button
import matplotlib.gridspec as gridspec
import os
import glob# ---------------------- 全局配置 ----------------------
plt.rcParams["font.sans-serif"] = ["SimHei", "WenQuanYi Micro Hei", "Heiti TC"]
plt.rcParams["axes.unicode_minus"] = False # 正确显示负号
plt.rcParams["mathtext.fontset"] = "cm" # 设置数学字体# ====================== 概率分布可视化 ======================
def plot_probability_distributions(dist_type="normal",params=None,x_range=None,title=None,output_path="probability_distribution.png",
):"""可视化常见概率分布"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()params = params or {}if dist_type == "normal":mu = params.get("mu", 0)sigma = params.get("sigma", 1)if x_range is None:x_range = (mu - 4 * sigma, mu + 4 * sigma)x = np.linspace(x_range[0], x_range[1], 500)y = norm.pdf(x, mu, sigma)plt.plot(x, y, "b-", linewidth=2.5, label=f"正态分布: μ={mu}, σ={sigma}")# 填充1σ区间x_fill = np.linspace(mu - sigma, mu + sigma, 100)plt.fill_between(x_fill,norm.pdf(x_fill, mu, sigma),color="skyblue",alpha=0.5,label="68.27%区间",)# 标记关键点plt.axvline(mu, color="k", linestyle="--", alpha=0.7)plt.text(mu, norm.pdf(mu, mu, sigma) * 1.05, "μ", ha="center", fontsize=12)# 添加2σ和3σ标记for i in [2, 3]:x_i = mu + i * sigmaplt.axvline(x_i, color="gray", linestyle=":", alpha=0.5)plt.text(x_i, norm.pdf(x_i, mu, sigma) * 0.9, f"{i}σ", ha="center", fontsize=10)elif dist_type == "poisson":mu = params.get("mu", 5)if x_range is None:x_range = (0, max(10, 2 * mu))x = np.arange(int(x_range[0]), int(x_range[1]) + 1)y = poisson.pmf(x, mu)plt.stem(x,y,linefmt="b-",markerfmt="bo",basefmt="none",label=f"泊松分布: λ={mu}",)elif dist_type == "binomial":n = params.get("n", 10)p = params.get("p", 0.5)if x_range is None:x_range = (0, n)x = np.arange(0, n + 1)y = binom.pmf(x, n, p)plt.stem(x,y,linefmt="b-",markerfmt="bo",basefmt="none",label=f"二项分布: n={n}, p={p:.2f}",)elif dist_type == "exponential":beta = params.get("beta", 1)if x_range is None:x_range = (0, 4 * beta)x = np.linspace(x_range[0], x_range[1], 500)y = expon.pdf(x, scale=beta)plt.plot(x, y, "b-", linewidth=2.5, label=f"指数分布: β={beta}")elif dist_type == "t":df = params.get("df", 5)if x_range is None:x_range = (-4, 4)x = np.linspace(x_range[0], x_range[1], 500)y = t.pdf(x, df)plt.plot(x, y, "b-", linewidth=2.5, label=f"t分布: df={df}")# 对比正态分布y_norm = norm.pdf(x)plt.plot(x, y_norm, "r--", alpha=0.7, label="标准正态分布")plt.title(title or f"{dist_type.capitalize()}分布可视化", fontsize=16)plt.xlabel("变量值" if dist_type != "binomial" else "成功次数")plt.ylabel("概率密度" if dist_type != "binomial" else "概率")plt.legend(loc="best", fontsize=11)plt.grid(True, linestyle="--", alpha=0.4)plt.savefig(output_path, bbox_inches="tight")plt.close()print(f"已生成{dist_type}分布图: {output_path}")# ====================== 假设检验可视化 ======================
def interactive_hypothesis_test(output_path="hypothesis_test.gif"):mu0 = 100 # 零假设均值sigma = 15 # 标准差alpha = 0.05 # 显著性水平n = 30 # 样本量fig, ax = plt.subplots(figsize=(12, 7), dpi=130)plt.subplots_adjust(bottom=0.25)se = sigma / np.sqrt(n)x = np.linspace(mu0 - 4 * se, mu0 + 4 * se, 500)plt.plot(x, norm.pdf(x, mu0, se), "b-", linewidth=2.5, label=f"零假设分布 (μ={mu0})")plt.axvline(mu0, color="k", linestyle="--", alpha=0.7, label="零假设均值 H0") # 修改为H0sample_mean = 105sample_line = plt.axvline(sample_mean,color="r",linestyle="-",linewidth=2,label=f"样本均值: {sample_mean}",)z_crit = norm.ppf(1 - alpha / 2)crit_lower = mu0 - z_crit * secrit_upper = mu0 + z_crit * seplt.axvline(crit_lower,color="g",linestyle=":",linewidth=2,label=f"临界值: {crit_lower:.1f}, {crit_upper:.1f}",)plt.title("假设检验P值演示 (双侧检验)", fontsize=16)plt.xlabel("样本均值")plt.ylabel("概率密度")plt.legend(loc="upper left", fontsize=10)plt.grid(True, linestyle="--", alpha=0.4)means = np.linspace(mu0 - 3 * se, mu0 + 3 * se, 50)means = np.concatenate([means, means[::-1]])def update(frame):current_mean = means[frame]sample_line.set_xdata([current_mean, current_mean])p_value = 2 * (1 - norm.cdf(abs(current_mean - mu0) / se))decision = "拒绝H0" if p_value < alpha else "不拒绝H0" # 修改为H0plt.title(f"假设检验P值演示 (双侧检验) - {decision}, P值={p_value:.4f}", fontsize=16)return (sample_line,)anim = FuncAnimation(fig, update, frames=len(means), interval=100)anim.save(output_path, writer="pillow", fps=10)plt.close()print(f"假设检验动画已保存: {output_path}")def plot_type_errors(output_path="type_errors.png"):mu0 = 100 # 零假设均值mu1 = 110 # 备择假设均值sigma = 15 # 总体标准差n = 30 # 样本量alpha = 0.05 # 显著性水平se = sigma / np.sqrt(n)fig, ax = plt.subplots(figsize=(12, 8), dpi=130)# 零假设分布（修改为H0）x0 = np.linspace(mu0 - 4 * se, mu0 + 4 * se, 500)y0 = norm.pdf(x0, mu0, se)plt.plot(x0, y0, "b-", linewidth=2.5, label=f"H0分布 (μ={mu0})")# 备择假设分布（修改为H1）x1 = np.linspace(mu1 - 4 * se, mu1 + 4 * se, 500)y1 = norm.pdf(x1, mu1, se)plt.plot(x1, y1, "r-", linewidth=2.5, label=f"H1分布 (μ={mu1})")z_crit = norm.ppf(1 - alpha / 2)crit_lower = mu0 - z_crit * secrit_upper = mu0 + z_crit * seplt.axvline(crit_lower,color="g",linestyle=":",linewidth=2,label=f"临界值: {crit_lower:.1f}, {crit_upper:.1f}",)# 填充I型错误区域x_alpha1 = np.linspace(mu0 - 4 * se, crit_lower, 100)y_alpha1 = norm.pdf(x_alpha1, mu0, se)plt.fill_between(x_alpha1, y_alpha1, color="blue", alpha=0.3, label=f"I型错误 (α={alpha})")x_alpha2 = np.linspace(crit_upper, mu0 + 4 * se, 100)y_alpha2 = norm.pdf(x_alpha2, mu0, se)plt.fill_between(x_alpha2, y_alpha2, color="blue", alpha=0.3)# 填充II型错误区域x_beta = np.linspace(crit_lower, crit_upper, 100)y_beta = norm.pdf(x_beta, mu1, se)plt.fill_between(x_beta, y_beta, color="red", alpha=0.3, label="II型错误 (β)")# 计算并添加统计信息（移至右上角）beta = norm.cdf(crit_upper, mu1, se) - norm.cdf(crit_lower, mu1, se)power = 1 - betastats_text = (f"统计量:\n"f"α = {alpha} (显著性水平)\n"f"β = {beta:.3f} (II型错误概率)\n"f"检验功效 = {power:.3f}")plt.text(0.95,0.95,stats_text,transform=ax.transAxes,fontsize=12,bbox=dict(facecolor="white", alpha=0.8),ha="right",) # 右上角plt.title("假设检验中的I型和II型错误", fontsize=16)plt.xlabel("样本均值")plt.ylabel("概率密度")plt.legend(loc="upper left", fontsize=10)plt.grid(True, linestyle="--", alpha=0.4)plt.savefig(output_path, bbox_inches="tight")plt.close()print(f"错误类型图已保存: {output_path}")# ====================== 分子运动模型 ======================
def maxwell_boltzmann_distribution(T=300, m=4.65e-26, output_path="maxwell_boltzmann.png"
):"""可视化麦克斯韦-玻尔兹曼速率分布"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 物理常数k = 1.38e-23 # 玻尔兹曼常数# 计算最概然速率vp = np.sqrt(2 * k * T / m)# 设置速率范围v_max = 4 * vpv = np.linspace(0, v_max, 500)# 计算分布函数f = (4* np.pi* (m / (2 * np.pi * k * T)) ** (3 / 2)* v**2* np.exp(-m * v**2 / (2 * k * T)))# 绘制分布曲线plt.plot(v, f, "b-", linewidth=2.5, label=f"T={T}K")# 标记最概然速率plt.axvline(vp, color="r", linestyle="--", alpha=0.7)plt.text(vp, f.max() * 0.9, f"最概然速率: {vp:.1f} m/s", ha="center", fontsize=11)# 计算并标记平均速率和方均根速率v_avg = np.sqrt(8 * k * T / (np.pi * m))v_rms = np.sqrt(3 * k * T / m)plt.axvline(v_avg, color="g", linestyle="--", alpha=0.7)plt.text(v_avg, f.max() * 0.8, f"平均速率: {v_avg:.1f} m/s", ha="center", fontsize=11)plt.axvline(v_rms, color="m", linestyle="--", alpha=0.7)plt.text(v_rms, f.max() * 0.7, f"方均根速率: {v_rms:.1f} m/s", ha="center", fontsize=11)# 设置标题和标签plt.title(f"麦克斯韦-玻尔兹曼速率分布 (分子质量={m:.2e} kg)", fontsize=16)plt.xlabel("速率 (m/s)")plt.ylabel("概率密度")# 添加网格和图例plt.grid(True, linestyle="--", alpha=0.4)plt.legend(loc="best", fontsize=11)# 添加物理公式formula_text = r"$f(v) = 4\pi \left(\frac{m}{2\pi k T}\right)^{3/2} v^2 \exp\left(-\frac{m v^2}{2kT}\right)$"plt.text(0.05,0.95,formula_text,transform=ax.transAxes,fontsize=14,bbox=dict(facecolor="white", alpha=0.8),)# 保存图像plt.savefig(output_path, bbox_inches="tight")plt.close()print(f"分子速率分布图已保存: {output_path}")def plot_temperature_effect(m=4.65e-26, output_path="temperature_effect.png"):"""可视化温度对分子速率分布的影响"""plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()# 物理常数k = 1.38e-23 # 玻尔兹曼常数# 不同温度temperatures = [100, 300, 1000] # Kcolors = ["b", "g", "r"]# 计算最概然速率确定范围v_p_max = np.sqrt(2 * k * max(temperatures) / m)v = np.linspace(0, 1.5 * v_p_max, 500)for T, color in zip(temperatures, colors):# 计算分布函数f = (4* np.pi* (m / (2 * np.pi * k * T)) ** (3 / 2)* v**2* np.exp(-m * v**2 / (2 * k * T)))# 绘制分布曲线plt.plot(v, f, color=color, linewidth=2.5, label=f"T={T}K")# 标记最概然速率v_p = np.sqrt(2 * k * T / m)plt.axvline(v_p, color=color, linestyle="--", alpha=0.5)# 设置标题和标签plt.title("温度对分子速率分布的影响", fontsize=16)plt.xlabel("速率 (m/s)")plt.ylabel("概率密度")# 添加网格和图例plt.grid(True, linestyle="--", alpha=0.4)plt.legend(loc="best", fontsize=11)# 添加物理解释plt.text(0.6,0.8,"温度升高:\n- 分布变宽\n- 峰值右移\n- 峰值降低",transform=ax.transAxes,fontsize=12,bbox=dict(facecolor="white", alpha=0.8),)# 保存图像plt.savefig(output_path, bbox_inches="tight")plt.close()print(f"温度效应图已保存: {output_path}")def plot_molecular_speeds(output_path="molecular_speeds.png"):plt.figure(figsize=(12, 8), dpi=130)ax = plt.gca()k = 1.38e-23 # 玻尔兹曼常数T = 298 # 室温 (298K)# 修改为普通文本，如CO2替代CO₂gases = {"氢气 (H2)": 3.34e-27,"氦气 (He)": 6.64e-27,"氮气 (N2)": 4.65e-26,"氧气 (O2)": 5.31e-26,"二氧化碳 (CO2)": 7.31e-26,}colors = plt.cm.viridis(np.linspace(0, 1, len(gases)))heaviest_mass = max(gases.values())v_p_max = np.sqrt(2 * k * T / heaviest_mass)v = np.linspace(0, 1.5 * v_p_max, 500)for (name, mass), color in zip(gases.items(), colors):f = (4* np.pi* (mass / (2 * np.pi * k * T)) ** (3 / 2)* v**2* np.exp(-mass * v**2 / (2 * k * T)))plt.plot(v, f, color=color, linewidth=2.5, label=name)plt.title(f"不同气体的分子速率分布 (T={T}K)", fontsize=16)plt.xlabel("速率 (m/s)")plt.ylabel("概率密度")plt.grid(True, linestyle="--", alpha=0.4)plt.legend(loc="best", fontsize=11)plt.text(0.6,0.8,"分子质量增加:\n- 分布变窄\n- 峰值左移\n- 峰值升高",transform=ax.transAxes,fontsize=12,bbox=dict(facecolor="white", alpha=0.8),)plt.savefig(output_path, bbox_inches="tight")plt.close()print(f"气体分子速率比较图已保存: {output_path}")# ====================== 使用示例 ======================
if __name__ == "__main__":# 清理旧文件for f in ["*.png", "*.gif"]:for file in glob.glob(f):os.remove(file)# 1. 概率分布可视化plot_probability_distributions("normal",{"mu": 100, "sigma": 15},title="智商(IQ)分布: μ=100, σ=15",output_path="normal_distribution.png",)plot_probability_distributions("poisson",{"mu": 3.5},title="泊松分布: λ=3.5 (单位时间事件发生次数)",output_path="poisson_distribution.png",)# 2. 假设检验演示interactive_hypothesis_test()plot_type_errors()# 3. 分子运动速率分布# 氮气分子在室温下的速率分布m_n2 = 4.65e-26 # 氮气分子质量 (kg)maxwell_boltzmann_distribution(T=298, m=m_n2, output_path="n2_velocity_distribution.png")# 温度效应plot_temperature_effect(m=m_n2)# 不同气体比较plot_molecular_speeds()print("\n所有图形已成功生成！文件列表：")for f in os.listdir():if any(f.endswith(ext) for ext in [".png", ".gif"]):print(f"- {f}")