/
Ellipsoid.R
159 lines (137 loc) · 5.66 KB
/
Ellipsoid.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
# was an internal function in heplot3d
# adapted from the shapes3d demo in the rgl package and from the Rcmdr package
# modified to return the bbox of the ellipsoid
# TODO:
# Handle greater, less than 3D
# Use label.ellipse?
#' @name Ellipsoid
#' @aliases Ellipsoid Ellipsoid.data.frame Ellipsoid.default
#'
#' @title Draw an Ellipsoid in an rgl Scene
#'
#' @description
#' This is an experimental function designed to separate internal code in \code{link{heplot3d}}.
#'
#' @param x An object. In the default method the parameter x should be a square positive definite matrix at least 3x3 in size. It will be treated as the correlation or covariance of a multivariate normal
#' distribution. For the \code{data.frame} method, it should be a numeric data frame with at
#' least 3 columns.
#' @param center center of the ellipsoid, a vector of length 3, typically the mean vector of data
#' @param which This parameter selects which variables from the object will be plotted. The default is the first 3.
#' @param method For the \code{data.frame} method, a character string to specify the covariance method to
#' be used: \emph{classical} product-moment (\code{"classical"}), or \emph{minimum volume ellipsoid}
#' (\code{"mve"}), or \emph{minimum covariance determinant} (\code{"mcd"}).
#' @param radius size of the ellipsoid
#' @param df degrees of freedom associated with the covariance matrix, used to calculate the appropriate F statistic
#' @param label label for the ellipsoid
#' @param cex.label text size of label
#' @param col color of the ellipsoid
#' @param lwd line with for the wire-frame version
#' @param segments number of segments composing each ellipsoid; defaults to \code{40}.
#' @param shade logical; should the ellipsoid be smoothly shaded?
#' @param alpha transparency of the shaded ellipsoid
#' @param wire logical; should the ellipsoid be drawn as a wire frame?
#' @param verbose logical; for debugging
#' @param warn.rank logical; warn if the ellipsoid is less than rank 3?
#' @param ... Other arguments
#'
#' @return returns the bounding box of the ellipsoid invisibly; otherwise used for it's side effect of
#' drawing the ellipsoid
#' @export
#'
#' @examples
#' # none yet
Ellipsoid <-
function(x, ...) UseMethod("Ellipsoid")
#' @param method the covariance method to be used: classical product-moment (\code{"classical"}),
#' or minimum volume ellipsoid (\code{"mve"}), or
#' minimum covariance determinant (\code{"mcd"}
#' @rdname Ellipsoid
#' @exportS3Method Ellipsoid data.frame
Ellipsoid.data.frame <- function(
x,
which = 1:3,
method = c("classical", "mve", "mcd"),
...) {
method <- match.arg(method)
if (length(which) != 3L) stop("'which' must be a vector of length 3, not ", which)
if (any(which) > ncol(x) |
any(which) < 0) stop("unavailable variables selected in ", which)
x <- x[, which]
if (nrow(x) < 4) stop("at least 4 cases are needed")
rcov <- MASS::cov.rob(x, method=method)
cov <- rcov$cov
means<- rcov$center
Ellipsoid.default(cov, center = means, ...)
}
# Ellipsoid.matrix <- Ellipsoid.data.frame
#' @rdname Ellipsoid
#' @exportS3Method Ellipsoid default
Ellipsoid.default <- function(
x,
center = c(0,0,0),
which = 1:3,
radius = 1,
df = Inf,
label = "",
cex.label = 1.5,
col = "pink",
lwd = 1,
segments = 40, # line segments in each ellipse
shade = TRUE,
alpha = 0.1,
wire = TRUE,
verbose = FALSE,
warn.rank = FALSE,
...
){
degvec <- seq(0, 2*pi, length=segments)
ecoord2 <- function(p) c(cos(p[1])*sin(p[2]), sin(p[1])*sin(p[2]), cos(p[2]))
# v <- t(apply(expand.grid(degvec,degvec), 1, ecoord2)) # modified to make smoother
v <- t(apply(expand.grid(degvec,degvec/2), 1, ecoord2))
if (!warn.rank){
warn <- options(warn=-1)
on.exit(options(warn))
}
dimx <- dim(x)
if (!length(dimx) == 2 & dimx[1] == dimx[2]) stop("'x' must be a square matrix")
if (length(which) != 3) stop("'which' must be a vector of length 3, not ", which)
if (any(which) > dimx[1]) stop("unavailable variables selected in ", which)
x <- x[which, which]
shape <- x
# TODO: select which
Q <- chol(shape, pivot=TRUE)
lwd <- if (df < 2) lwd[2] else lwd[1]
order <- order(attr(Q, "pivot"))
v <- center + radius * t(v %*% Q[, order])
v <- rbind(v, rep(1,ncol(v)))
e <- expand.grid(1:(segments-1), 1:segments)
i1 <- apply(e, 1, function(z) z[1] + segments*(z[2] - 1))
i2 <- i1 + 1
i3 <- (i1 + segments - 1) %% segments^2 + 1
i4 <- (i2 + segments - 1) %% segments^2 + 1
i <- rbind(i1, i2, i4, i3)
x <- rgl::asEuclidean(t(v))
ellips <- rgl::qmesh3d(v, i)
# override settings for 1 df line
if (df<2) {
wire <- TRUE
shade <- FALSE
}
back <- if (df < 3) "culled" else "filled"
depth_mask <- if (alpha <.8) FALSE else TRUE
if (verbose) cat(paste("df=", df, "col:", col, " shade:", shade, " alpha:", alpha,
" wire:", wire, "back:", back, "depth_mask:", depth_mask,
sep=" "), "\n")
if(shade) rgl::shade3d(ellips,
col=col, alpha=alpha, lit=TRUE, back=back, depth_mask=depth_mask)
if(wire) rgl::wire3d(ellips,
col=col, size=lwd, lit=FALSE, back=back, depth_mask=depth_mask)
# get bounding box of object
bbox <- matrix(rgl::par3d("bbox"), nrow=2)
ranges <- apply(bbox, 2, diff)
if (!is.null(label) && label !="")
rgl::texts3d(x[which.max(x[,2]),] + offset*ranges, adj=0, texts=label, color=col, lit=FALSE)
rownames(bbox) <- c("min", "max")
colnames(bbox) <- names(center)
invisible(bbox)
}