Fitting ellipse with 95% CI to points and calculating area within specific line segments

Question

I have a set of points where

• x is in [-1, 1],
• y is between [-1, 1], and
• the maximum sum between abs(x) + abs(y) = 1.

I want to fit an ellipsis (95% CI) and then calculate the are of the ellipsis that lies inside the segments A[(0,0);(1,0);(0.5,0.5)] and B[(0,0);(0,1);(0.5,0.5)]. How can I implement this?

Reproducible example

library(dplyr)
# Set seed for reproducibility
set.seed(30)
#Generate  points
n_points <- 1000
x <- runif(n=n_points, -1, 1)
y <- runif(n=n_points, -1, 1)
df<-cbind(x,y)%>%as.data.frame
df%>%
  mutate(sum=abs(x)+abs(x),
         ratio_x=x/abs(y),
         ratio_y=y/abs(x))%>%
  #Create skewed data
  filter(sum<=1,
         ratio_x>0,
         ratio_y<1)->df
# Visualize
plot(x~y,data=df)
#
# Fit an ellipse
df %>%
  select(x,y) %>%
  as.matrix -> matDat
#
matCovLS <- cov(matDat)
vecMeans <- colMeans(matDat)
vecMeans <- colMeans(matDat)
### get 95% CI ellipse
d2.95 <- qchisq(0.95, df = 2)
cluster::ellipsoidPoints(matCovLS, d2.95, loc = vecMeans) -> matEllipseHull95
plot(matDat, asp = 1, xlim = c(-1, 1))
lines(matEllipseHull95, col="blue")

Answer 1

{sf} is written for operations like this. In this approach, we use a variety of functions: st_cast(), st_multipoint(), st_multilinestring(), st_intersection(), st_union(), and st_area(). There are certainly more elegant ways.

library(sf)
elli <- st_cast(st_multipoint(matEllipseHull95), "POLYGON")
tri <- st_cast(st_multilinestring(
    list(matrix(c(0, 0, 1, 0, 0, 1, 0, 0) ,, 2 , byrow = TRUE))), "POLYGON")
plot(tri)
plot(elli, add = TRUE)
plot(st_intersection(st_union(elli), st_union(tri)), add = TRUE, col = "red")

which gives

st_area(st_intersection(st_union(elli), st_union(tri)))
[1] 0.2001515

A more complete plot:

library(scales)
plot(matDat, asp = 1, ylim = c(-1.5, 1), pch = 3)
plot(tri, add = TRUE)
plot(elli, border = "blue", add = TRUE)
plot(st_intersection(st_union(elli), st_union(tri)), add = TRUE, col = alpha("red", .25))
segments(x0 = 0 , y0 = 0, x1 = .5, y1 = .5, lty = 2)