I have a set of points where
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")
{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)