groupSpread <-
function(xy, center=FALSE, plots=TRUE, CEPlevel=0.5, CIlevel=0.95,
CEPtype="CorrNormal", bootCI="none", dstTarget, conversion) {
UseMethod("groupSpread")
}
groupSpread.data.frame <-
function(xy, center=FALSE, plots=TRUE, CEPlevel=0.5, CIlevel=0.95,
CEPtype="CorrNormal", bootCI="none", dstTarget, conversion) {
## distance to target from override or from data
if(missing(dstTarget)) {
dstTarget <- if(hasName(xy, "distance")) {
xy[["distance"]]
} else {
NA_real_
}
}
## determine conversion factor from data if override is not given
if(missing(conversion)) {
conversion <- determineConversion(xy)
}
xy <- getXYmat(xy, center=center)
center <- FALSE # centering was done in getXYmat()
groupSpread(xy, center=center, plots=plots, CEPlevel=CEPlevel,
CIlevel=CIlevel, CEPtype=CEPtype, bootCI=bootCI,
dstTarget=dstTarget, conversion=conversion)
# NextMethod("groupSpread")
}
groupSpread.default <-
function(xy, center=FALSE, plots=TRUE, CEPlevel=0.5, CIlevel=0.95,
CEPtype="CorrNormal", bootCI="none", dstTarget, conversion) {
if(!is.matrix(xy)) { stop("xy must be a matrix") }
if(!is.numeric(xy)) { stop("xy must be numeric") }
if(ncol(xy) != 2L) { stop("xy must have two columns") }
if(!is.numeric(CEPlevel)) { stop("CEPlevel must be numeric") }
if(CEPlevel <= 0) { stop("CEPlevel must be > 0") }
if(!is.numeric(CIlevel)) { stop("CIlevel must be numeric") }
if(CIlevel <= 0) { stop("CIlevel must be > 0") }
if(center) {
warning("Centering only works for data frames, ignored here")
}
bootCI <- match.arg(bootCI,
choices=c("none", "norm", "basic", "perc", "bca"),
several.ok=TRUE)
## check if CEP / CI level is given in percent
if(CEPlevel >= 1) {
while(CEPlevel >= 1) { CEPlevel <- CEPlevel / 100 }
warning(paste0("CEPlevel must be in (0,1) and was set to ", CEPlevel))
}
if(CIlevel >= 1) {
while(CIlevel >= 1) { CIlevel <- CIlevel / 100 }
warning(paste0("CIlevel must be in (0,1) and was set to ", CIlevel))
}
dstTarget <- if(missing(dstTarget) ||
all(is.na(dstTarget)) ||
(length(unique(dstTarget)) > 1L)) {
NA_real_
} else {
mean(dstTarget)
}
conversion <- if(missing(conversion) ||
all(is.na(conversion)) ||
(length(unique(conversion)) > 1L)) {
NA_character_
} else {
unique(conversion)
}
#####-----------------------------------------------------------------------
## prepare data
X <- xy[ , 1] # x-coords
Y <- xy[ , 2] # y-coords
Npts <- nrow(xy) # number of observations
res <- vector("list", 0) # empty list to later collect the results
## can we do robust estimation?
haveRobustbase <- requireNamespace("robustbase", quietly=TRUE)
haveRob <- if(haveRobustbase && (Npts >= 4L)) {
TRUE
} else {
if(Npts < 4L) {
warning("We need >= 4 points for robust estimations")
}
if(!haveRobustbase) {
warning("Please install package 'robustbase' for robust estimations")
}
FALSE
}
## to determine axis limits later, collect all results in a vector
axesCollX <- numeric(0)
axesCollY <- numeric(0)
#####-----------------------------------------------------------------------
## non-parametric bootstrap-CIs (basic and BCa)
## for all parameters where a CI is later reported
if(!("none" %in% bootCI)) { # do bootstrap CIs
NrplMin <- 1499 # minimum number of replications
Nrpl <- if("bca" %in% bootCI) { # number of replications
max(NrplMin, Npts+1) # BCa needs at least number of points
} else {
NrplMin
}
## sdX, sdY, sigma, RSD, MR for one replication
getSdSigRSDMR <- function(x, idx) {
sdXY <- sqrt(diag(cov(x[idx, ])))
rayParam <- getRayParam(x[idx, ], level=CIlevel, doRob=FALSE)
sigmaHat <- rayParam$sigma["sigma"]
RSDhat <- rayParam$RSD["RSD"]
MRhat <- rayParam$MR["MR"]
return(c(sdXY, sigmaHat, RSDhat, MRhat))
}
bs <- boot::boot(xy, statistic=getSdSigRSDMR, R=Nrpl) # bootstrap
## extract CIs for all returned parameters
sdXciBoot <- boot::boot.ci(bs, conf=CIlevel, type=bootCI, index=1)
sdYciBoot <- boot::boot.ci(bs, conf=CIlevel, type=bootCI, index=2)
sigCIboot <- boot::boot.ci(bs, conf=CIlevel, type=bootCI, index=3)
RSDciBoot <- boot::boot.ci(bs, conf=CIlevel, type=bootCI, index=4)
MRciBoot <- boot::boot.ci(bs, conf=CIlevel, type=bootCI, index=5)
## CI type names in output structure of boot.ci()
CInames <- c(basic="basic", norm="normal", perc="percent", bca="bca")
CItype <- CInames[bootCI]
sdXciBoot <- c(vapply(CItype, function(x) {
len <- length(sdXciBoot[[x]])
sdXciBoot[[x]][(len-1):len] }, numeric(2)))
sdYciBoot <- c(vapply(CItype, function(x) {
len <- length(sdYciBoot[[x]])
sdYciBoot[[x]][(len-1):len] }, numeric(2)))
sigCIboot <- c(vapply(CItype, function(x) {
len <- length(sigCIboot[[x]])
sigCIboot[[x]][(len-1):len] }, numeric(2)))
RSDciBoot <- c(vapply(CItype, function(x) {
len <- length(RSDciBoot[[x]])
RSDciBoot[[x]][(len-1):len] }, numeric(2)))
MRciBoot <- c(vapply(CItype, function(x) {
len <- length(MRciBoot[[x]])
MRciBoot[[x]][(len-1):len] }, numeric(2)))
names(sdXciBoot) <- paste("sdX", rep(bootCI, each=2), c("(", ")"))
names(sdYciBoot) <- paste("sdY", rep(bootCI, each=2), c("(", ")"))
names(sigCIboot) <- paste("sigma", rep(bootCI, each=2), c("(", ")"))
names(RSDciBoot) <- paste("RSD", rep(bootCI, each=2), c("(", ")"))
names(MRciBoot) <- paste("MR", rep(bootCI, each=2), c("(", ")"))
} else {
sdXciBoot <- numeric(0)
sdYciBoot <- numeric(0)
sigCIboot <- numeric(0)
RSDciBoot <- numeric(0)
MRciBoot <- numeric(0)
}
#####-----------------------------------------------------------------------
## standard deviations of x- and y-coords
varXY <- diag(cov(xy))
sdXY <- sqrt(varXY)
res$sdXY <- makeMOA(sdXY, dst=dstTarget, conversion=conversion)
## parametric CIs for true sd
alpha <- 1-CIlevel
sdXci <- sqrt((Npts-1)*varXY[1] / qchisq(c(1-(alpha/2), alpha/2), Npts-1))
sdYci <- sqrt((Npts-1)*varXY[2] / qchisq(c(1-(alpha/2), alpha/2), Npts-1))
names(sdXci) <- c("sdX (", "sdX )")
names(sdYci) <- c("sdY (", "sdY )")
## combine parametric and bootstrap CIs and add to results
res$sdXci <- makeMOA(c(sdXci, sdXciBoot), dst=dstTarget, conversion=conversion)
res$sdYci <- makeMOA(c(sdYci, sdYciBoot), dst=dstTarget, conversion=conversion)
## robust standard deviations of x- and y-coords
res$sdXYrob <- if(haveRob) {
rob <- robustbase::covMcd(xy, cor=FALSE)
sdXYrob <- sqrt(diag(rob$cov))
makeMOA(sdXYrob, dst=dstTarget, conversion=conversion)
} else {
NULL
} # if(haveRob)
## (robust) center and covariance-matrix
ctr <- colMeans(xy) # group center
ctrRob <- if(haveRob) {
rob$center
} else {
NULL
} # if(haveRob)
res$covXY <- cov(xy) # covariance matrix
res$covXYrob <- if(haveRob) {
rob$cov
} else {
NULL
} # if(haveRob)
#####-----------------------------------------------------------------------
## mean distance to group center and associated parameterss
dstCtr <- getDistToCtr(xy)
dstCtrMean <- mean(dstCtr)
dstCtrMed <- median(dstCtr)
dstCtrMax <- max(dstCtr)
## radial standard deviation
## http://ballistipedia.com/index.php?title=Describing_Precision
rayParam <- getRayParam(xy, level=CIlevel, doRob=FALSE)
## sigma, RSD, MR estimates with parametric confidence intervals
sigma <- rayParam$sigma
RSD <- rayParam$RSD
MR <- rayParam$MR
names(sigma) <- c("sigma", "sigma (", "sigma )")
names(RSD) <- c("RSD", "RSD (", "RSD )")
names(MR) <- c("MR", "MR (", "MR )")
mDTCsigRSDmr <- c(mean=dstCtrMean, median=dstCtrMed, max=dstCtrMax,
sigma["sigma"], RSD["RSD"], MR["MR"])
res$distToCtr <- makeMOA(mDTCsigRSDmr, dst=dstTarget, conversion=conversion)
## combine parametric and bootstrap CIs and add to results
sigmaCI <- c(sigma[c("sigma (", "sigma )")], sigCIboot)
res$sigmaCI <- makeMOA(sigmaCI, dst=dstTarget, conversion=conversion)
RSDci <- c(RSD[c("RSD (", "RSD )")], RSDciBoot)
res$RSDci <- makeMOA(RSDci, dst=dstTarget, conversion=conversion)
MRci <- c(MR[c("MR (", "MR )")], MRciBoot)
res$MRci <- makeMOA(MRci, dst=dstTarget, conversion=conversion)
#####-----------------------------------------------------------------------
## maximum pairwise distance
maxPD <- getMaxPairDist(xy)
res$maxPairDist <- makeMOA(maxPD$d, dst=dstTarget, conversion=conversion)
#####-----------------------------------------------------------------------
## width, height, FoM, diagonal of (minimum) bounding box
bb <- getBoundingBox(xy) # bounding box
bbMin <- getMinBBox(xy) # minimum bounding box
groupRect <- c(width=bb$width, height=bb$height, FoM=bb$FoM, diag=bb$diag)
groupRectMin <- c(width=bbMin$width, height=bbMin$height, FoM=bbMin$FoM, diag=bbMin$diag)
res$groupRect <- makeMOA(groupRect, dst=dstTarget, conversion=conversion)
res$groupRectMin <- makeMOA(groupRectMin, dst=dstTarget, conversion=conversion)
## for axis limits
axesCollX <- c(axesCollX, bb$pts[c(1, 3)], bbMin$pts[ , 1])
axesCollY <- c(axesCollY, bb$pts[c(2, 4)], bbMin$pts[ , 2])
#####-----------------------------------------------------------------------
## radius of minimum enclosing circle
minCirc <- getMinCircle(xy) # minimum enclosing circle
res$minCircleRad <- makeMOA(minCirc$rad, dst=dstTarget, conversion=conversion)
## for axis limits
axesCollX <- c(axesCollX, minCirc$ctr[1] + minCirc$rad,
minCirc$ctr[1] - minCirc$rad)
axesCollY <- c(axesCollY, minCirc$ctr[2] + minCirc$rad,
minCirc$ctr[2] - minCirc$rad)
#####-----------------------------------------------------------------------
## radii of minimum enclosing ellipse
minEll <- getMinEllipse(xy) # minimum enclosing circle
res$minEll <- t(rbind(semi_major=makeMOA(minEll$size[1], dst=dstTarget, conversion=conversion),
semi_minor=makeMOA(minEll$size[2], dst=dstTarget, conversion=conversion)))
#####-----------------------------------------------------------------------
## confidence ellipse measures
confEll <- getConfEll(xy, CEPlevel, dstTarget=dstTarget,
conversion=conversion, doRob=haveRob)
res$confEll <- confEll$size
## for axis limits
axesCollX <- c(axesCollX, confEll$ctr[1] + confEll$size["unit", "semi-major"],
confEll$ctr[1] - confEll$size["unit", "semi-major"])
axesCollY <- c(axesCollY, confEll$ctr[2] + confEll$size["unit", "semi-major"],
confEll$ctr[2] - confEll$size["unit", "semi-major"])
if(haveRob) {
res$confEllRob <- confEll$sizeRob
## for axis limits
axesCollX <- c(axesCollX, confEll$ctrRob[1] + confEll$sizeRob["unit", "semi-major"],
confEll$ctrRob[1] - confEll$sizeRob["unit", "semi-major"])
axesCollY <- c(axesCollY, confEll$ctrRob[2] + confEll$sizeRob["unit", "semi-major"],
confEll$ctrRob[2] - confEll$sizeRob["unit", "semi-major"])
} else {
res$confEllRob <- NULL
} # if(haveRob)
res$confEllShape <- confEll$shape
res$confEllShapeRob <- if(haveRob) {
confEll$shapeRob
} else {
NULL
} # if(haveRob)
#####-----------------------------------------------------------------------
## circular error probable
CEP <- getCEP(xy, CEPlevel=CEPlevel, type=CEPtype, dstTarget=dstTarget,
conversion=conversion, doRob=FALSE)
res$CEP <- CEP$CEP
#####-----------------------------------------------------------------------
## plotting
if(plots) {
## infer (x,y)-coord units from conversion
unitXY <- na.omit(getUnits(conversion, first=FALSE))
unitDst <- na.omit(getUnits(conversion, first=TRUE))
devNew <- getDevice() # platform-dependent window open
## distance to target may be heterogeneous
dstTargetPlot <- paste(unique(round(na.omit(dstTarget))), collapse=", ")
#####-------------------------------------------------------------------
## diagram: histogram for distances to group center
## Rayleigh fit and kernel density estimate
xRange <- range(dstCtr)
xPts <- seq(xRange[1], xRange[2], length.out=200)
yRayleigh <- dRayleigh(xPts, sigma["sigma"])
yKDE <- density(dstCtr)
devNew() # open new diagram
h <- hist(dstCtr, breaks="FD", plot=FALSE)
plot(h, freq=FALSE,
main="Histogram distances to center w/ kernel density estimate",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("distance [", unitXY, "]"),
ylim=c(0, max(c(h$density, yRayleigh, yKDE$y))))
rug(jitter(dstCtr)) # show single values
## add Rayleigh fit and kernel density estimate
lines(xPts, yRayleigh, lwd=2, col="blue")
lines(yKDE$x, yKDE$y, lwd=2, col="red")
legend(x="topright",
legend=c("Rayleigh distribution", "kernel density estimate"),
col=c("blue", "red"), lty=1, lwd=2, bg=rgb(1, 1, 1, 0.7))
#####-------------------------------------------------------------------
## diagram: 2D-scatter plot for the (x,y)-distribution
## determine axis limits
xLims <- range(c(X, axesCollX))
yLims <- range(c(Y, axesCollY))
devNew() # open new diagram
plot(Y ~ X, asp=1, xlim=xLims, ylim=yLims, pch=20,
main="Group (x,y)-coordinates",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("X [", unitXY, "]"), ylab=paste0("Y [", unitXY, "]"))
abline(v=0, h=0, col="lightgray") # add point of aim
## add group center and robust estimate for group center
points(ctr[1], ctr[2], col="red", pch=4, lwd=2, cex=1.5)
if(haveRob) {
points(ctrRob[1], ctrRob[2], col="blue", pch=4, lwd=2, cex=1.5)
} # if(haveRob)
## add confidence ellipses (parametric, robust),
## and a circle with mean distance to center
drawEllipse(confEll, pch=4, fg="red", lwd=2)
if(haveRob) {
drawEllipse(ctrRob, res$covXYrob, radius=confEll$magFac,
pch=4, fg="blue", lwd=2)
} # if(haveRob)
drawCircle(ctr, radius=mean(dstCtr), fg="black", lwd=2)
## add legend
legend(x="bottomleft", legend=c("center", "center (robust)",
paste(100*CEPlevel, "% confidence ellipse", sep=""),
paste(100*CEPlevel, "% confidence ellipse (robust)", sep=""),
"mean distance to center"),
col=c("red", "blue", "red", "blue", "black"),
pch=c(4, 4, NA, NA, NA),
lty=c(NA, NA, 1, 1, 1), lwd=2, bg=rgb(1, 1, 1, 0.7))
#####-------------------------------------------------------------------
## diagram: 2D-scatter plot for the (x,y)-distribution
devNew() # open new diagram
plot(Y ~ X, asp=1, xlim=xLims, ylim=yLims, pch=20,
main="Group (x,y)-coordinates",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("X [", unitXY, "]"), ylab=paste0("Y [", unitXY, "]"))
abline(v=0, h=0, col="lightgray") # add point of aim
points(ctr[1], ctr[2], col="gray40", pch=4, lwd=4, cex=2.5) # add group center
## add bounding box, minimum bounding box, minimum enclosing circle,
## and maximum group spread
drawBox(bb, fg="magenta", lwd=2)
drawBox2(bbMin, fg="red", lwd=2)
drawCircle(minCirc, fg="blue", lwd=2)
segments(x0=xy[maxPD$idx[1], 1], y0=xy[maxPD$idx[1], 2],
x1=xy[maxPD$idx[2], 1], y1=xy[maxPD$idx[2], 2],
col="green3", lwd=2)
## add legend
legend(x="bottomleft", legend=c("group center", "bounding box",
"minimum bounding box",
"minimum enclosing circle", "maximum group spread"),
col=c("gray40", "magenta", "red", "blue", "green3"),
pch=c(4, NA, NA, NA, NA), lty=c(NA, 1, 1, 1, 1), lwd=2,
bg=rgb(1, 1, 1, 0.7))
} # if(plots)
#####-----------------------------------------------------------------------
## return all the collected numerical results and tests
return(res)
}
groupSpreadPlot <-
function(xy, which=1L, center=FALSE, CEPlevel=0.5, CIlevel=0.95,
dstTarget, conversion) {
if(!is.data.frame(xy)) { stop("xy must be a data.frame") }
xy <- getXYmat(xy, center=center)
if(!is.numeric(CEPlevel)) { stop("CEPlevel must be numeric") }
if(CEPlevel <= 0) { stop("CEPlevel must be > 0") }
if(!is.numeric(CIlevel)) { stop("CIlevel must be numeric") }
if(CIlevel <= 0) { stop("CIlevel must be > 0") }
which <- match.arg(as.character(which), choices=1:4)
## check if CEP / CI level is given in percent
if(CIlevel >= 1) {
while(CIlevel >= 1) { CIlevel <- CIlevel / 100 }
warning(c("CIlevel must be in (0,1) and was set to ", CIlevel))
}
if(CEPlevel >= 1) {
while(CEPlevel >= 1) { CEPlevel <- CEPlevel / 100 }
warning(c("CEPlevel must be in (0,1) and was set to ", CEPlevel))
}
dstTarget <- if(missing(dstTarget) ||
all(is.na(dstTarget)) ||
(length(unique(dstTarget)) > 1L)) {
NA_real_
} else {
mean(dstTarget)
}
conversion <- if(missing(conversion) ||
all(is.na(conversion)) ||
(length(unique(conversion)) > 1L)) {
NA_character_
} else {
unique(conversion)
}
#####-----------------------------------------------------------------------
## prepare data
X <- xy[ , 1] # x-coords
Y <- xy[ , 2] # y-coords
Npts <- nrow(xy) # number of observations
res <- vector("list", 0) # empty list to later collect the results
## can we do robust estimation?
haveRobustbase <- requireNamespace("robustbase", quietly=TRUE)
haveRob <- if(haveRobustbase && (Npts >= 4L)) {
TRUE
} else {
if(Npts < 4L) {
warning("We need >= 4 points for robust estimations")
}
if(!haveRobustbase) {
warning("Please install package 'robustbase' for robust estimations")
}
FALSE
}
## to determine axis limits later, collect all results in a vector
axesCollX <- numeric(0)
axesCollY <- numeric(0)
#####-----------------------------------------------------------------------
## standard deviations of x- and y-coords
varXY <- diag(cov(xy))
sdXY <- sqrt(varXY)
## robust standard deviations of x- and y-coords
if(haveRob) {
rob <- robustbase::covMcd(xy, cor=FALSE)
sdXYrob <- sqrt(diag(rob$cov))
} # if(haveRob)
## (robust) center and covariance-matrix
ctr <- colMeans(xy) # group center
if(haveRob) {
ctrRob <- rob$center
} # if(haveRob)
res$covXY <- cov(xy) # covariance matrix
if(haveRob) {
res$covXYrob <- rob$cov
} # if(haveRob)
#####-----------------------------------------------------------------------
## mean distance to group center and associated parameterss
dstCtr <- getDistToCtr(xy)
## radial standard deviation
## http://ballistipedia.com/index.php?title=Describing_Precision
rayParam <- getRayParam(xy, level=CIlevel, doRob=FALSE)
## sigma, RSD, MR estimates with parametric confidence intervals
sigma <- rayParam$sigma
names(sigma) <- c("sigma", "sigma (", "sigma )")
#####-----------------------------------------------------------------------
## maximum pairwise distance
maxPD <- getMaxPairDist(xy)
res$maxPairDist <- makeMOA(maxPD$d, dst=dstTarget, conversion=conversion)
#####-----------------------------------------------------------------------
## width, height, FoM, diagonal of (minimum) bounding box
bb <- getBoundingBox(xy) # bounding box
bbMin <- getMinBBox(xy) # minimum bounding box
## for axis limits
axesCollX <- c(axesCollX, bb$pts[c(1, 3)], bbMin$pts[ , 1])
axesCollY <- c(axesCollY, bb$pts[c(2, 4)], bbMin$pts[ , 2])
#####-----------------------------------------------------------------------
## radius of minimum enclosing circle
minCirc <- getMinCircle(xy) # minimum enclosing circle
## for axis limits
axesCollX <- c(axesCollX, minCirc$ctr[1] + minCirc$rad,
minCirc$ctr[1] - minCirc$rad)
axesCollY <- c(axesCollY, minCirc$ctr[2] + minCirc$rad,
minCirc$ctr[2] - minCirc$rad)
#####-----------------------------------------------------------------------
## confidence ellipse measures
confEll <- getConfEll(xy, CEPlevel, dstTarget=dstTarget,
conversion=conversion, doRob=haveRob)
## for axis limits
axesCollX <- c(axesCollX, confEll$ctr[1] + confEll$size["unit", "semi-major"],
confEll$ctr[1] - confEll$size["unit", "semi-major"])
axesCollY <- c(axesCollY, confEll$ctr[2] + confEll$size["unit", "semi-major"],
confEll$ctr[2] - confEll$size["unit", "semi-major"])
#####-----------------------------------------------------------------------
## minimum enclosing ellipse measures
minEll <- getMinEllipse(xy)
## for axis limits
axesCollX <- c(axesCollX,
minEll$ctr[1] + minEll$size[1],
minEll$ctr[1] - minEll$size[1])
axesCollY <- c(axesCollY,
minEll$ctr[2] + minEll$size[1],
minEll$ctr[2] - minEll$size[1])
if(haveRob) {
res$confEllRob <- confEll$sizeRob
## for axis limits
axesCollX <- c(axesCollX, confEll$ctrRob[1] + confEll$sizeRob["unit", "semi-major"],
confEll$ctrRob[1] - confEll$sizeRob["unit", "semi-major"])
axesCollY <- c(axesCollY, confEll$ctrRob[2] + confEll$sizeRob["unit", "semi-major"],
confEll$ctrRob[2] - confEll$sizeRob["unit", "semi-major"])
} # if(haveRob)
#####-----------------------------------------------------------------------
## plotting
## infer (x,y)-coord units from conversion
unitXY <- na.omit(getUnits(conversion, first=FALSE))
unitDst <- na.omit(getUnits(conversion, first=TRUE))
## determine axis limits
xLims <- range(c(X, axesCollX))
yLims <- range(c(Y, axesCollY))
## distance to target may be heterogeneous
dstTargetPlot <- paste(unique(round(na.omit(dstTarget))), collapse=", ")
if(which == 1L) {
#####-------------------------------------------------------------------
## diagram: histogram for distances to group center
xRange <- range(dstCtr)
xPts <- seq(xRange[1], xRange[2], length.out=200)
yRayleigh <- dRayleigh(xPts, sigma["sigma"])
yKDE <- density(dstCtr)
h <- hist(dstCtr, breaks="FD", plot=FALSE)
plot(h, freq=FALSE,
main="Histogram distances to center w/ kernel density estimate",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("distance [", unitXY, "]"),
ylim=c(0, max(c(h$density, yRayleigh, yKDE$y))))
rug(jitter(dstCtr)) # show single values
## add Rayleigh fit and kernel density estimate
lines(xPts, yRayleigh, lwd=2, col="blue")
lines(yKDE$x, yKDE$y, lwd=2, col="red")
legend(x="topright",
legend=c("Rayleigh distribution", "kernel density estimate"),
col=c("blue", "red"), lty=1, lwd=2, bg=rgb(1, 1, 1, 0.7))
}
if(which == 2L) {
#####-------------------------------------------------------------------
## diagram: 2D-scatter plot for the (x,y)-distribution
plot(Y ~ X, asp=1, xlim=xLims, ylim=yLims, pch=20,
main="Group (x,y)-coordinates",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("X [", unitXY, "]"), ylab=paste0("Y [", unitXY, "]"))
abline(v=0, h=0, col="lightgray") # add point of aim
## add group center and robust estimate for group center
points(ctr[1], ctr[2], col="red", pch=4, lwd=2, cex=1.5)
if(haveRob) {
points(ctrRob[1], ctrRob[2], col="blue", pch=4, lwd=2, cex=1.5)
} # if(haveRob)
## add confidence ellipses (parametric, robust),
## and a circle with mean distance to center
drawEllipse(confEll, pch=4, fg="red", lwd=2)
if(haveRob) {
drawEllipse(ctrRob, res$covXYrob, radius=confEll$magFac,
pch=4, fg="blue", lwd=2)
} # if(haveRob)
drawCircle(ctr, radius=mean(dstCtr), fg="black", lwd=2)
## add legend
legend(x="bottomleft", legend=c("center", "center (robust)",
paste(100*CEPlevel, "% confidence ellipse", sep=""),
paste(100*CEPlevel, "% confidence ellipse (robust)", sep=""),
"mean distance to center"),
col=c("red", "blue", "red", "blue", "black"),
pch=c(4, 4, NA, NA, NA),
lty=c(NA, NA, 1, 1, 1), lwd=2, bg=rgb(1, 1, 1, 0.7))
}
if(which == 3L) {
#####-------------------------------------------------------------------
## diagram: 2D-scatter plot for the (x,y)-distribution
plot(Y ~ X, asp=1, xlim=xLims, ylim=yLims, pch=20,
main="Group (x,y)-coordinates",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("X [", unitXY, "]"), ylab=paste0("Y [", unitXY, "]"))
abline(v=0, h=0, col="lightgray") # add point of aim
points(ctr[1], ctr[2], col="gray40", pch=4, lwd=4, cex=2.5) # add group center
## add bounding box, minimum bounding box, and maximum group spread
drawBox(bb, fg="magenta", lwd=2)
drawBox2(bbMin, fg="red", lwd=2)
segments(x0=xy[maxPD$idx[1], 1], y0=xy[maxPD$idx[1], 2],
x1=xy[maxPD$idx[2], 1], y1=xy[maxPD$idx[2], 2],
col="blue", lwd=2)
## add legend
legend(x="bottomleft",
legend=c("group center",
"bounding box",
"minimum bounding box",
"maximum group spread"),
col=c("gray40", "magenta", "red", "blue"),
pch=c(4, NA, NA, NA),
lty=c(NA, 1, 1, 1),
lwd=2,
bg=rgb(1, 1, 1, 0.7))
} # if(plots)
if(which == 4L) {
#####-------------------------------------------------------------------
## diagram: 2D-scatter plot for the (x,y)-distribution
plot(Y ~ X, asp=1, xlim=xLims, ylim=yLims, pch=20,
main="Group (x,y)-coordinates",
sub=paste("distance:", dstTargetPlot, unitDst),
xlab=paste0("X [", unitXY, "]"), ylab=paste0("Y [", unitXY, "]"))
abline(v=0, h=0, col="lightgray") # add point of aim
points(ctr[1], ctr[2], col="gray40", pch=4, lwd=4, cex=2.5) # add group center
## add minimum enclosing circle, minimum enclosing ellipse,
## and maximum group spread
drawCircle(minCirc, fg="magenta", lwd=2)
drawEllipse(minEll, fg="red", lwd=2)
segments(x0=xy[maxPD$idx[1], 1], y0=xy[maxPD$idx[1], 2],
x1=xy[maxPD$idx[2], 1], y1=xy[maxPD$idx[2], 2],
col="blue", lwd=2)
## add legend
legend(x="bottomleft",
legend=c("group center",
"minimum enclosing circle",
"minimum enclosing ellipse",
"maximum group spread"),
col=c("gray40", "magenta", "red", "blue"),
pch=c(4, NA, NA, NA),
lty=c(NA, 1, 1, 1),
lwd=2,
bg=rgb(1, 1, 1, 0.7))
} # if(plots)
#####-----------------------------------------------------------------------
## return all the collected numerical results and tests
return(invisible(NULL))
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.