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
159
160
161
162
163
164
165
166
167
168
169
|
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))
}
|
