Skip to content

Commit bd306cf

Browse files
committed
pulled calculation out into separate function
1 parent 30a72a6 commit bd306cf

File tree

1 file changed

+36
-33
lines changed

1 file changed

+36
-33
lines changed

R/stat-ellipse.R

Lines changed: 36 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,6 @@ stat_ellipse <- function(mapping = NULL, data = NULL, geom = "path", position =
4848
StatEllipse$new(mapping = mapping, data = data, geom = geom, position = position, type = type, level = level, segments = segments, na.rm = na.rm, ...)
4949
}
5050

51-
5251
StatEllipse <- proto(Stat, {
5352
objname <- "ellipse"
5453

@@ -60,40 +59,44 @@ StatEllipse <- proto(Stat, {
6059
}
6160
calculate <- function(., data, scales, type = "t", level = 0.95, segments = 51, na.rm = FALSE, ...){
6261
data <- remove_missing(data, na.rm, vars = c("x","y"), name = "stat_ellipse", finite = TRUE)
62+
ellipse <- calculate_ellipse(data=data, vars= c("x","y"), type=type, level=level, segments=segments)
63+
return(ellipse)
64+
}
65+
})
6366

64-
dfn <- 2
65-
dfd <- length(data$x) - 1
67+
calculate_ellipse <- function(data, vars, type, level, segments){
68+
dfn <- 2
69+
dfd <- nrow(data) - 1
6670

67-
if (!type %in% c("t", "norm", "euclid")){
68-
message("Unrecognized ellipse type")
69-
ellipse <- rbind(as.numeric(c(NA, NA)))
70-
} else if (dfd < 3){
71-
message("Too few points to calculate an ellipse")
72-
ellipse <- rbind(as.numeric(c(NA, NA)))
71+
if (!type %in% c("t", "norm", "euclid")){
72+
message("Unrecognized ellipse type")
73+
ellipse <- rbind(as.numeric(c(NA, NA)))
74+
} else if (dfd < 3){
75+
message("Too few points to calculate an ellipse")
76+
ellipse <- rbind(as.numeric(c(NA, NA)))
77+
} else {
78+
if (type == "t"){
79+
v <- cov.trob(data[,vars])
80+
} else if (type == "norm"){
81+
v <- cov.wt(data[,vars])
82+
} else if (type == "euclid"){
83+
v <- cov.wt(data[,vars])
84+
v$cov <- diag(rep(min(diag(v$cov)), 2))
85+
}
86+
shape <- v$cov
87+
center <- v$center
88+
chol_decomp <- chol(shape)
89+
if (type == "euclid"){
90+
radius <- level/max(chol_decomp)
7391
} else {
74-
if (type == "t"){
75-
v <- cov.trob(cbind(data$x, data$y))
76-
} else if (type == "norm"){
77-
v <- cov.wt(cbind(data$x, data$y))
78-
} else if (type == "euclid"){
79-
v <- cov.wt(cbind(data$x, data$y))
80-
v$cov <- diag(rep(min(diag(v$cov)), 2))
81-
}
82-
shape <- v$cov
83-
center <- v$center
84-
chol_decomp <- chol(shape)
85-
if (type == "euclid"){
86-
radius <- level/max(chol_decomp)
87-
} else {
88-
radius <- sqrt(dfn * qf(level, dfn, dfd))
89-
}
90-
angles <- (0:segments) * 2 * pi/segments
91-
unit.circle <- cbind(cos(angles), sin(angles))
92-
ellipse <- t(center + radius * t(unit.circle %*% chol_decomp))
92+
radius <- sqrt(dfn * qf(level, dfn, dfd))
9393
}
94-
95-
ellipse <- as.data.frame(ellipse)
96-
colnames(ellipse) <- c("x", "y")
97-
return(ellipse)
94+
angles <- (0:segments) * 2 * pi/segments
95+
unit.circle <- cbind(cos(angles), sin(angles))
96+
ellipse <- t(center + radius * t(unit.circle %*% chol_decomp))
9897
}
99-
})
98+
99+
ellipse <- as.data.frame(ellipse)
100+
colnames(ellipse) <- vars
101+
return(ellipse)
102+
}

0 commit comments

Comments
 (0)