概率与统计

概率与统计是处理机器人学中不确定性的数学基础。本章将讲解机器人学所需的概率统计知识。

学习目标

完成本章后，你将能够：

• 理解概率论的基本概念
• 掌握常用概率分布
• 理解贝叶斯推断
• 应用概率统计解决机器人学问题

1. 概率基础

1.1 概率公理

概率的三条公理：
1. 非负性: P(A) ≥ 0
2. 规范性: P(Ω) = 1 (Ω 是样本空间)
3. 可加性: 对于互斥事件 A₁, A₂, ...
   P(A₁ ∪ A₂ ∪ ...) = P(A₁) + P(A₂) + ...

常用公式：
- 条件概率: P(A|B) = P(A∩B) / P(B)
- 乘法公式: P(A∩B) = P(A|B) * P(B)
- 全概率公式: P(A) = Σ P(A|Bᵢ) * P(Bᵢ)
- 贝叶斯公式: P(B|A) = P(A|B) * P(B) / P(A)

1.2 概率计算

import numpy as np

# 概率质量函数 (PMF) - 离散随机变量
def pmf_binomial(n, k, p):
    """二项分布 PMF"""
    from math import comb
    return comb(n, k) * (p ** k) * ((1 - p) ** (n - k))

# 概率密度函数 (PDF) - 连续随机变量
def pdf_gaussian(x, mu, sigma):
    """高斯分布 PDF"""
    return (1 / (sigma * np.sqrt(2 * np.pi))) * \
           np.exp(-0.5 * ((x - mu) / sigma) ** 2)

# 累积分布函数 (CDF)
def cdf_gaussian(x, mu, sigma):
    """高斯分布 CDF"""
    from scipy import stats
    return stats.norm.cdf(x, mu, sigma)

# 示例
# 二项分布：抛硬币 10 次，正面朝上 7 次的概率
n, k, p = 10, 7, 0.5
prob = pmf_binomial(n, k, p)
print(f"P(X=7) = {prob:.4f}")

# 高斯分布：x=1.5 处的概率密度
mu, sigma = 0, 1
x = 1.5
density = pdf_gaussian(x, mu, sigma)
print(f"f(1.5) = {density:.4f}")

# 高斯分布：x<1.5 的概率
prob_less = cdf_gaussian(x, mu, sigma)
print(f"P(X < 1.5) = {prob_less:.4f}")

2. 常用概率分布

2.1 离散分布

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

# 伯努利分布
# X ∈ {0, 1}, P(X=1) = p
p = 0.3
bernoulli = stats.bernoulli(p)
print("伯努利分布:")
print(f"  均值: {bernoulli.mean()}")
print(f"  方差: {bernoulli.var()}")

# 二项分布
# X ~ Bin(n, p)
n, p = 10, 0.5
binomial = stats.binom(n, p)
print("\n二项分布:")
print(f"  均值: {binomial.mean()}")
print(f"  方差: {binomial.var()}")

# 泊松分布
# X ~ Pois(λ)
lam = 3
poisson = stats.poisson(lam)
print("\n泊松分布:")
print(f"  均值: {poisson.mean()}")
print(f"  方差: {poisson.var()}")

# 绘制分布
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 伯努利分布
x = [0, 1]
axes[0].bar(x, bernoulli.pmf(x))
axes[0].set_title('Bernoulli Distribution')
axes[0].set_xlabel('x')
axes[0].set_ylabel('P(x)')

# 二项分布
x = np.arange(0, n+1)
axes[1].bar(x, binomial.pmf(x))
axes[1].set_title('Binomial Distribution')
axes[1].set_xlabel('x')
axes[1].set_ylabel('P(x)')

# 泊松分布
x = np.arange(0, 10)
axes[2].bar(x, poisson.pmf(x))
axes[2].set_title('Poisson Distribution')
axes[2].set_xlabel('x')
axes[2].set_ylabel('P(x)')

plt.tight_layout()
plt.show()

2.2 连续分布

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

