用Python的NumPy和SciPy库，5分钟搞定概率论期末试卷里的分布律与期望计算

概率论作为数学的重要分支，在理工科学习中占据着核心地位。然而，传统的手工计算不仅耗时耗力，还容易出错。本文将展示如何利用Python的NumPy和SciPy科学计算库，快速解决概率论试卷中的典型问题，包括离散/连续分布计算、期望方差求解以及联合分布分析。

1. 环境准备与基础工具介绍

在开始之前，我们需要确保安装了必要的Python库。这三个核心库将构成我们的主要工具集：

pip install numpy scipy matplotlib

NumPy提供高效的数组运算能力，SciPy包含了丰富的统计分布函数，而Matplotlib则用于可视化展示。这三个库的组合，能够覆盖概率论计算中的绝大多数需求。

SciPy的stats模块特别值得关注，它内置了超过100种概率分布，每种分布都提供了如下方法：

• pmf/pdf: 概率质量/密度函数
• cdf: 累积分布函数
• rvs: 随机变量生成
• stats: 计算期望、方差等统计量
• fit: 参数估计

2. 离散型随机变量的快速求解

试卷中常见的离散分布问题，如二项分布、泊松分布等，都可以用几行代码解决。让我们看一个典型例子：

问题：设离散型随机变量X的分布律为：

X   -2   -1    0    1    2
P   2a    a  1/8  a/2   5a

求：(1)常数a (2)P{-1<X<2} (3)Y=2X²+1的分布律

import numpy as np
from scipy import stats

# (1) 求解常数a
probs = np.array([2, 1, 1/8, 0.5, 5])  # 系数数组
a = 1 / np.sum(probs)  # 利用概率总和为1的性质
print(f"a = {a:.4f}")

# (2) 计算P{-1<X<2}
x_values = np.array([-2, -1, 0, 1, 2])
p_values = np.array([2*a, a, 1/8, a/2, 5*a])
mask = (x_values > -1) & (x_values < 2)
prob = np.sum(p_values[mask])
print(f"P{{-1<X<2}} = {prob:.4f}")

# (3) 求Y=2X²+1的分布律
y_values = 2 * x_values**2 + 1
unique_y = np.unique(y_values)
y_probs = [np.sum(p_values[y_values == y]) for y in unique_y]

print("\nY的分布律：")
for y, p in zip(unique_y, y_probs):
    print(f"P(Y={int(y)}) = {p:.4f}")

这种方法不仅快速准确，还能自动处理更复杂的变换问题。对于常见的离散分布，SciPy已经提供了现成的实现：

# 二项分布示例
n, p = 10, 0.3
binom_dist = stats.binom(n, p)
print(f"P(X=3) = {binom_dist.pmf(3):.4f}")
print(f"E(X) = {binom_dist.mean():.2f}, Var(X) = {binom_dist.var():.2f}")

# 泊松分布示例
lambda_ = 3
poisson_dist = stats.poisson(lambda_)
print(f"P(X=5) = {poisson_dist.pmf(5):.4f}")

3. 连续型随机变量的高效处理

连续型随机变量的计算通常涉及积分运算，手工计算尤为繁琐。Python可以轻松解决这类问题。

典型问题：设随机变量X的密度函数为：

