给点集拟合95%置信区间椭圆并计算特定线段内的椭圆面积

1. 修正点集生成代码 #

原代码存在笔误（sum=abs(x)+abs(x)应为abs(x)+abs(y)），且计算比值时会出现除以0的情况，以下是修正后的点集生成代码：

library(dplyr)
set.seed(30)
n_points <- 1000

# 生成初始均匀分布点
x <- runif(n = n_points, -1, 1)
y <- runif(n = n_points, -1, 1)
df <- cbind(x, y) %>% as.data.frame()

# 筛选符合条件的点：|x|+|y|≤1，x/y>0（同号），y/x<1（|y|<|x|）
df <- df %>%
  mutate(
    sum_abs = abs(x) + abs(y),
    # 避免除以0，用ifelse处理x或y为0的情况
    ratio_x = ifelse(y == 0, NA, x / abs(y)),
    ratio_y = ifelse(x == 0, NA, y / abs(x))
  ) %>%
  filter(
    sum_abs <= 1,
    ratio_x > 0,  # x和y同号
    ratio_y < 1   # |y| < |x|
  ) %>%
  select(x, y)  # 保留需要的列

# 可视化点集
plot(y ~ x, data = df, asp = 1, main = "符合条件的点集")

2. 拟合95%置信区间椭圆 #

基于样本的均值和协方差矩阵，结合卡方分布的95%临界值生成椭圆：

library(cluster)

# 转换为矩阵格式
matDat <- as.matrix(df)

# 计算样本均值和协方差
vecMeans <- colMeans(matDat)
matCovLS <- cov(matDat)

# 95%置信区间对应的卡方临界值（自由度为2）
d2.95 <- qchisq(0.95, df = 2)

# 生成椭圆的轮廓点
matEllipseHull95 <- cluster::ellipsoidPoints(matCovLS, d2.95, loc = vecMeans)

# 可视化点集和椭圆
plot(y ~ x, data = df, asp = 1, main = "点集与95%CI椭圆")
lines(matEllipseHull95, col = "blue", lwd = 2)

3. 定义目标区域（线段组对应的三角形） #

线段组A和B分别对应两个三角形区域，我们用sf包将其定义为空间多边形：

library(sf)

# 定义线段组A的三角形：(0,0)、(1,0)、(0.5,0.5)
poly_A <- st_polygon(list(rbind(
  c(0, 0), c(1, 0), c(0.5, 0.5), c(0, 0)
))) %>% st_sfc()

# 定义线段组B的三角形：(0,0)、(0,1)、(0.5,0.5)
poly_B <- st_polygon(list(rbind(
  c(0, 0), c(0, 1), c(0.5, 0.5), c(0, 0)
))) %>% st_sfc()

# 可视化目标区域
plot(y ~ x, data = df, asp = 1, main = "点集、椭圆与目标区域")
lines(matEllipseHull95, col = "blue", lwd = 2)
plot(poly_A, col = rgb(1, 0, 0, 0.3), add = TRUE)
plot(poly_B, col = rgb(0, 1, 0, 0.3), add = TRUE)

4. 计算椭圆在目标区域内的面积 #

将椭圆轮廓转换为多边形，然后与两个三角形求交集，最后计算交集的面积：

# 将椭圆轮廓转换为sf多边形
ellipse_poly <- st_polygon(list(rbind(matEllipseHull95, matEllipseHull95[1, ]))) %>% st_sfc()

# 计算椭圆与三角形A的交集面积
intersect_A <- st_intersection(ellipse_poly, poly_A)
area_A <- st_area(intersect_A)

# 计算椭圆与三角形B的交集面积
intersect_B <- st_intersection(ellipse_poly, poly_B)
area_B <- st_area(intersect_B)

# 输出结果
cat("椭圆在线段组A区域内的面积：", as.numeric(area_A), "\n")
cat("椭圆在线段组B区域内的面积：", as.numeric(area_B), "\n")

说明 #

• 椭圆拟合的原理：95%置信区间椭圆基于多元正态分布的性质，以样本均值为中心，协方差矩阵缩放qchisq(0.95, 2)倍（约5.991）。
• 空间计算依赖sf包，若未安装可执行：install.packages(c("dplyr", "cluster", "sf"))

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