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
|
ad.pval<-
function (tx,m,version=1)
{
# This function "adk.pval" evaluates the p-value of the observed value
# tx of the T_m statistic in "K-Sample Anderson-Darling Tests" by F.W. Scholz
# and M.A. Stephens (1987), Journal of the American Statistical Association,
# Vol 82, No. 399, pp 918-924. Thus this p-value is P(T_m >= tx).
#
# Input: tx = observed value of T_m, tx > 0.
# m = the index of T_m, m >= 1.
# version = 1 (default) uses the first version of the AD statistic,
# otherwise the second version is used.
#
# Output: a list with components
# p0 = p-value of tx, i.e., p0 = P(T_m >= tx)
# extrap = a logical indicator
# extrap = TRUE means that linear extrapolation took place
# extrap = FALSE means that quadratic interpolation was used.
#
# Computational Details:
#
# This function uses the upper T_m quantiles as obtained via simulation of
# the Anderson-Darling test statistics (Nsim = 2*10^6) with sample sizes n=500
# for each sample, and after standardization, in order to emulate the Table 1
# values given in the above reference. However, here we estimate p-quantiles
# for p = .00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,
# .1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,
# .99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999
# First the appropriate p-quantiles are determined from those simulated
# for ms = 1,2,3,4,6,8,10, Inf, interpolating to the given value of m
# by fitting a smoothing spline in 1/sqrt(ms) to the simulated quantiles.
# Visual inspection of plots of fitted splines versus 1/sqrt(ms) shows
# good agreement for the used smoothing parameter spar = .4.
#
# Next a smoothing spline is used to fit the log((1-p)/p) to
# these interpolated quantiles and the value fitted to tx is
# obtained (extrapolating linearly) beyond p = .00001 and .99999.
#
# The p-values from Table 1 were reproduced mainly with relative error
# bounded by 1%, in 6 cases with respective relative errors of 1.25%, 1.2%,
# 1.3%, 1.5%, 2.7% and 2.8% as can be seen from the relative error table below.
# The columns correspond to Table 1 p-values of .25, .10, .05, .025, .01 and
# the rows correspond to m = 1, 2 3, 4, 6, 8, 10, Inf.
#
# [,1] [,2] [,3] [,4] [,5]
# [1,] -0.0012 0.0011 0.0038 0.0123 0.0277
# [2,] -0.0105 -0.0109 -0.0022 0.0129 0.0271
# [3,] -0.0077 -0.0052 0.0042 0.0106 0.0152
# [4,] -0.0050 -0.0031 0.0035 0.0067 0.0057
# [5,] -0.0041 -0.0033 0.0055 0.0014 -0.0050
# [6,] -0.0048 -0.0026 0.0002 -0.0021 -0.0082
# [7,] -0.0070 -0.0042 -0.0023 -0.0036 -0.0119
# [8,] 0.0010 0.0005 0.0013 0.0017 0.0023
#
#
#
#
# Fritz Scholz, March 2012
#=========================================================================================
if(version==1){
table1.adk <- structure(c(1, 2, 3, 4, 6, 8, 10, Inf, -1.1954, -1.5806, -1.8172,
-2.0032, -2.2526, -2.4204, -2.5283, -4.2649, -1.1786, -1.5394,
-1.7728, -1.9426, -2.1685, -2.3288, -2.4374, -3.8906, -1.166,
-1.5193, -1.7462, -1.9067, -2.126, -2.2818, -2.3926, -3.719,
-1.1407, -1.4659, -1.671, -1.8105, -2.0048, -2.1356, -2.2348,
-3.2905, -1.1253, -1.4371, -1.6314, -1.7619, -1.9396, -2.0637,
-2.1521, -3.0902, -1.0777, -1.3503, -1.5102, -1.6177, -1.761,
-1.8537, -1.9178, -2.5758, -1.0489, -1.2984, -1.4415, -1.5355,
-1.6625, -1.738, -1.7936, -2.3263, -0.9978, -1.2098, -1.3251,
-1.4007, -1.4977, -1.5555, -1.5941, -1.96, -0.9417, -1.1187,
-1.209, -1.2671, -1.3382, -1.379, -1.405, -1.6449, -0.8981, -1.0491,
-1.1235, -1.1692, -1.2249, -1.2552, -1.2755, -1.4395, -0.8598,
-0.9904, -1.0513, -1.0879, -1.1317, -1.155, -1.1694, -1.2816,
-0.7258, -0.7938, -0.8188, -0.8312, -0.8435, -0.8471, -0.8496,
-0.8416, -0.5966, -0.617, -0.6177, -0.6139, -0.6073, -0.5987,
-0.5941, -0.5244, -0.4572, -0.4383, -0.419, -0.4033, -0.3834,
-0.3676, -0.3587, -0.2533, -0.2966, -0.2428, -0.2078, -0.1844,
-0.1548, -0.1346, -0.1224, 0, -0.1009, -0.0169, 0.0304, 0.0596,
0.0933, 0.1156, 0.1294, 0.2533, 0.1571, 0.2635, 0.3169, 0.348,
0.3823, 0.4038, 0.4166, 0.5244, 0.5357, 0.6496, 0.6992, 0.7246,
0.7528, 0.7683, 0.7771, 0.8416, 1.2255, 1.2989, 1.3202, 1.3254,
1.3305, 1.3286, 1.3257, 1.2816, 1.5262, 1.5677, 1.5709, 1.5663,
1.5561, 1.5449, 1.5356, 1.4395, 1.9633, 1.943, 1.919, 1.8975,
1.8641, 1.8389, 1.8212, 1.6449, 2.7314, 2.5899, 2.5, 2.4451,
2.3664, 2.3155, 2.2823, 1.96, 3.7825, 3.4425, 3.2582, 3.1423,
3.0036, 2.9101, 2.8579, 2.3263, 4.1241, 3.716, 3.4984, 3.3651,
3.2003, 3.0928, 3.0311, 2.4324, 4.6044, 4.0847, 3.8348, 3.6714,
3.4721, 3.3453, 3.2777, 2.5758, 5.409, 4.7223, 4.4022, 4.1791,
3.9357, 3.7809, 3.6963, 2.807, 6.4954, 5.5823, 5.1456, 4.8657,
4.5506, 4.3275, 4.2228, 3.0902, 6.8279, 5.8282, 5.3658, 5.0749,
4.7318, 4.4923, 4.3642, 3.1747, 7.2755, 6.197, 5.6715, 5.3642,
4.9991, 4.7135, 4.5945, 3.2905, 8.1885, 6.8537, 6.2077, 5.8499,
5.4246, 5.1137, 4.9555, 3.4808, 9.3061, 7.6592, 6.85, 6.4806,
5.9919, 5.6122, 5.5136, 3.719, 9.6132, 7.9234, 7.1025, 6.6731,
6.1549, 5.8217, 5.7345, 3.7911, 10.0989, 8.2395, 7.4326, 6.9567,
6.3908, 6.011, 5.9566, 3.8906, 10.8825, 8.8994, 7.8934, 7.4501,
6.9009, 6.4538, 6.2705, 4.0556, 11.8537, 9.5482, 8.5568, 8.0283,
7.4418, 6.9524, 6.6195, 4.2649), .Dim = c(8L, 36L), .Dimnames = list(
NULL, NULL))}else{
table1.adk <- structure(c(1, 2, 3, 4, 6, 8, 10, Inf, -1.1976, -1.5824, -1.8195,
-2.005, -2.2546, -2.422, -2.5292, -4.2649, -1.1806, -1.5416,
-1.7747, -1.9434, -2.1687, -2.3301, -2.438, -3.8906, -1.1681,
-1.5212, -1.7479, -1.9078, -2.1268, -2.2827, -2.3937, -3.719,
-1.1427, -1.4677, -1.6724, -1.8115, -2.0059, -2.1363, -2.2359,
-3.2905, -1.1272, -1.4387, -1.6325, -1.7629, -1.9405, -2.0649,
-2.1527, -3.0902, -1.0794, -1.3518, -1.5112, -1.6187, -1.7617,
-1.8545, -1.9182, -2.5758, -1.0504, -1.2997, -1.4425, -1.5362,
-1.6632, -1.7387, -1.7943, -2.3263, -0.999, -1.2109, -1.3259,
-1.4014, -1.4981, -1.5561, -1.5945, -1.96, -0.9428, -1.1196,
-1.2098, -1.2677, -1.3386, -1.3795, -1.4054, -1.6449, -0.8991,
-1.05, -1.1241, -1.1697, -1.2253, -1.2557, -1.2758, -1.4395,
-0.8607, -0.9911, -1.0518, -1.0883, -1.1321, -1.1555, -1.1698,
-1.2816, -0.7264, -0.7944, -0.8192, -0.8315, -0.8437, -0.8473,
-0.8498, -0.8416, -0.597, -0.6173, -0.6179, -0.6141, -0.6074,
-0.5989, -0.5942, -0.5244, -0.4574, -0.4385, -0.4191, -0.4034,
-0.3835, -0.3677, -0.3588, -0.2533, -0.2966, -0.2427, -0.2078,
-0.1844, -0.1548, -0.1347, -0.1224, 0, -0.1007, -0.0168, 0.0305,
0.0596, 0.0934, 0.1157, 0.1295, 0.2533, 0.1573, 0.2638, 0.3171,
0.3482, 0.3825, 0.404, 0.4168, 0.5244, 0.5363, 0.6501, 0.6996,
0.7249, 0.753, 0.7685, 0.7773, 0.8416, 1.2263, 1.2997, 1.3209,
1.3258, 1.3309, 1.329, 1.326, 1.2816, 1.5274, 1.5686, 1.5716,
1.5667, 1.5565, 1.5453, 1.536, 1.4395, 1.9644, 1.944, 1.92, 1.8983,
1.8647, 1.8396, 1.8216, 1.6449, 2.7334, 2.5915, 2.5013, 2.4457,
2.3671, 2.3162, 2.2827, 1.96, 3.7851, 3.4443, 3.2595, 3.1436,
3.0046, 2.9111, 2.8585, 2.3263, 4.1255, 3.7175, 3.4997, 3.3661,
3.2011, 3.0939, 3.0318, 2.4324, 4.6067, 4.0869, 3.8363, 3.6724,
3.4729, 3.3463, 3.278, 2.5758, 5.4121, 4.7248, 4.4032, 4.1812,
3.9369, 3.7819, 3.6977, 2.807, 6.5, 5.5856, 5.1469, 4.8683, 4.552,
4.3284, 4.2229, 3.0902, 6.8324, 5.8302, 5.3678, 5.0769, 4.7332,
4.4933, 4.3654, 3.1747, 7.278, 6.1999, 5.674, 5.3661, 5.0001,
4.7147, 4.5956, 3.2905, 8.1926, 6.8586, 6.2082, 5.8524, 5.4265,
5.115, 4.9571, 3.4808, 9.3096, 7.6673, 6.8522, 6.4825, 5.9934,
5.6135, 5.5147, 3.719, 9.6207, 7.929, 7.1051, 6.6763, 6.1548,
5.8229, 5.7358, 3.7911, 10.1076, 8.2437, 7.4349, 6.9593, 6.3923,
6.0136, 5.9573, 3.8906, 10.8874, 8.9034, 7.8991, 7.4543, 6.9017,
6.4568, 6.2723, 4.0556, 11.8602, 9.5499, 8.5596, 8.0315, 7.4425,
6.9537, 6.6213, 4.2649), .Dim = c(8L, 36L), .Dimnames = list(
NULL, NULL))
}
extrap <- FALSE
mt <- table1.adk[, 1]
sqm1 <- 1/sqrt(mt)
sqm2 <- sqm1^2
tm <- NULL
p <- 1 - c(.00001,.00005,.0001,.0005,.001,.005,.01,.025,.05,.075,
.1,.2,.3,.4,.5,.6,.7,.8,.9,.925,.95,.975,.99,.9925,.995,.9975,.999,
.99925,.9995,.99975,.9999,.999925,.99995,.999975,.99999)
np <- length(p)
lp <- log(p/(1 - p))
sqm0 <- 1/sqrt(m)
for (i in 1:np) {
out <- smooth.spline(sqm1,table1.adk[,i+1],spar=.4)
y <- predict(out,sqm0)$y
tm <- c(tm, y)
}
out <- smooth.spline(tm,lp,spar=.25)
lp0 <- predict(out,tx)$y
p0 <- exp(lp0)/(1 + exp(lp0))
names(p0) <- NULL
p0
}
|