p(x) = { 2x² + a, 0 < x < 1
       { 0, 其他

求：(1)常数a (2)E(2X+1) (3)分布函数F(x)

from scipy.integrate import quad

# (1) 求常数a
def integrand(x):
    return 2*x**2

total_prob, _ = quad(integrand, 0, 1)
a = 1 - total_prob
print(f"a = {a:.4f}")

# (2) 计算E(2X+1)
def integrand(x):
    return (2*x + 1) * (2*x**2 + a)

expectation, _ = quad(integrand, 0, 1)
print(f"E(2X+1) = {expectation:.4f}")

# (3) 分布函数F(x)
def cdf(x):
    if x <= 0:
        return 0
    elif x >= 1:
        return 1
    else:
        result, _ = quad(lambda t: 2*t**2 + a, 0, x)
        return result

# 验证几个点
for x in [0.2, 0.5, 0.8]:
    print(f"F({x}) = {cdf(x):.4f}")

对于常见的连续分布，SciPy.stats提供了更简洁的使用方式：

# 正态分布示例
mu, sigma = 1, 2
norm_dist = stats.norm(mu, sigma)

# 计算P(1≤X≤5)
prob = norm_dist.cdf(5) - norm_dist.cdf(1)
print(f"P(1≤X≤5) = {prob:.4f}")

# 均匀分布示例
a, b = 0, 2
uniform_dist = stats.uniform(a, b-a)
print(f"E(X) = {uniform_dist.mean():.2f}, D(X) = {uniform_dist.var():.2f}")

# 指数分布示例
lambda_ = 1
exp_dist = stats.expon(scale=1/lambda_)
print(f"P(X>2) = {1 - exp_dist.cdf(2):.4f}")

4. 多维随机变量与联合分布分析

处理二维随机变量问题时，Python的矩阵运算能力可以大显身手。让我们看一个联合分布的例子：

问题：设二维随机变量(X,Y)的分布律为：

Y\X  1    2    3
0   0.2  0.1  0.1
-1  0.15 0.2  0.25

求：(1)边缘分布 (2)判断独立性 (3)Z=X+Y的分布

# 构建联合分布表
joint_prob = np.array([[0.2, 0.1, 0.1],
                      [0.15, 0.2, 0.25]])

# (1) 边缘分布
x_marginal = np.sum(joint_prob, axis=0)
y_marginal = np.sum(joint_prob, axis=1)

print("X的边缘分布:")
for i, prob in enumerate(x_marginal, 1):
    print(f"P(X={i}) = {prob:.2f}")

print("\nY的边缘分布:")
for j, prob in zip([0, -1], y_marginal):
    print(f"P(Y={j}) = {prob:.2f}")

# (2) 判断独立性
independent = True
for i in range(joint_prob.shape[0]):
    for j in range(joint_prob.shape[1]):
        if not np.isclose(joint_prob[i,j], y_marginal[i]*x_marginal[j]):
            independent = False
            break
print(f"\nX与Y是否独立: {'是' if independent else '否'}")

# (3) Z = X + Y 的分布
z_values = []
z_probs = []

x = [1, 2, 3]
y = [0, -1]

for xi in x:
    for yi in y:
        z = xi + yi
        prob = joint_prob[y.index(yi), x.index(xi)]
        if z in z_values:
            idx = z_values.index(z)
            z_probs[idx] += prob
        else:
            z_values.append(z)
            z_probs.append(prob)

print("\nZ的分布律:")
for z, p in sorted(zip(z_values, z_probs)):
    print(f"P(Z={z}) = {p:.2f}")

对于连续型二维随机变量，我们可以类似地处理：

# 二维均匀分布示例
def joint_pdf(x, y):
    if 0 < x < 2 and 2 < y < 4:
        return (6 - x - y)/8
    return 0

# (1) 边缘密度p_X(x)
def marginal_px(x):
    if 0 < x < 2:
        result, _ = quad(lambda y: joint_pdf(x, y), 2, 4)
        return result
    return 0

# (2) P(X+Y ≤ 4)
# 积分区域: 2 < y < 4, 0 < x < 4-y
def integrand(x, y):
    return joint_pdf(x, y)

from scipy.integrate import dblquad
prob, _ = dblquad(integrand, 2, 4, lambda y: 0, lambda y: 4-y)
print(f"P(X+Y≤4) = {prob:.4f}")

5. 参数估计与假设检验

试卷中常见的参数估计问题，Python也能轻松应对：

问题：设总体X的密度函数为：

p(x) = { θx^(θ-1), 0 < x < 1
       { 0, 其他

X1,...,Xn为样本，求θ的矩估计和极大似然估计

# 生成模拟数据
theta_true = 2.5
np.random.seed(42)
sample = np.random.power(theta_true, size=100)

# (1) 矩估计
sample_mean = np.mean(sample)
theta_moment = sample_mean / (1 - sample_mean)
print(f"矩估计值: {theta_moment:.4f}")

# (2) 极大似然估计
def neg_log_likelihood(theta, data):
    return -np.sum(np.log(theta) + (theta-1)*np.log(data))

from scipy.optimize import minimize
res = minimize(neg_log_likelihood, x0=1, args=(sample,), bounds=[(0.1, None)])
theta_mle = res.x[0]
print(f"极大似然估计值: {theta_mle:.4f}")

# 与理论值比较
print(f"真实theta值: {theta_true}")

对于正态分布参数的估计与检验，SciPy提供了更便捷的方法：

# 正态分布参数估计
mu, sigma = 1, 2
data = stats.norm(mu, sigma).rvs(100)

# 参数估计
mu_hat, sigma_hat = stats.norm.fit(data)
print(f"μ估计: {mu_hat:.4f}, σ估计: {sigma_hat:.4f}")

# t检验示例
t_stat, p_value = stats.ttest_1samp(data, popmean=0)
print(f"t统计量: {t_stat:.4f}, p值: {p_value:.4f}")

6. 可视化验证与结果分析

可视化是理解概率分布的有力工具，Matplotlib可以帮我们直观验证计算结果：

import matplotlib.pyplot as plt

# 绘制正态分布密度函数
mu, sigma = 1, 2
x = np.linspace(mu-4*sigma, mu+4*sigma, 100)
pdf = stats.norm.pdf(x, mu, sigma)

plt.figure(figsize=(10, 5))
plt.plot(x, pdf, label=f'N({mu},{sigma**2})')
plt.title('正态分布概率密度函数')
plt.xlabel('x')
plt.ylabel('概率密度')
plt.legend()
plt.grid(True)
plt.show()

# 绘制二项分布与泊松分布比较
n, p = 100, 0.02
lambda_ = n * p

x = np.arange(0, 10)
binom_pmf = stats.binom.pmf(x, n, p)
poisson_pmf = stats.poisson.pmf(x, lambda_)

plt.figure(figsize=(10, 5))
plt.bar(x-0.1, binom_pmf, width=0.2, label='二项分布')
plt.bar(x+0.1, poisson_pmf, width=0.2, label='泊松分布')
plt.title('二项分布与泊松分布比较(n=100,p=0.02)')
plt.xlabel('k')
plt.ylabel('P(X=k)')
plt.legend()
plt.grid(True)
plt.show()

对于二维随机变量，我们可以绘制联合分布和边缘分布：

from mpl_toolkits.mplot3d import Axes3D

# 二维正态分布示例
mu = [0, 1]
cov = [[1, 0.5], [0.5, 2]]
x, y = np.random.multivariate_normal(mu, cov, 5000).T

# 绘制散点图
plt.figure(figsize=(10, 5))
plt.scatter(x, y, alpha=0.1)
plt.title('二维正态分布样本')
plt.xlabel('X')
plt.ylabel('Y')
plt.grid(True)
plt.show()

# 绘制联合密度函数
fig = plt.figure(figsize=(12, 5))
ax = fig.add_subplot(121, projection='3d')

X, Y = np.mgrid[-3:3:.1, -3:4:.1]
pos = np.dstack((X, Y))
rv = stats.multivariate_normal(mu, cov)
Z = rv.pdf(pos)

ax.plot_surface(X, Y, Z, cmap='viridis')
ax.set_title('联合概率密度函数')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_zlabel('密度')

# 绘制边缘分布
ax2 = fig.add_subplot(122)
ax2.plot(X[0,:], stats.norm.pdf(X[0,:], mu[0], np.sqrt(cov[0][0])), label='X边缘分布')
ax2.plot(Y[:,0], stats.norm.pdf(Y[:,0], mu[1], np.sqrt(cov[1][1])), label='Y边缘分布')
ax2.set_title('边缘分布')
ax2.legend()
ax2.grid(True)

plt.tight_layout()
plt.show()

7. 实战技巧与性能优化

在处理大规模概率计算时，效率变得尤为重要。以下是几个提升计算效率的技巧：

向量化计算：利用NumPy的向量化运算替代循环

# 低效方式
n = 1000000
slow_result = np.zeros(n)
for i in range(n):
    slow_result[i] = stats.norm.pdf(i/100000)

# 高效向量化方式
x = np.linspace(0, 10, n)
fast_result = stats.norm.pdf(x)

并行计算：对于独立重复实验，可以使用并行加速

from concurrent.futures import ProcessPoolExecutor

def simulate(n):
    return np.mean(stats.norm.rvs(size=n))

n_sim = 1000
sample_size = 1000

# 串行计算
%time serial = [simulate(sample_size) for _ in range(n_sim)]

# 并行计算
with ProcessPoolExecutor() as executor:
    %time parallel = list(executor.map(simulate, [sample_size]*n_sim))

内存优化：对于大型数组，注意内存管理

# 创建大型数组时使用更高效的数据类型
large_array = np.zeros(int(1e7), dtype=np.float32)  # 比默认float64节省一半内存

# 及时释放不再需要的大数组
del large_array

随机数生成优化：

# 一次性生成大量随机数比多次生成小批量更高效
%time slow = [stats.norm.rvs() for _ in range(10000)]
%time fast = stats.norm.rvs(size=10000)

# 使用更快的随机数生成器
from numpy.random import Generator, PCG64
rng = Generator(PCG64())
%time fastest = rng.normal(size=10000)

概率计算缓存：对于需要重复计算的分布，可以预先计算并存储

from functools import lru_cache

@lru_cache(maxsize=None)
def cached_pdf(x, mu=0, sigma=1):
    return stats.norm.pdf(x, mu, sigma)

# 第一次计算会实际执行
%time first = cached_pdf(1.5)
# 第二次相同参数调用会直接返回缓存结果
%time cached = cached_pdf(1.5)