用Python实战推导泊松分布与正态分布的特征函数：从数学公式到代码实现

在概率论与统计学的学习中，特征函数是一个强大但常被忽视的工具。它不仅能完整描述随机变量的概率分布，还能简化许多复杂计算。本文将带你用Python代码一步步推导泊松分布和正态分布的特征函数，并通过可视化验证其性质。

1. 特征函数：概念与Python基础实现

特征函数是概率分布的一种变换表示，定义为随机变量X的复指数期望：φ(t) = E[e^(itX)]。这个看似简单的定义蕴含着丰富的信息：

• 唯一性 ：特征函数与概率分布一一对应
• 矩生成 ：可通过导数计算各阶矩
• 卷积简化 ：独立随机变量和的特征函数等于特征函数的乘积

让我们先用Python定义特征函数的基础计算框架：

import numpy as np
import matplotlib.pyplot as plt
from sympy import symbols, exp, I, integrate, oo, sqrt, pi

def characteristic_function(dist, t):
    """
    计算给定分布的特征函数
    :param dist: 分布对象，需实现expectation方法
    :param t: 特征函数参数
    :return: 特征函数值
    """
    return dist.expectation(lambda x: np.exp(1j * t * x))

2. 泊松分布的特征函数推导与验证

泊松分布是描述稀有事件发生次数的经典模型。设X~Poisson(λ)，其概率质量函数为：

P(X=k) = (λ^k e^{-λ})/k!, k=0,1,2,...

2.1 理论推导

特征函数的数学推导过程：

φ(t) = E[e^(itX)] = Σ_{k=0}^∞ e^(itk) (λ^k e^{-λ})/k!
= e^{-λ} Σ_{k=0}^∞ (λe^{it})^k /k!
= e^{-λ} e^{λe^{it}}
= exp(λ(e^{it}-1))

2.2 Python实现与验证

from scipy.stats import poisson

class PoissonDistribution:
    def __init__(self, lam):
        self.lam = lam
        
    def expectation(self, func):
        # 数值计算期望，截断到足够大的k值
        max_k = int(self.lam * 10)  # 确保覆盖主要概率质量
        ks = np.arange(0, max_k+1)
        probs = poisson.pmf(ks, self.lam)
        return np.sum(func(ks) * probs)

# 验证特征函数
lam = 3.0
poisson_dist = PoissonDistribution(lam)
t_values = np.linspace(-5, 5, 500)

# 数值计算的特征函数
numeric_phi = [characteristic_function(poisson_dist, t) for t in t_values]

# 理论特征函数
theoretical_phi = np.exp(lam * (np.exp(1j * t_values) - 1))

# 可视化比较
plt.figure(figsize=(12, 6))
plt.plot(t_values, np.real(numeric_phi), 'b-', label='数值实部')
plt.plot(t_values, np.real(theoretical_phi), 'r--', label='理论实部')
plt.plot(t_values, np.imag(numeric_phi), 'g-', label='数值虚部')
plt.plot(t_values, np.imag(theoretical_phi), 'm--', label='理论虚部')
plt.title('泊松分布特征函数验证(λ=3)')
plt.xlabel('t')
plt.ylabel('φ(t)')
plt.legend()
plt.grid(True)
plt.show()

运行结果将显示数值计算与理论公式完美吻合，验证了我们的推导。

3. 正态分布的特征函数推导与验证

标准正态分布N(0,1)的特征函数推导是概率论中的经典案例。让我们用Python再现这一过程。

3.1 理论推导

φ(t) = ∫_{-∞}^∞ e^{itx} (1/√(2π)) e^{-x^2/2} dx
= e^{-t^2/2} ∫_{-∞}^∞ (1/√(2π)) e^{-(x-it)^2/2} dx
= e^{-t^2/2}

关键步骤利用了复变函数中的围道积分技巧。

3.2 Python实现与验证

from scipy.stats import norm

class NormalDistribution:
    def __init__(self, mu=0, sigma=1):
        self.mu = mu
        self.sigma = sigma
        
    def expectation(self, func):
        # 使用数值积分计算期望
        x = np.linspace(-10, 10, 10000)  # 足够大的范围
        pdf = norm.pdf(x, self.mu, self.sigma)
        return np.trapz(func(x) * pdf, x)

# 验证标准正态分布的特征函数
normal_dist = NormalDistribution()
t_values = np.linspace(-5, 5, 500)

# 数值计算的特征函数
numeric_phi = [characteristic_function(normal_dist, t) for t in t_values]

