我想我有一个需要两个功能的解决方案,cov.rob
来自MASS
包装和ellipsoidhull
来自cluster
包裹。cov.rob(xy, quantile.used = 50, method = "mve")
从 2d 点总数中找到大约“最佳”50 个点xy
包含在最小体积椭圆中。然而,cov.rob
不直接返回this椭圆,而是从最佳点估计的其他一些椭圆(目标是稳健地估计协方差矩阵)。为了找到实际的最小椭圆,我们可以给出最好的点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")
看起来不错!您还可以检查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
这是一个方便的函数,可以将椭圆添加到现有的基本图形中:
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)
}
让我们在更多的点上使用它:
x <- runif(100)
y <- runif(100)
xy <- cbind(x, y)
plot(xy)
plot_min_ellipse(xy, points_in_ellipse = 50)