高斯分布乘积的可视化探索:用Python与Matlab构建概率直觉
在机器人定位和传感器融合领域,理解高斯分布的乘积行为是掌握卡尔曼滤波核心思想的关键。当我们面对来自不同传感器的测量数据时,如何融合这些信息并量化不确定性?传统教材往往从数学推导入手,而本文将带你通过代码绘图直接观察这一过程的几何意义。
想象一下,你正在设计一个自动驾驶汽车的定位系统。GPS提供的位置估计可能具有较大的方差(低精度),而激光雷达的测量则方差较小(高精度)。当这两个分布相乘时,融合后的结果会如何变化?通过可视化,我们不仅能验证数学公式,更能培养对概率分布的"直觉"——这种直觉在实际工程决策中往往比纯理论更有价值。
1. 环境准备与基础分布绘制
1.1 工具链配置
我们推荐使用Python生态进行实验,因其丰富的可视化库和交互性。以下是基础环境配置步骤:
# 创建并激活虚拟环境(推荐) python -m venv gauss_visual source gauss_visual/bin/activate # Linux/Mac gauss_visual\Scripts\activate # Windows # 安装核心依赖 pip install numpy matplotlib ipython scipy对于习惯Matlab的用户,可直接使用内置的统计与机器学习工具箱。两种工具的核心函数对比如下:
| 功能 | Python (SciPy) | Matlab |
|---|---|---|
| 高斯分布PDF | scipy.stats.norm.pdf | normpdf |
| 数组生成 | numpy.linspace | linspace |
| 图形绘制 | matplotlib.pyplot.plot | plot |
| 图形显示 | plt.show() | 自动显示 |
1.2 绘制单个高斯分布
让我们从基础开始,定义一个高斯分布并绘制其曲线。高斯分布由两个参数完全确定:
- μ(均值):分布的中心位置
- σ²(方差):分布的离散程度
import numpy as np import matplotlib.pyplot as plt from scipy.stats import norm def plot_gaussian(mean, variance, label=None): sigma = np.sqrt(variance) x = np.linspace(mean - 4*sigma, mean + 4*sigma, 200) y = norm.pdf(x, mean, sigma) plt.plot(x, y, label=label) plt.figure(figsize=(10, 6)) plot_gaussian(mean=2, variance=1.5, label="N(2, 1.5)") plt.xlabel("Value") plt.ylabel("Probability Density") plt.title("Single Gaussian Distribution") plt.legend() plt.grid(True) plt.show()执行这段代码将显示一个中心在2.0、方差1.5的高斯钟形曲线。调整mean和variance参数,观察曲线如何变化——这是建立分布参数直觉的第一步。
2. 两个高斯分布的乘积可视化
2.1 定义并绘制原始分布
现在我们考虑两个不同的高斯分布相乘的情况。这在传感器融合中非常常见——每个传感器提供的数据都可以看作是一个高斯分布。
# 定义两个高斯分布的参数 mu1, var1 = 2.0, 1.5 # 分布1:均值2.0,方差1.5 mu2, var2 = 3.5, 0.8 # 分布2:均值3.5,方差0.8 # 绘制两个原始分布 plt.figure(figsize=(10, 6)) plot_gaussian(mu1, var1, label=f"N({mu1}, {var1})") plot_gaussian(mu2, var2, label=f"N({mu2}, {var2})") plt.title("Two Original Gaussian Distributions") plt.legend() plt.grid(True) plt.show()你会看到两个重叠的钟形曲线,一个较宽(方差大),一个较窄(方差小)。思考一个问题:它们的乘积曲线会是什么形状?会比原始曲线更宽还是更窄?均值会偏向哪一边?
2.2 计算并绘制乘积分布
高斯分布乘积的数学表达式告诉我们,两个高斯分布N(μ₁,σ₁²)和N(μ₂,σ₂²)的乘积仍然是高斯分布,其参数为:
- 新方差:σ² = (σ₁² * σ₂²) / (σ₁² + σ₂²)
- 新均值:μ = (μ₁σ₂² + μ₂σ₁²) / (σ₁² + σ₂²)
让我们用代码实现这个计算:
def multiply_gaussians(mu1, var1, mu2, var2): """计算两个高斯分布的乘积""" new_var = (var1 * var2) / (var1 + var2) new_mean = (mu1 * var2 + mu2 * var1) / (var1 + var2) return new_mean, new_var # 计算乘积分布 mu_prod, var_prod = multiply_gaussians(mu1, var1, mu2, var2) # 绘制三个分布对比 plt.figure(figsize=(10, 6)) plot_gaussian(mu1, var1, label=f"N({mu1}, {var1})") plot_gaussian(mu2, var2, label=f"N({mu2}, {var2})") plot_gaussian(mu_prod, var_prod, label=f"Product N({mu_prod:.2f}, {var_prod:.2f})") plt.title("Product of Two Gaussian Distributions") plt.legend() plt.grid(True) plt.show()观察结果图,你会发现几个关键现象:
- 乘积分布的方差比两个原始分布都小——这意味着融合后的不确定性降低了
- 乘积分布的均值位于两个原始均值之间,但更靠近方差较小的分布(本例中N(3.5,0.8))
提示:在传感器融合中,这意味着更精确的测量(方差小)在融合结果中会获得更大的权重
3. 参数变化对乘积分布的影响
3.1 方差比例的影响
固定两个分布的均值,只改变它们的方差比例,观察乘积分布如何变化:
# 固定均值,变化方差比例 mu1, mu2 = 2.0, 3.0 var_ratios = [0.1, 0.5, 1, 2, 10] # var2/var1的比例 plt.figure(figsize=(12, 8)) for ratio in var_ratios: var1 = 1.0 var2 = var1 * ratio mu_prod, var_prod = multiply_gaussians(mu1, var1, mu2, var2) # 绘制乘积分布 plot_gaussian(mu_prod, var_prod, label=f"var2/var1={ratio:.1f}, μ={mu_prod:.2f}, σ²={var_prod:.2f}") plt.title("Product Distribution under Different Variance Ratios") plt.legend() plt.grid(True) plt.show()这个实验揭示了几个重要规律:
- 当var2远小于var1(ratio→0)时,乘积分布几乎完全由N(μ₂,σ₂²)决定
- 当var1=var2时,乘积分布的均值正好是两个原始均值的中间点
- 随着ratio增大,乘积均值逐渐向μ₁移动
3.2 均值分离程度的影响
现在固定方差,改变两个均值之间的距离,观察乘积分布的变化:
# 固定方差,变化均值距离 var1, var2 = 1.0, 1.0 mean_distances = [0.5, 1.0, 2.0, 3.0, 4.0] # μ2 - μ1的距离 plt.figure(figsize=(12, 8)) for dist in mean_distances: mu1 = 0.0 mu2 = mu1 + dist mu_prod, var_prod = multiply_gaussians(mu1, var1, mu2, var2) # 绘制乘积分布 plot_gaussian(mu_prod, var_prod, label=f"μ2-μ1={dist}, μ={mu_prod:.2f}, σ²={var_prod:.2f}") plt.title("Product Distribution under Different Mean Separations") plt.legend() plt.grid(True) plt.show()实验结果展示了:
- 当两个均值接近时,乘积分布的方差显著减小
- 随着均值距离增加,乘积分布的峰值降低(因为分布重叠区域减少)
- 当均值距离足够大时,乘积分布几乎为零——这意味着两个测量结果矛盾,不能简单融合
4. 交互式可视化与卡尔曼滤波联系
4.1 创建交互式可视化
为了更直观地理解参数变化的影响,我们可以创建交互式可视化。使用IPython的交互式控件:
from ipywidgets import interact, FloatSlider def interactive_gaussian_product(mu1=0.0, var1=1.0, mu2=2.0, var2=1.0): plt.figure(figsize=(10, 6)) # 绘制原始分布 plot_gaussian(mu1, var1, label=f"N({mu1}, {var1})") plot_gaussian(mu2, var2, label=f"N({mu2}, {var2})") # 计算并绘制乘积分布 mu_prod, var_prod = multiply_gaussians(mu1, var1, mu2, var2) plot_gaussian(mu_prod, var_prod, label=f"Product N({mu_prod:.2f}, {var_prod:.2f})") plt.title("Interactive Gaussian Product Visualization") plt.legend() plt.grid(True) plt.show() # 创建交互式控件 interact(interactive_gaussian_product, mu1=FloatSlider(min=-5, max=5, step=0.1, value=0), var1=FloatSlider(min=0.1, max=3, step=0.1, value=1), mu2=FloatSlider(min=-5, max=5, step=0.1, value=2), var2=FloatSlider(min=0.1, max=3, step=0.1, value=1))拖动滑块实时观察分布变化,你会对高斯乘积的行为建立更牢固的直觉。特别注意:
- 调整方差时,乘积分布如何"偏向"方差较小的分布
- 当两个分布相距很远时,乘积几乎消失(矛盾信息)
4.2 与卡尔曼滤波的联系
卡尔曼滤波的测量更新步骤本质上就是两个高斯分布的乘积:
- 预测状态分布:N(μ_pred, σ_pred²)
- 测量分布:N(μ_meas, σ_meas²)
- 更新后的状态就是这两个分布的乘积
# 卡尔曼滤波更新步骤的简化实现 def kalman_update(prior_mean, prior_var, meas_mean, meas_var): # 计算卡尔曼增益 K = prior_var / (prior_var + meas_var) # 更新均值和方差 new_mean = prior_mean + K * (meas_mean - prior_mean) new_var = (1 - K) * prior_var return new_mean, new_var # 验证与高斯乘积等价 prior_mean, prior_var = 2.0, 1.5 meas_mean, meas_var = 3.5, 0.8 # 两种方法计算结果应该相同 print("Kalman update:", kalman_update(prior_mean, prior_var, meas_mean, meas_var)) print("Gaussian product:", multiply_gaussians(prior_mean, prior_var, meas_mean, meas_var))这个联系解释了为什么卡尔曼滤波能有效融合不确定信息——它在数学上等价于寻找两个高斯分布的重叠部分,这正是我们可视化展示的内容。