Nothing
rss2d <- function(design, lower, upper, gof.test.type="greenwood", gof.test.stat=NULL, transform=NULL, n.angle=360, graphics=1, trace=TRUE, lines.lwd=1, lines.lty="dotted", ...) {
design <- as.matrix(design)
n <- dim(design)[1]
d <- dim(design)[2]
# some arguments checks
if (!is.element(gof.test.type, c("greenwood", "qm", "ks", "V", "cvm"))) stop("The goodness-of-fit test must be one of: Greenwood, Quesenberry-Miller, Kolmogorov-Smirnov, V = (D+) + (D-), or Cramer-Von Mises")
if ((length(lower)!=d) & (length(upper)!=d)) stop("lower and upper must be d-dimensional vectors")
# domain transformation to [-1,1]^d
for (j in 1:d) {
design.col <- design[,j]
design.col.min <- min(design.col)
if (design.col.min < lower[j]) stop('the minimum of design values is not compatible with lower')
design.col.max <- max(design.col)
if (design.col.max > upper[j]) stop('the maximum of design values is not compatible with upper')
design.col <- 2*((design.col - lower[j])/(upper[j] - lower[j]) - 0.5)
design[,j] <- design.col
}
# compute the subdivision of the (half) circle in cartesian coordinates
theta.degree <- seq(from=0, to=180, length=n.angle+1)
theta.degree <- theta.degree[1:n.angle]
theta.degree <- as.matrix(theta.degree)
theta <- theta.degree*2*pi/360
cos.theta <- cos(theta)
sin.theta <- sin(theta)
n.theta <- length(theta)
subdiv.halfcircle <- cbind(cos.theta, sin.theta)
# loop over all pairs of dimensions
global.stat <- matrix(NA, d, d)
global.stat.max <- 0
if (trace) {
cat("\n2D Radial Scanning Statistic (RSS) with ", toupper(gof.test.type), " statistic\n", sep="")
cat("Discretization step (in degree) : ", 180/n.angle, sep="")
cat("\n\nMaximum of RS statistic values (global statistic) per pair of dimensions")
}
print.out <- data.frame(global.stat = rep(NA, (d*(d-1))%/%2))
meter <- 0
for (i in 1:(d-1)) {
x <- design[,i]
for (j in ((i+1):d)) {
y <- design[,j]
# compute anglewise statistic
# 1st step : compute the matrix of the F(projected points onto Vect(cos.theta, sin.theta))
out <- .C("C_rss2Dloop", as.double(x), as.double(y), as.double(cos.theta), as.double(sin.theta), as.integer(n), as.integer(n.theta), ans=double(n * n.theta), PACKAGE="DiceDesign")
F.projections <- matrix(out$ans, n, n.theta)
# 2nd step : for each angle, compute the statistic
# In future version, should be computed inside the C loop
anglewise.stat.ij <- matrix(NA, n.theta, 1)
for (angle in 1:n.theta) {
anglewise.stat.ij[angle] <- unif.test.statistic(x=F.projections[,angle], type=gof.test.type, transform=transform)
}
# compute the worst value over all angles and store it
global.stat.ij <- max(anglewise.stat.ij)
global.stat[i,j] <- global.stat[j,i] <- global.stat.ij
if (global.stat.ij > global.stat.max) {
global.stat.max <- global.stat.ij
pair.worst <- c(i,j)
anglewise.stat <- anglewise.stat.ij
}
meter <- meter + 1
print.out[meter, 1] <- global.stat.ij
name.current <- paste("(", i, ",", j, ")", sep="")
row.names(print.out)[meter] <- name.current
if (trace) cat("\n", name.current, " ", global.stat.ij, sep="")
} # end loop j
} # end loop i
if (trace) cat("\n\n")
# rss curve
rss.curve.x <- anglewise.stat*subdiv.halfcircle[,1]
rss.curve.y <- anglewise.stat*subdiv.halfcircle[,2]
# statistic upper tail percentage points
# see D'Agostino and Stephens "Goodness-of-fit techniques", 1986
if (is.null(gof.test.stat)) {
gof.test.stat <- unif.test.quantile(type=gof.test.type, n=n, alpha=0.05)
}
# compute worst direction for the worst case
design.names <- names(as.data.frame(design))
design <- design[ , pair.worst]
design.names <- design.names[pair.worst]
anglewise.stat.max <- max(anglewise.stat)
index.max <- which.max(anglewise.stat)
cos.theta.max <- subdiv.halfcircle[index.max, 1]
sin.theta.max <- subdiv.halfcircle[index.max, 2]
dir.max <- c(cos.theta.max, sin.theta.max)
# --------------------
# graphics if required
# --------------------
if (graphics>=0) {
if (is.element(graphics, c(0,1))) {
par(mfrow=c(1,2+graphics))
plot(design, xlim=c(-1,1), ylim=c(-1,1),
xlab=design.names[1], ylab=design.names[2])
}
# draw the rss curve
if (graphics>0) {
rx <- c(rss.curve.x, -rss.curve.x, rss.curve.x[1])
ry <- c(rss.curve.y, -rss.curve.y, rss.curve.y[1])
graph.size <- max(abs((anglewise.stat.max)*dir.max), gof.test.stat)*1.05
plot(rx, ry, xlim=c(-graph.size,graph.size), ylim=c(-graph.size, graph.size),
xlab="", ylab="", ...)
# draw the circle with radius equal to the threshold at significance level 5%
theta_aux <- seq(from=0, to=2*pi+0.1,by=0.1)
lines(gof.test.stat*cos(theta_aux), gof.test.stat*sin(theta_aux))
# draw the coordinate axis in dotted lines
abline(h=0, v=0, lty=lines.lty, col="black", lwd=lines.lwd)
}
if (is.element(graphics, c(0,1))) {
plot(design, xlim=c(-1,1), ylim=c(-1,1),
xlab=design.names[1], ylab=design.names[2])
projections <- design %*% dir.max
points(projections*dir.max[1], projections*dir.max[2], pch=20, col="red")
if (cos.theta.max==0) {
lines(c(0,0), c(-1,1), col="red")
} else lines(c(-1,1), c(-1,1)*sin.theta.max/cos.theta.max, col="red")
for (i in 1:n) lines(c(design[i,1], projections[i]*cos.theta.max), c(design[i,2], projections[i]*sin.theta.max), lty=lines.lty, lwd=lines.lwd)
}
par(mfrow=c(1,1))
}
return(list(global.stat=global.stat, worst.case=pair.worst, worst.dir=dir.max, stat=as.numeric(anglewise.stat), angle=as.numeric(theta), curve=cbind(c(rss.curve.x,-rss.curve.x), c(rss.curve.y,-rss.curve.y)), gof.test.stat=gof.test.stat))
}
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.