我提到的两篇文献在这里:
https://www.dropbox.com/s/1g9xjzv87t7otdh/1268584.pdf
https://www.dropbox.com/s/fov1m4mm3tlqv4a/2334520.pdf
那本书确实是Simultaneous Statistical Inference。亚马逊上有预览,不过预览里关键的几页看不到,我们图书馆有这本书,我懒得爬出去找了。
今天我再想了一下这个问题,意识到它其实比我原先想象的要简单得多,因为我比较懒,能模拟的情况下先考虑模拟,但模拟其实是最劣等的方案,最理想的方案当然是找到解析解,对这个问题来说解析解不太可能。折衷一下,半计算半解析:
若X服从t分布,累积分布函数是F(x),密度函数是f(x),令Y = |X|,那么密度g(y) = 2f(y),y > 0;分布G(y) = 2F(y) - 1。对r个独立的t分布随机变量X1, ..., Xn,令Z = max(|X1|, ..., |Xn|),那么Z分布H(z) = (2F(z) - 1)^r。这都是些比较初等的概率推导。
既然分布函数能以很简单的形式写出来,代码也就大大简化了,可以直接用t分布的分布函数pt()
:
psmm = function(x, r, nu) {<br />
ifelse(x > 0, (2 * pt(x, nu) - 1)^r, 0)<br />
}<br />
# 分位数函数仍然用H(x) - q = 0求解x<br />
qsmm = function(q, r, nu) {<br />
res = uniroot(function(x, r, nu, q) {<br />
psmm(x, r = r, nu = nu) - q<br />
}, c(0, 100), r = r, nu = nu, q = q)<br />
res$root<br />
}</p>
<p># 算0.99分位数<br />
qsmm(.99, r=3, nu=5)</p>
<p>ks = 3:20<br />
vs = c(5, 7, 10, 12, 16, 20, 24, 30, 40, 60, 120, Inf)<br />
q0.01 = matrix(nrow = length(ks), ncol = length(vs), dimnames = list(ks, vs))<br />
for (k in ks) {<br />
for (v in vs) {<br />
q0.01[as.character(k), as.character(v)] = qsmm(.99, r = k*(k-1)/2, nu = v)<br />
}<br />
}<br />
cbind('r*' = choose(ks, 2), round(q0.01, 3))</p>
<p> r* 5 7 10 12 16 20 24 30 40 60 120 Inf<br />
3 3 5.243 4.353 3.825 3.647 3.443 3.329 3.257 3.188 3.121 3.056 2.994 2.934<br />
4 6 6.133 4.941 4.256 4.029 3.771 3.629 3.539 3.453 3.370 3.291 3.215 3.143<br />
5 10 6.862 5.404 4.584 4.315 4.013 3.848 3.744 3.644 3.549 3.459 3.372 3.289<br />
6 15 7.491 5.791 4.852 4.547 4.206 4.021 3.905 3.794 3.689 3.589 3.493 3.402<br />
7 21 8.051 6.126 5.079 4.742 4.368 4.165 4.038 3.918 3.803 3.695 3.591 3.493<br />
8 28 8.558 6.424 5.278 4.912 4.506 4.288 4.152 4.023 3.900 3.784 3.673 3.569<br />
9 36 9.024 6.693 5.454 5.061 4.628 4.396 4.251 4.114 3.984 3.861 3.744 3.634<br />
10 45 9.457 6.939 5.614 5.196 4.737 4.491 4.339 4.194 4.058 3.929 3.807 3.691<br />
11 55 9.862 7.166 5.760 5.319 4.835 4.577 4.417 4.267 4.124 3.989 3.862 3.742<br />
12 66 10.244 7.378 5.894 5.431 4.925 4.656 4.489 4.332 4.184 4.044 3.912 3.787<br />
13 78 10.606 7.576 6.019 5.535 5.008 4.728 4.555 4.392 4.238 4.094 3.957 3.829<br />
14 91 10.950 7.762 6.135 5.632 5.084 4.794 4.615 4.447 4.288 4.139 3.999 3.866<br />
15 105 11.279 7.939 6.245 5.722 5.156 4.856 4.672 4.498 4.335 4.181 4.037 3.901<br />
16 120 11.595 8.107 6.348 5.807 5.223 4.914 4.724 4.546 4.378 4.221 4.073 3.933<br />
17 136 11.898 8.267 6.446 5.888 5.286 4.968 4.773 4.590 4.418 4.257 4.106 3.963<br />
18 153 12.190 8.419 6.538 5.964 5.345 5.020 4.820 4.632 4.456 4.291 4.137 3.991<br />
19 171 12.472 8.566 6.627 6.037 5.402 5.068 4.863 4.672 4.492 4.324 4.166 4.018<br />
20 190 12.745 8.707 6.711 6.106 5.455 5.114 4.905 4.709 4.526 4.354 4.193 4.042
</p>
文献中的分位数表是下面这个,可以看出自由度无穷大的时候能对得上(正态分布),但自由度比较小的时候他的分位数值总体偏小(相应置信区间的长度也会小):

要么是 a) 参考文献13(Stoline等人)的Fortran程序写的有问题,要么是 b) R的t分布函数写的有问题,要么是 c) 那1966年的书推导有问题。现在看来 a) 和 c) 应该是一致的,因为他们的结果我用暴力积分的办法可以重复出来。那么要么R错了,要么公式错了。希望不是R搞错了,不然感觉真的不会再爱了。
20分钟后更新
===
再看了一遍公式2.2,发现它真的错了耶。

那个幂r应该提到积分的外面去。我了个去,我这榆木脑袋居然能发现一处错误的公式,生平第一次哎。叩谢楼主的这份圣诞礼物。