# 高斯分布 (正态分布)
# X ~ N(μ, σ²)
mu, sigma = 0, 1
gaussian = stats.norm(mu, sigma)
print("高斯分布:")
print(f"  均值: {gaussian.mean()}")
print(f"  方差: {gaussian.var()}")
print(f"  标准差: {gaussian.std()}")

# 均匀分布
# X ~ U(a, b)
a, b = 0, 1
uniform = stats.uniform(a, b-a)
print("\n均匀分布:")
print(f"  均值: {uniform.mean()}")
print(f"  方差: {uniform.var()}")

# 指数分布
# X ~ Exp(λ)
lam = 0.5
exponential = stats.expon(lam)
print("\n指数分布:")
print(f"  均值: {exponential.mean()}")
print(f"  方差: {exponential.var()}")

# 绘制分布
fig, axes = plt.subplots(1, 3, figsize=(15, 4))

# 高斯分布
x = np.linspace(-4, 4, 100)
axes[0].plot(x, gaussian.pdf(x))
axes[0].set_title('Gaussian Distribution')
axes[0].set_xlabel('x')
axes[0].set_ylabel('f(x)')

# 均匀分布
x = np.linspace(-0.5, 1.5, 100)
axes[1].plot(x, uniform.pdf(x))
axes[1].set_title('Uniform Distribution')
axes[1].set_xlabel('x')
axes[1].set_ylabel('f(x)')

# 指数分布
x = np.linspace(0, 10, 100)
axes[2].plot(x, exponential.pdf(x))
axes[2].set_title('Exponential Distribution')
axes[2].set_xlabel('x')
axes[2].set_ylabel('f(x)')

plt.tight_layout()
plt.show()

2.3 多元高斯分布

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

# 多元高斯分布
# X ~ N(μ, Σ)
mu = np.array([0, 0])  # 均值向量
Sigma = np.array([[1, 0.5],   # 协方差矩阵
                  [0.5, 2]])

# 创建分布
multivariate_gaussian = stats.multivariate_normal(mu, Sigma)

# 采样
samples = multivariate_gaussian.rvs(1000)

# 绘制
plt.figure(figsize=(8, 6))
plt.scatter(samples[:, 0], samples[:, 1], alpha=0.5)
plt.xlabel('X1')
plt.ylabel('X2')
plt.title('Multivariate Gaussian Distribution')
plt.grid(True)
plt.axis('equal')
plt.show()

# 概率密度
x = np.array([1, 1])
density = multivariate_gaussian.pdf(x)
print(f"在点 {x} 处的概率密度: {density:.4f}")

3. 期望与方差

3.1 期望

import numpy as np

# 离散随机变量的期望
# E[X] = Σ x * P(X=x)
def expected_value_discrete(x, p):
    """计算离散随机变量的期望"""
    return np.sum(x * p)

# 连续随机变量的期望
# E[X] = ∫ x * f(x) dx
def expected_value_continuous(f, a, b):
    """计算连续随机变量的期望"""
    from scipy import integrate
    result, _ = integrate.quad(lambda x: x * f(x), a, b)
    return result

# 示例
# 离散：骰子的期望值
x = np.array([1, 2, 3, 4, 5, 6])
p = np.array([1/6, 1/6, 1/6, 1/6, 1/6, 1/6])
E_X = expected_value_discrete(x, p)
print(f"骰子的期望值: {E_X}")

# 连续：高斯分布的期望
from scipy.stats import norm
mu, sigma = 0, 1
E_X = expected_value_continuous(lambda x: norm.pdf(x, mu, sigma), -10, 10)
print(f"高斯分布的期望: {E_X:.4f}")

3.2 方差

import numpy as np

# 方差
# Var(X) = E[(X - μ)²] = E[X²] - (E[X])²
def variance_discrete(x, p):
    """计算离散随机变量的方差"""
    E_X = np.sum(x * p)
    E_X2 = np.sum(x**2 * p)
    return E_X2 - E_X**2

# 标准差
def std_discrete(x, p):
    """计算离散随机变量的标准差"""
    return np.sqrt(variance_discrete(x, p))

