如何用Python基于椭圆参数生成指定标准差与角度的二元正态分布

核心参数对应关系 #

先明确椭圆参数和二元正态分布（BNPD）的严格对应规则：

• 椭圆长轴长度 = 3 × 对应方向的标准差（σ₁）→ σ₁ = 长轴长度 / 3
• 椭圆短轴长度 = 3 × 对应方向的标准差（σ₂）→ σ₂ = 短轴长度 / 3
• 椭圆旋转角θ = BNPD的旋转角（即协方差矩阵特征向量的方向角）

正确构造协方差矩阵的步骤 #

旋转角偏差的根源是协方差矩阵构造错误，以下是标准实现流程：

1. 将角度转为弧度：numpy/scipy的三角函数默认使用弧度，直接传角度会导致计算错误
2. 构造逆时针旋转矩阵：确保旋转方向与椭圆旋转角一致
3. 生成对角方差矩阵：存储两个主轴方向的方差（标准差的平方）
4. 通过旋转矩阵得到最终协方差矩阵：旋转对角矩阵得到对应方向的协方差

完整可运行代码 #

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

# 椭圆参数（按需修改）
target_theta = 45  # 预期旋转角（°）
major_axis = 6     # 长轴长度 → σ₁ = 6/3 = 2
minor_axis = 3     # 短轴长度 → σ₂ = 3/3 = 1
mean = [0, 0]      # 分布均值（可自定义）

# 1. 构造协方差矩阵
theta_rad = np.radians(target_theta)
# 逆时针旋转矩阵
rot_matrix = np.array([
    [np.cos(theta_rad), -np.sin(theta_rad)],
    [np.sin(theta_rad), np.cos(theta_rad)]
])
# 对角方差矩阵（主轴方向的方差）
diag_cov = np.diag([(major_axis/3)**2, (minor_axis/3)**2])
# 最终协方差矩阵
cov_matrix = rot_matrix @ diag_cov @ rot_matrix.T

# 2. 创建二元正态分布
bnpd = multivariate_normal(mean=mean, cov=cov_matrix)

# 3. 验证旋转角是否符合预期
eigenvalues, eigenvectors = np.linalg.eig(cov_matrix)
# 找到对应长轴的特征向量（最大特征值对应长轴）
major_idx = np.argmax(eigenvalues)
major_dir = eigenvectors[:, major_idx]
# 计算实际旋转角（转成0-360°范围）
calculated_theta = np.degrees(np.arctan2(major_dir[1], major_dir[0])) % 360
print(f"预期旋转角: {target_theta}°, 实际旋转角: {calculated_theta:.1f}°")

# 4. 可视化验证（绘制3σ椭圆和样本点）
x, y = np.mgrid[-5:5:.01, -5:5:.01]
pos = np.dstack((x, y))
# 绘制3σ等概率线（概率密度为均值处的exp(-9/2)倍）
plt.contour(x, y, bnpd.pdf(pos), levels=[bnpd.pdf(mean) * np.exp(-9/2)])
# 生成样本点验证分布形态
samples = bnpd.rvs(size=1000)
plt.scatter(samples[:,0], samples[:,1], s=1, alpha=0.5)
plt.axis('equal')
plt.title(f"二元正态分布（旋转角{target_theta}°）")
plt.show()

常见错误排查 #

• 旋转矩阵方向错误：如果用顺时针旋转矩阵，会导致旋转角与预期相反，必须使用上述逆时针旋转矩阵
• 混淆标准差与方差：协方差矩阵存储的是方差值，必须用(轴长/3)²而非直接用轴长或标准差
• 角度未转弧度：直接传入角度值到三角函数，会导致旋转角计算严重偏差
• 长轴短轴对应错误：若颠倒σ₁和σ₂的位置，椭圆的长短轴会互换，但旋转角仍正确，需注意参数对应关系

内容的提问来源于stack exchange，提问作者Wade Wang