How to find the smallest ellipse covering a given fraction of a set of points in R?

Question

I'm wondering: Is there some function/clever way to find the smallest ellipse covering a given fraction of a set of 2d points in R? With smallest I mean the ellipse with the smallest area.

Clarification: I'm fine with an approximately correct solution if the number of points are large (as I guess an exact solution would have to try all combinations of subsets of points)

This question might sound like a duplicate of the question Ellipse containing percentage of given points in R but the way that question is phrased the resulting answer does not result in the smallest ellipse. For example, using the solution given to Ellipse containing percentage of given points in R:

require(car)
x <- runif(6)
y <- runif(6)
dataEllipse(x,y, levels=0.5)

The resulting ellipse is clearly not the smallest ellipse containing half of the points, Which, I guess, would be a small ellipse covering the three points up in the top-left corner.

Answer 1

I think I have a solution which requires two functions, cov.rob from the MASS package and ellipsoidhull from the cluster package. cov.rob(xy, quantile.used = 50, method = "mve") finds approximately the "best" 50 points out of the total number of 2d points in xy that are contained in the minimum volume ellipse. However, cov.rob does not directly return this ellipse but rather some other ellipse estimated from the best points (the goal being to robustly estimate the covariance matrix). To find the actuall minimum ellipse we can give the best points to ellipsoidhull which finds the minimum ellipse, and we can use predict.ellipse to get out the coordinates of the path defining the hull of the elllipse.

I'm not 100% certain this method is the easiest and/or that it works 100% (It feels like it should be possible to avoid the seconds step of using ellipsoidhull but I havn't figured out how.). It seems to work on my toy example at least....

Enough talking, here is the code:

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")

Looks pretty good! You can also inspect the ellipse object for more info

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 

Here is a handy function which adds an ellipse to an existing base graphics plot:

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)
}

Let's use it on a larger number of points:

x <- runif(100)
y <- runif(100)
xy <- cbind(x, y)
plot(xy)
plot_min_ellipse(xy, points_in_ellipse = 50)