# 示例
x = np.array([1, 2, 3, 4, 5, 6])
p = np.array([1/6, 1/6, 1/6, 1/6, 1/6, 1/6])

Var_X = variance_discrete(x, p)
Std_X = std_discrete(x, p)

print(f"方差: {Var_X:.4f}")
print(f"标准差: {Std_X:.4f}")

3.3 协方差

import numpy as np

# 协方差
# Cov(X, Y) = E[(X - μₓ)(Y - μᵧ)]
def covariance(X, Y):
    """计算协方差"""
    E_X = np.mean(X)
    E_Y = np.mean(Y)
    return np.mean((X - E_X) * (Y - E_Y))

# 相关系数
# ρ = Cov(X, Y) / (σₓ * σᵧ)
def correlation(X, Y):
    """计算相关系数"""
    cov = covariance(X, Y)
    std_X = np.std(X)
    std_Y = np.std(Y)
    return cov / (std_X * std_Y)

# 协方差矩阵
def covariance_matrix(data):
    """计算协方差矩阵"""
    return np.cov(data.T)

# 示例
np.random.seed(42)
X = np.random.randn(100)
Y = 2 * X + np.random.randn(100) * 0.5

cov_XY = covariance(X, Y)
corr_XY = correlation(X, Y)

print(f"协方差: {cov_XY:.4f}")
print(f"相关系数: {corr_XY:.4f}")

# 协方差矩阵
data = np.column_stack([X, Y])
cov_matrix = covariance_matrix(data)
print(f"\n协方差矩阵:")
print(cov_matrix)

4. 贝叶斯推断

4.1 贝叶斯定理

贝叶斯定理：
P(A|B) = P(B|A) * P(A) / P(B)

其中：
- P(A|B): 后验概率 (posterior)
- P(B|A): 似然 (likelihood)
- P(A): 先验概率 (prior)
- P(B): 证据 (evidence)

import numpy as np

def bayes_theorem(prior, likelihood, evidence):
    """贝叶斯定理"""
    posterior = (likelihood * prior) / evidence
    return posterior

# 示例：疾病检测
# 假设：
# - 疾病患病率: 1%
# - 检测准确率: 99% (真阳性)
# - 假阳性率: 5%

# 先验：患病概率
P_disease = 0.01

# 似然：检测为阳性 | 患病
P_positive_given_disease = 0.99

# 证据：检测为阳性的总概率
# P(positive) = P(positive|disease) * P(disease) + P(positive|no_disease) * P(no_disease)
P_positive = P_positive_given_disease * P_disease + 0.05 * (1 - P_disease)

# 后验：患病 | 检测为阳性
P_disease_given_positive = bayes_theorem(P_disease, P_positive_given_disease, P_positive)

print(f"患病概率 (先验): {P_disease:.4f}")
print(f"检测为阳性后患病概率 (后验): {P_disease_given_positive:.4f}")

4.2 贝叶斯更新

import numpy as np

def bayesian_update(prior, likelihood, data):
    """贝叶斯更新"""
    # 计算后验
    posterior = prior * likelihood

    # 归一化
    posterior = posterior / np.sum(posterior)

    return posterior

# 示例：估计硬币正面概率
# 假设：正面概率 θ 可能是 0.3, 0.5, 0.7

# 先验：均匀分布
theta_values = np.array([0.3, 0.5, 0.7])
prior = np.array([1/3, 1/3, 1/3])

# 数据：抛硬币 10 次，7 次正面
n_trials = 10
n_heads = 7

# 似然：P(data|θ) = C(n,k) * θ^k * (1-θ)^(n-k)
from scipy.stats import binom
likelihood = binom.pmf(n_heads, n_trials, theta_values)

# 贝叶斯更新
posterior = bayesian_update(prior, likelihood, None)

print("先验概率:", prior)
print("似然:", likelihood)
print("后验概率:", posterior)

# 多次更新
data = [7, 8, 6, 9, 7]  # 每次 10 次试验中的正面次数

current_prior = prior
for n_heads in data:
    likelihood = binom.pmf(n_heads, n_trials, theta_values)
    current_prior = bayesian_update(current_prior, likelihood, None)
    print(f"\n数据: {n_heads}/10 正面")
    print(f"后验: {current_prior}")

