Nothing
# Essentially, this file contains three functions (two exported, one internal).
# confidenceEllipse and bvNormalContour both take the eigen decomposition of a (covariance) matrix to plot
# the appropriate ellipse. These functions calculate the appropriate lengths of the major and minor axes
# and the angle of rotation of the ellipse; these values are then passed to eigenEllipseHelper to actually draw
# the ellipse.
# To do:
# make an add=TRUE parameter (at least for bvNormalContour so that one can make a contour plot with a loop)
# This might be easiest to just do by accepting a vector for alpha. Then, when it comes time to calculate lengths,
# check the length of clevel and if > 1 just loop through the tail end of the function.
#
# change X.mean for confidenceEllipse, allow Sigma = for confidence ellipse
#
# perhaps add options to graph the length of the axes (labels). something like c/lambda_1 = value, near the corresponding axis
# This functions graphs a confidence ellipse for mu based on the eigenvalues and eigenvectors for the covariance matrix S.
confidenceEllipse <- function(X.mean = c(0,0),
eig,
n,
p,
xl = NULL, yl = NULL, # the x and y axis limits; calculated dynamically if not specified
axes = TRUE, # if TRUE, the major and minor axes of the ellipse are graphed
center = FALSE, # if TRUE, a dot at the center of the ellipse and dashed lines to the axes are shown
lim.adj = 0.02, # an axis adjustment factor used if xl or yl is NULL (to adjust the calculated limits)
alpha = 0.05, # the alpha-level used to determine which ellipse is drawn
...){
# eigenEllipseHelper is expecting a matrix... might be better to eventually change that
X.mean <- matrix(X.mean, ncol=1)
# Calculate the angle of rotation.
# Because we are taking this as output of eigen(),
# the eigenvector in column 1 will correspond with the largest eigenvalue,
# which will thus be the eigenvector associated with the major axis,
# which is what is used for determining the angle of rotation.
# Actual calculation is just application of SOHCAHTOA to this vector.
# The result of this calculation is in radians.
angle <- atan(eig$vectors[2,1] / eig$vectors[1,1])
# Calculate the half-length of the major axis
# This is an application of result (5-19) in Johnson & Wichern (2007)
axis1 <- sqrt(eig$values[1]) *
sqrt((p * (n - 1)) /
(n * (n - p)) *
qf(1 - alpha, p, n - p))
# Calculate the half-length of the minor axis
# This is an application of result (5-19) in Johnson & Wichern (2007)
axis2 <- sqrt(eig$values[2]) *
sqrt((p * (n - 1)) /
(n * (n - p)) *
qf(1 - alpha, p, n - p))
# create a vector of axis lengths we can use
lengths <- c(axis1,axis2)
eigenEllipseHelper(mu = X.mean, lengths = lengths, angle = angle,
xl = xl, yl = yl,
axes = axes, center = center,
lim.adj = lim.adj, ...)
}
bvNormalContour <- function(mu = c(0,0), Sigma=NULL, eig=NULL,
xl = NULL, yl = NULL,
axes = TRUE, center = FALSE,
lim.adj = 0.02, alpha = 0.05, ...){
# eigenEllipseHelper is expecting a matrix... might be better to eventually change that
mu <- matrix(mu, ncol=1)
# Critical value for the constant density contour (always df=2 because this is bivariate normal)
# Johnson & Wichern (2008) result (4-8)
clevel <- qchisq(1 - alpha, df = 2)
# User needs to supply either eig or Sigma
if (!is.null(Sigma)){
eig <- eigen(Sigma)
}
# From Johnson & Wichern (2008) result (4-8)
# The half-lengths are c*sqrt(lambda_i); the eigenvectors determine the angle
#lengths <- c(clevel * sqrt(s11+s12),clevel*sqrt(s11-s12))
lengths <- c(clevel * sqrt(eig$values[1]), clevel * sqrt(eig$values[2]))
angle <- atan(eig$vectors[2,1]/eig$vectors[1,1]) # sohcahtoa
# Given the center, axis lengths, and angle of rotation, draw the ellipse.
eigenEllipseHelper(mu = mu, lengths = lengths, angle = angle,
xl = xl, yl = yl,
axes = axes, center = center,
lim.adj = lim.adj, ...)
}
# There seems to be a problem with the lim.adj for the xl yl calculation
eigenEllipseHelper <- function(mu, lengths, angle, xl, yl, lim.adj, axes, center, ...){
axis1 <- lengths[1]
axis2 <- lengths[2]
# Calculate x and y components for each axis.
# These are used if xl or y1 is NULL or if center = TRUE
axis1.x <- cos(angle) * axis1
axis1.y <- sin(angle) * axis1
axis2.x <- cos(angle + pi / 2) * axis2
axis2.y <- sin(angle + pi / 2) * axis2
# if no x limits are specified, calculate some that will ensure the entirety of the ellipse will fit on the x axis
if (is.null(xl)){
xl1 <- mu[1,1] - abs(axis1.x) - abs(axis2.x)
xl2 <- mu[1,1] + abs(axis1.x) + abs(axis2.x)
xl <- c(xl1*(1-lim.adj),xl2*(1+lim.adj))
}
# if no y limits are specified, calculate some that will ensure the entirety of the ellipse will fit on the y axis
if (is.null(yl)){
yl1 <- mu[2,1] - abs(axis1.y) - abs(axis2.y)
yl2 <- mu[2,1] + abs(axis1.y) + abs(axis2.y)
yl <- c(yl1*(1-lim.adj),yl2*(1+lim.adj))
}
# Call an empty graph that we can add the normal contour ellipseto
plot(0,pch='',ylab='',xlab='',xlim=xl,ylim=yl)
# draw the ellipse
# draw.ellipse is from package plotrix
plotrix::draw.ellipse(x = mu[1,1], y = mu[2,1], # center of the ellipse
a = lengths[1], b = lengths[2], # major and minor axis half-lengths
angle = angle, deg = FALSE) # calculated angle earlier was in radians
# if TRUE draw the major and minor axes for the normal contour ellipse
# These are just line segments using the offsets from the mean calculcated earlier.
if (axes){
segments(x1 = mu[1,1], y1 = mu[2,1], x0 = mu[1,1] + axis1.x, y0 = mu[2,1] + axis1.y)
segments(x1 = mu[1,1], y1 = mu[2,1], x0 = mu[1,1] - axis1.x, y0 = mu[2,1] - axis1.y)
segments(x1 = mu[1,1], y1 = mu[2,1], x0 = mu[1,1] + axis2.x, y0 = mu[2,1] + axis2.y)
segments(x1 = mu[1,1], y1 = mu[2,1], x0 = mu[1,1] - axis2.x, y0 = mu[2,1] - axis2.y)
}
# if TRUE show the center as a point and draw dashed lines orthogonal to the x and y axes
# These are just line vertical and horizontal segments from the origin to the x and y axes, respectively
if (center){
segments(x1 = mu[1,1], y1 = mu[2,1], y0 = mu[2,1], x0 = 0, lty=2)
segments(x1 = mu[1,1], y1 = mu[2,1], x0 = mu[1,1], y0 = 0, lty=2)
points(x = mu[1,1], y = mu[2,1], pch = 19)
}
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.