如何找到覆盖R中一组点的给定部分的最小椭圆?

Ras*_*åth 8 r ellipse minimum

我想知道:是否有一些功能/巧妙的方法来找到覆盖R中一组2d点的给定部分的最小椭圆?随着最小的我的意思是面积最小的椭圆形.

澄清:如果点数很大,我可以使用近似正确的解决方案(因为我猜一个确切的解决方案必须尝试点的子集的所有组合)

这个问题可能听起来像是包含R中给定点百分比的椭圆问题的副本,但问题的表达形式是得到的答案不会产生最小的椭圆.例如,使用给予Ellipse的解决方案,其中包含R中给定点的百分比:

require(car)
x <- runif(6)
y <- runif(6)
dataEllipse(x,y, levels=0.5)
Run Code Online (Sandbox Code Playgroud)

得到的椭圆显然不是包含一半点的最小椭圆,我猜,这是一个覆盖左上角三个点的小椭圆.

在此输入图像描述

Ras*_*åth 4

我想我有一个解决方案,需要两个函数,cov.rob来自MASS包和ellipsoidhull来自cluster包。从包含在最小体积椭圆中cov.rob(xy, quantile.used = 50, method = "mve")的 2d 点总数中找到大约“最佳”50 个点。xy但是,cov.rob不会直接返回椭圆,而是返回从最佳点估计的其他一些椭圆(目标是稳健地估计协方差矩阵)。ellipsoidhull为了找到实际的最小椭圆,我们可以给出找到最小椭圆的最佳点,并且我们可以predict.ellipse使用它来获取定义椭圆外壳的路径的坐标。

我不是 100% 确定这个方法是最简单的和/或它 100% 有效(感觉应该可以避免使用的秒步骤,ellipsoidhull但我还没弄清楚如何。)。它似乎至少适用于我的玩具示例......

说得够多了,这是代码:

library(MASS)
library(cluster)

# Using the same six points as in the question
xy <- cbind(x, y)
# Finding the 3 points in the smallest ellipse (not finding 
# the actual ellipse though...)
fit <- cov.rob(xy, quantile.used = 3, method = "mve")
# Finding the minimum volume ellipse that contains these three points
best_ellipse <- ellipsoidhull( xy[fit$best,] )
plot(xy)
# The predict() function returns a 2d matrix defining the coordinates of
# the hull of the ellipse 
lines(predict(best_ellipse), col="blue")
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述

看起来不错!您还可以检查该ellipse对象以获取更多信息

best_ellipse
## 'ellipsoid' in 2 dimensions:
##  center = ( 0.36 0.65 ); squared ave.radius d^2 =  2 
##  and shape matrix =
##         x      y
## x 0.00042 0.0065
## y 0.00654 0.1229
##   hence, area  =  0.018 
Run Code Online (Sandbox Code Playgroud)

这是一个方便的函数,可以将椭圆添加到现有的基本图形中:

plot_min_ellipse <- function(xy, points_in_ellipse, color = "blue") {
  fit <- cov.rob(xy, quantile.used = points_in_ellipse, method = "mve")
  best_ellipse <- ellipsoidhull( xy[fit$best,] )
  lines(predict(best_ellipse), col=color)
}
Run Code Online (Sandbox Code Playgroud)

让我们在更多的点上使用它:

x <- runif(100)
y <- runif(100)
xy <- cbind(x, y)
plot(xy)
plot_min_ellipse(xy, points_in_ellipse = 50)
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述