4.3 贝叶斯滤波

import numpy as np

class BayesianFilter:
    """贝叶斯滤波器"""

    def __init__(self, state_space, transition_model, observation_model):
        """
        state_space: 状态空间
        transition_model: 状态转移模型 P(x_t | x_{t-1})
        observation_model: 观测模型 P(z_t | x_t)
        """
        self.state_space = state_space
        self.transition_model = transition_model
        self.observation_model = observation_model

        # 初始概率（均匀分布）
        self.belief = np.ones(len(state_space)) / len(state_space)

    def predict(self):
        """预测步骤"""
        new_belief = np.zeros_like(self.belief)

        for i, x_new in enumerate(self.state_space):
            for j, x_old in enumerate(self.state_space):
                new_belief[i] += self.transition_model(x_new, x_old) * self.belief[j]

        self.belief = new_belief

    def update(self, observation):
        """更新步骤"""
        for i, x in enumerate(self.state_space):
            self.belief[i] *= self.observation_model(observation, x)

        # 归一化
        self.belief /= np.sum(self.belief)

    def get_estimate(self):
        """获取估计值"""
        return np.sum(self.state_space * self.belief)

# 示例：机器人定位
# 状态空间：一维位置 [0, 1, 2, 3, 4]
state_space = np.array([0, 1, 2, 3, 4])

# 转移模型：随机移动
def transition_model(x_new, x_old):
    if abs(x_new - x_old) <= 1:
        return 1/3
    return 0

# 观测模型：带噪声的位置观测
def observation_model(z, x):
    if z == x:
        return 0.7
    elif abs(z - x) == 1:
        return 0.15
    return 0

# 创建滤波器
bf = BayesianFilter(state_space, transition_model, observation_model)

# 模拟观测
observations = [1, 2, 2, 3]

print("贝叶斯滤波:")
print(f"初始置信度: {bf.belief}")

for z in observations:
    bf.predict()
    bf.update(z)
    print(f"观测 {z} 后: {bf.belief}")
    print(f"位置估计: {bf.get_estimate():.2f}")

5. 统计推断

5.1 参数估计

import numpy as np
from scipy import stats

# 最大似然估计 (MLE)
def mle_gaussian(data):
    """高斯分布的最大似然估计"""
    mu_mle = np.mean(data)
    sigma_mle = np.std(data, ddof=0)  # 总体标准差
    return mu_mle, sigma_mle

# 矩估计
def moment_estimation(data):
    """矩估计"""
    mu = np.mean(data)
    sigma2 = np.var(data, ddof=0)
    return mu, sigma2

# 示例
np.random.seed(42)
data = np.random.normal(2, 1.5, 100)  # 真实参数: μ=2, σ=1.5

# MLE
mu_mle, sigma_mle = mle_gaussian(data)
print("最大似然估计:")
print(f"  μ = {mu_mle:.4f}")
print(f"  σ = {sigma_mle:.4f}")

# 矩估计
mu_mom, sigma2_mom = moment_estimation(data)
print("\n矩估计:")
print(f"  μ = {mu_mom:.4f}")
print(f"  σ² = {sigma2_mom:.4f}")

5.2 假设检验

import numpy as np
from scipy import stats

# t 检验
def t_test(sample1, sample2, alpha=0.05):
    """双样本 t 检验"""
    t_stat, p_value = stats.ttest_ind(sample1, sample2)

    print(f"t 统计量: {t_stat:.4f}")
    print(f"p 值: {p_value:.4f}")

    if p_value < alpha:
        print(f"拒绝原假设 (α={alpha})")
    else:
        print(f"接受原假设 (α={alpha})")

    return t_stat, p_value

# 卡方检验
def chi2_test(observed, expected, alpha=0.05):
    """卡方检验"""
    chi2_stat, p_value = stats.chisquare(observed, expected)

    print(f"卡方统计量: {chi2_stat:.4f}")
    print(f"p 值: {p_value:.4f}")

    if p_value < alpha:
        print(f"拒绝原假设 (α={alpha})")
    else:
        print(f"接受原假设 (α={alpha})")

    return chi2_stat, p_value