# 理论特征函数
theoretical_phi = np.exp(-t_values**2 / 2)

# 可视化比较
plt.figure(figsize=(12, 6))
plt.plot(t_values, np.real(numeric_phi), 'b-', label='数值实部')
plt.plot(t_values, np.real(theoretical_phi), 'r--', label='理论实部')
plt.plot(t_values, np.imag(numeric_phi), 'g-', label='数值虚部')
plt.title('标准正态分布特征函数验证')
plt.xlabel('t')
plt.ylabel('φ(t)')
plt.legend()
plt.grid(True)
plt.show()

4. 特征函数的应用：矩计算与分布再生性

特征函数不仅是一个理论概念，更有强大的实用价值。让我们通过Python实现两个典型应用。

4.1 矩的计算

特征函数的k阶导数在0点的值与随机变量的k阶矩相关：

E[X^k] = i^{-k} φ^{(k)}(0)

from scipy.misc import derivative

def compute_moment(dist, k, delta=1e-5):
    """
    使用特征函数计算随机变量的k阶矩
    """
    def phi(t):
        return characteristic_function(dist, t)
    
    # 计算k阶导数
    dk_phi = derivative(phi, 0, dx=delta, n=k, order=2*k+1)
    return (1j)**(-k) * dk_phi

# 计算泊松分布(λ=3)的前三阶矩
lam = 3.0
poisson_dist = PoissonDistribution(lam)

moments = [compute_moment(poisson_dist, k) for k in range(1, 4)]
print(f"泊松分布(λ={lam})的矩：")
print(f"一阶矩(期望): {moments[0]:.4f} (理论值: {lam})")
print(f"二阶矩: {moments[1]:.4f} (理论值: {lam + lam**2})")
print(f"三阶矩: {np.real(moments[2]):.4f} (理论值: {lam + 3*lam**2 + lam**3})")

4.2 分布再生性的验证

泊松分布具有再生性：X~Poisson(λ1), Y~Poisson(λ2)独立，则X+Y~Poisson(λ1+λ2)。我们可以用特征函数验证：

# 生成两个独立的泊松随机变量
lam1, lam2 = 2.0, 3.0
poisson_dist1 = PoissonDistribution(lam1)
poisson_dist2 = PoissonDistribution(lam2)

# 计算各自特征函数的乘积
t_values = np.linspace(-5, 5, 500)
phi1 = np.array([characteristic_function(poisson_dist1, t) for t in t_values])
phi2 = np.array([characteristic_function(poisson_dist2, t) for t in t_values])
phi_sum = phi1 * phi2

# 理论上的泊松(λ1+λ2)特征函数
theoretical_phi_sum = np.exp((lam1 + lam2) * (np.exp(1j * t_values) - 1))

# 可视化比较
plt.figure(figsize=(12, 6))
plt.plot(t_values, np.real(phi_sum), 'b-', label='乘积实部')
plt.plot(t_values, np.real(theoretical_phi_sum), 'r--', label='理论实部')
plt.plot(t_values, np.imag(phi_sum), 'g-', label='乘积虚部')
plt.plot(t_values, np.imag(theoretical_phi_sum), 'm--', label='理论虚部')
plt.title('泊松分布再生性验证(λ1=2, λ2=3)')
plt.xlabel('t')
plt.ylabel('φ(t)')
plt.legend()
plt.grid(True)
plt.show()

5. 高级应用：从特征函数反推分布

虽然实践中不常用，但理论上可以从特征函数反推概率分布。让我们实现一个简单的逆变换示例。

def inverse_transform(phi, t_values, x):
    """
    使用数值积分从特征函数反推概率密度
    """
    integrand = np.exp(-1j * t_values * x) * phi
    f_x = np.trapz(np.real(integrand), t_values) / (2 * np.pi)
    return f_x

# 对标准正态分布进行测试
t_samples = np.linspace(-50, 50, 10000)  # 需要足够大的范围
phi_normal = np.exp(-t_samples**2 / 2)   # 标准正态特征函数

# 反推几个点的概率密度
x_points = [-2, -1, 0, 1, 2]
for x in x_points:
    estimated = inverse_transform(phi_normal, t_samples, x)
    true_value = norm.pdf(x)
    print(f"x={x}: 估计值={estimated:.6f}, 真实值={true_value:.6f}")

这个简单的实现展示了特征函数与概率分布之间的对偶关系，虽然数值精度受限于积分范围和采样密度，但验证了理论的有效性。