如何在Matlab中绘制以加权均值为中心、含95%数据质量的椭圆

问题背景 #

我使用Matlab 2022a处理空间数据，同时也接受通用公式/算法或R语言的解决方案。现有一组包含x、y坐标及对应位置数据值（第三列）的模拟数据，生成代码如下：

num = 100;
P = mvnrnd([0,0,20], [ 11.3669 -2.9268 -2.0568;
-2.9268 1.3372 0.4483;
-2.0568 0.4483 51.9697],num);
CenterofMass=sum(P(:,1:2).*P(:,3))/sum(P(:,3)); % 加权均值（质心）
scatter(P(:,1),P(:,2),abs(P(:,3))), hold on
plot(CenterofMass(:,1),CenterofMass(:,2),'r+')

需要完成的任务：

• 计算并绘制以质心（加权均值）为中心的椭圆，使椭圆的轴能包含指定统计比例（如95%）的数据质量（对应第三列数据值）。

此前我仅找到无权重的椭圆绘制方案，在处理带权重的空间协方差时遇到困难。

通用核心算法 #

生成带权重的椭圆，核心是计算加权协方差矩阵，再结合卡方分布确定椭圆缩放因子，步骤如下：

1. 计算加权质心：公式为 $ \boldsymbol{\mu} = \frac{\sum w_i \mathbf{x}_i}{\sum w_i} $，其中$w_i$是第三列数据值，$\mathbf{x}_i=(x_i,y_i)$为坐标点。
2. 计算加权协方差矩阵：
   $$
   \boldsymbol{\Sigma} = \frac{1}{\sum w_i} \sum w_i (\mathbf{x}_i - \boldsymbol{\mu})(\mathbf{x}_i - \boldsymbol{\mu})^T
   $$
   若需无偏样本协方差，可将分母改为$\sum w_i - 1$；针对质量比例场景，用总质量作为分母更贴合需求。
3. 确定椭圆缩放因子：二维场景下，要包含$p$比例的总质量，需获取卡方分布的$p$分位数$\chi^{2_{2,p}$（如95%对应$\chi}2_{2,0.95}=5.991$）。椭圆半轴长度为$\sqrt{\lambda_i \cdot \chi^2_{2,p}}$，其中$\lambda_i$是协方差矩阵的特征值，特征向量对应椭圆的轴方向。
4. 生成椭圆轮廓：通过极坐标参数化，结合特征向量旋转生成椭圆的轮廓点。

Matlab 2022a 实现代码 #

num = 100;
P = mvnrnd([0,0,20], [ 11.3669 -2.9268 -2.0568;
-2.9268 1.3372 0.4483;
-2.0568 0.4483 51.9697],num);
w = P(:,3); % 提取第三列作为权重（数据质量）
total_w = sum(w);

% 1. 计算加权质心
CenterofMass = sum(P(:,1:2).*w)/total_w;

% 2. 计算加权协方差矩阵
centered_xy = P(:,1:2) - CenterofMass;
weighted_cov = (centered_xy' .* w) * centered_xy / total_w;

% 3. 提取特征值与特征向量（确定椭圆的轴长和方向）
[V, D] = eig(weighted_cov);
lambda = diag(D); % 特征值对应半轴平方的基础值

% 4. 设置质量包含比例，获取对应卡方分位数
conf_level = 0.95;
chi2_val = chi2inv(conf_level, 2); % 二维卡方分布分位数

% 5. 生成椭圆轮廓点
theta = linspace(0, 2*pi, 100);
ellipse_points = CenterofMass + V * diag(sqrt(lambda * chi2_val)) * [cos(theta); sin(theta)];

% 可视化
scatter(P(:,1), P(:,2), abs(w)), hold on
plot(CenterofMass(1), CenterofMass(2), 'r+', 'MarkerSize', 10)
plot(ellipse_points(1,:), ellipse_points(2,:), 'b-', 'LineWidth', 1.5)
hold off
xlabel('X坐标')
ylabel('Y坐标')
title(['包含', num2str(conf_level*100), '%数据质量的加权椭圆'])

R语言实现代码 #

library(MASS)
library(ggplot2)

# 生成模拟数据
num <- 100
cov_mat <- matrix(c(11.3669, -2.9268, -2.0568,
                    -2.9268, 1.3372, 0.4483,
                    -2.0568, 0.4483, 51.9697), nrow=3)
P <- mvrnorm(num, mu=c(0,0,20), Sigma=cov_mat)
w <- P[,3]
total_w <- sum(w)

# 1. 计算加权质心
center_of_mass <- colSums(P[,1:2] * w) / total_w

# 2. 计算加权协方差矩阵
centered_xy <- sweep(P[,1:2], 2, center_of_mass)
weighted_cov <- t(centered_xy * w) %*% centered_xy / total_w

# 3. 提取特征值与特征向量
eig <- eigen(weighted_cov)
lambda <- eig$values
V <- eig$vectors

# 4. 获取卡方分位数
conf_level <- 0.95
chi2_val <- qchisq(conf_level, df=2)

# 5. 生成椭圆轮廓点
theta <- seq(0, 2*pi, length.out=100)
ellipse_points <- center_of_mass + V %*% diag(sqrt(lambda * chi2_val)) %*% rbind(cos(theta), sin(theta))
ellipse_df <- data.frame(x=ellipse_points[1,], y=ellipse_points[2,])

# 可视化
ggplot() +
  geom_point(aes(x=P[,1], y=P[,2], size=abs(w)), alpha=0.6) +
  geom_point(aes(x=center_of_mass[1], y=center_of_mass[2]), color="red", shape="+", size=5) +
  geom_path(data=ellipse_df, aes(x=x, y=y), color="blue", linewidth=1) +
  labs(x="X坐标", y="Y坐标", size="数据值", title=paste0("包含", conf_level*100, "%数据质量的加权椭圆")) +
  theme_bw()

内容的提问来源于stack exchange，提问作者Rebecca Thomas - NOAA Federal