# 示例
np.random.seed(42)
sample1 = np.random.normal(0, 1, 50)
sample2 = np.random.normal(0.5, 1, 50)

print("t 检验:")
t_test(sample1, sample2)

5.3 置信区间

import numpy as np
from scipy import stats

def confidence_interval_mean(data, confidence=0.95):
    """均值的置信区间"""
    n = len(data)
    mean = np.mean(data)
    std = np.std(data, ddof=1)

    # t 分布的临界值
    t_critical = stats.t.ppf((1 + confidence) / 2, df=n-1)

    # 置信区间
    margin_of_error = t_critical * std / np.sqrt(n)
    ci_lower = mean - margin_of_error
    ci_upper = mean + margin_of_error

    return ci_lower, ci_upper

# 示例
np.random.seed(42)
data = np.random.normal(2, 1, 100)

ci_lower, ci_upper = confidence_interval_mean(data, confidence=0.95)
print(f"样本均值: {np.mean(data):.4f}")
print(f"95% 置信区间: [{ci_lower:.4f}, {ci_upper:.4f}]")

6. 应用示例

6.1 传感器融合

import numpy as np

def sensor_fusion(measurements, uncertainties):
    """传感器融合（加权平均）"""
    # 计算权重（方差的倒数）
    weights = 1 / np.array(uncertainties) ** 2

    # 归一化权重
    weights = weights / np.sum(weights)

    # 加权平均
    fused = np.sum(weights * measurements)

    # 融合后的不确定性
    fused_uncertainty = 1 / np.sqrt(np.sum(1 / np.array(uncertainties) ** 2))

    return fused, fused_uncertainty

# 示例：多个传感器测量同一物理量
measurements = [10.2, 10.5, 10.1, 10.3]  # 测量值
uncertainties = [0.5, 0.3, 0.8, 0.4]      # 不确定性

fused_value, fused_uncertainty = sensor_fusion(measurements, uncertainties)

print("传感器融合:")
for i, (m, u) in enumerate(zip(measurements, uncertainties)):
    print(f"  传感器 {i+1}: {m} ± {u}")
print(f"\n融合结果: {fused_value:.4f} ± {fused_uncertainty:.4f}")

6.2 卡尔曼滤波

import numpy as np

class KalmanFilter:
    """卡尔曼滤波器"""

    def __init__(self, A, B, H, Q, R, x0, P0):
        """
        A: 状态转移矩阵
        B: 控制输入矩阵
        H: 观测矩阵
        Q: 过程噪声协方差
        R: 观测噪声协方差
        x0: 初始状态
        P0: 初始协方差
        """
        self.A = A
        self.B = B
        self.H = H
        self.Q = Q
        self.R = R
        self.x = x0
        self.P = P0

    def predict(self, u=None):
        """预测步骤"""
        # 状态预测
        if u is not None:
            self.x = self.A @ self.x + self.B @ u
        else:
            self.x = self.A @ self.x

        # 协方差预测
        self.P = self.A @ self.P @ self.A.T + self.Q

    def update(self, z):
        """更新步骤"""
        # 计算卡尔曼增益
        S = self.H @ self.P @ self.H.T + self.R
        K = self.P @ self.H.T @ np.linalg.inv(S)

        # 状态更新
        y = z - self.H @ self.x  # 残差
        self.x = self.x + K @ y

        # 协方差更新
        I = np.eye(len(self.x))
        self.P = (I - K @ self.H) @ self.P

    def get_state(self):
        """获取状态估计"""
        return self.x

    def get_covariance(self):
        """获取协方差"""
        return self.P

# 示例：一维位置跟踪
dt = 1.0  # 时间步长

# 状态转移矩阵 [位置, 速度]
A = np.array([[1, dt],
              [0, 1]])

# 控制输入矩阵
B = np.array([[0.5 * dt**2],
              [dt]])

# 观测矩阵（只观测位置）
H = np.array([[1, 0]])

# 过程噪声
Q = np.array([[1, 0],
              [0, 1]])

# 观测噪声
R = np.array([[10]])

# 初始状态
x0 = np.array([0, 0])
P0 = np.array([[100, 0],
               [0, 100]])

# 创建卡尔曼滤波器
kf = KalmanFilter(A, B, H, Q, R, x0, P0)

# 模拟数据
np.random.seed(42)
true_positions = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]
measurements = [p + np.random.normal(0, np.sqrt(10)) for p in true_positions]

# 卡尔曼滤波
print("卡尔曼滤波:")
print(f"{'步骤':<5} {'真实位置':<10} {'测量值':<10} {'估计位置':<10} {'不确定性':<10}")

for i, (true_pos, z) in enumerate(zip(true_positions, measurements)):
    kf.predict()
    kf.update(np.array([z]))

    estimated_pos = kf.get_state()[0]
    uncertainty = np.sqrt(kf.get_covariance()[0, 0])

    print(f"{i:<5} {true_pos:<10.2f} {z:<10.2f} {estimated_pos:<10.2f} {uncertainty:<10.2f}")

6.3 粒子滤波

import numpy as np

class ParticleFilter:
    """粒子滤波器"""

    def __init__(self, n_particles, state_space, motion_model, measurement_model):
        """
        n_particles: 粒子数量
        state_space: 状态空间
        motion_model: 运动模型
        measurement_model: 测量模型
        """
        self.n_particles = n_particles
        self.state_space = state_space
        self.motion_model = motion_model
        self.measurement_model = measurement_model

        # 初始化粒子（均匀分布）
        self.particles = np.random.choice(state_space, n_particles)
        self.weights = np.ones(n_particles) / n_particles

    def predict(self):
        """预测步骤"""
        # 根据运动模型移动粒子
        self.particles = self.motion_model(self.particles)

    def update(self, measurement):
        """更新步骤"""
        # 计算权重
        for i in range(self.n_particles):
            self.weights[i] = self.measurement_model(measurement, self.particles[i])

        # 归一化权重
        self.weights /= np.sum(self.weights)

    def resample(self):
        """重采样步骤"""
        # 系统重采样
        indices = np.random.choice(
            self.n_particles,
            size=self.n_particles,
            p=self.weights
        )
        self.particles = self.particles[indices]
        self.weights = np.ones(self.n_particles) / self.n_particles

    def get_estimate(self):
        """获取估计值"""
        return np.average(self.particles, weights=self.weights)

# 示例
state_space = np.array([0, 1, 2, 3, 4, 5, 6, 7, 8, 9])

def motion_model(particles):
    """运动模型：随机移动"""
    return particles + np.random.choice([-1, 0, 1], size=len(particles))

def measurement_model(z, x):
    """测量模型：高斯噪声"""
    return np.exp(-0.5 * (z - x) ** 2 / 1)

# 创建粒子滤波器
pf = ParticleFilter(1000, state_space, motion_model, measurement_model)

# 模拟
true_position = 3
measurements = [3.2, 2.8, 3.5, 3.1, 2.9]

print("粒子滤波:")
for z in measurements:
    pf.predict()
    pf.update(z)
    pf.resample()

    estimate = pf.get_estimate()
    print(f"测量: {z:.1f}, 估计: {estimate:.2f}")

7. 练习题

基础练习

1. 计算伯努利分布的期望和方差
2. 使用贝叶斯定理计算条件概率
3. 计算样本的均值和方差

进阶练习

1. 实现多元高斯分布的采样
2. 实现贝叶斯更新
3. 计算置信区间

应用练习

1. 实现传感器融合算法
2. 实现卡尔曼滤波器
3. 实现粒子滤波器

下一步

• 优化理论 — 学习优化方法
• 控制理论 — 学习控制方法
• SLAM — 应用概率方法进行同时定位与建图

参考资源

• Probability and Statistics for Engineers
• Bayesian Reasoning and Machine Learning
• Kalman Filter Tutorial

---

← 返回索引 | 下一页：优化理论 →