JAPLAシンポジュウム資料(2004.12.11)

 

       ノンパラメトリック統計学とJ言語

 

                       帝京平成大学 鈴木 義一郎

 

   §1 適合度検定

 一般に、連続な分布関数F(x)に従う確率変数に対して、大きさの順番に並べ直した{ }を「順序統計量」という。さらに

で定義されるものを「経験分布関数」といい、

     Fn(x) → F(x) n → ∞ for all x

という性質がある。そこで

といった統計量が定義できる。これらはコルモゴルフの「D−統計量」と呼ばれているもので、特には「片側統計量」である。またこれらの統計量の分布は、F(x)という分布形に依存しないで定まることから「分布に依らない(distribution-free)統計量」とも呼ばれている。

 一般に、,i=1,2,……… ,nのように変換すると、XがF(x)に従うときには,これらは[0,1]上の一様分布からの順序統計量になる。そこで、当該データから、の経験分布関数([0,1]上の関数になる!)を描いてみて、これが対角線上から極端に離れる(つまりの値が大きくなる)ようならば、一様分布からの順序統計量ではない、したがって「観測データがF(x)には従わない」と判断できる。特にFとして正規分布関数を考えれば、正規分布か否かを識別することができるというになる。

 ところが問題は、この種の統計量が適合度を測っているのでなく、“非適合度”を測っているに過ぎないという点である。つまりデータ数が少ないと、どんな分布を仮定しても、情報不足の故に、その分布に従わないという結論が得られない、つまりその分布に従うことを容認せざるを得ないことになる。結局、どんな分布を仮定しても否定されないということは、「どんな分布を仮定すれば良いか?」という問の回答が得られないということになる。逆に、データ数が非常に多い場合には、どんな分布を仮定しても「その分布には従わない」という結論になるケースが多くなる。結局、どんな分布を仮定しても否定されるということで、やはり「どの分布を仮定すれば良いか?」に対する回答は得られない。要するに、データが少なすぎても多すぎても分布を特定できないということから、「適合度検定は使いものにならない」ケースが圧倒的に多いというのが結論である。

 それはともかくも、片側統計量の確率分布の求め方の概略を紹介してみよう。まず

という値を算出する。そこで

という多項式群を定義して

という関数を求めると

が求める確率分布を与えることになる。

 

 次のデータは、昭和35年の栃・若時代の幕内力士43人の体重を、大小の順に並べたものである。

    85  88  95  97 101 101 102 103 103 104 105 105 105 107 109

   109 110 110 111 111 111 112 112 114 115 115 116 116 116 120

   120 120 122 124 128 128 131 135 141 146 148 150 150

このデータを「WAT」という変数に入力して、これを「stand」という関数によって基準化すると

   0.3": 4 10 $ stand WAT

_1.604 _1.478 _1.257 _1.132 _0.943 _0.659 _0.659 _0.565 _0.502 _0.502

_0.471 _0.471 _0.439 _0.408 _0.408 _0.376 _0.345 _0.282 _0.219 _0.502

_0.156 _0.156 _0.156 _0.124 _0.093 _0.093  0.093  0.285  0.348  0.348

 0.379  0.379  0.474  0.474  0.505  1.229  1.481  1.607  1.922  4.125

0.153 1.179 _0.728 0.298 _0.738 _0.718 _0.515 2.399 _0.709 _0.621

この値に標準正規分布関数Φ(・)をかぶせてみると

   0.3": 4 10 $ E=:ndf"0 /:〜 stand WAT

0.054 0.070 0.104 0.129 0.173 0.255 0.255 0.286 0.308 0.308

0.319 0.319 0.330 0.342 0.342 0.353 0.365 0.389 0.413 0.438

0.438 0.438 0.438 0.451 0.463 0.463 0.538 0.612 0.636 0.636

0.648 0.648 0.682 0.682 0.693 0.890 0.931 0.946 0.973 1.000

といった結果が得られる。


一方、一様分布の理論値のほうは

     i/n,(t=i/n) ,i=1,2,…………,n

であるが、両側統計量の場合の近似に用いるためには

     (i-0.5)/n,(t=i/n) ,i=1,2,…………,n

を考えたほうがよい。

   4 10 $ T=:((i.+0.5"_)%])40

0.0125 0.0375 0.0625 0.0875 0.1125 0.1375 0.1625 0.1875 0.2125 0.2375

0.2625 0.2875 0.3125 0.3375 0.3625 0.3875 0.4125 0.4375 0.4625 0.4875

0.5125 0.5375 0.5625 0.5875 0.6125 0.6375 0.6625 0.6875 0.7125 0.7375

0.7625 0.7875 0.8125 0.8375 0.8625 0.8875 0.9125 0.9375 0.9625 0.9875

 そこで、{E−T}の絶対値をとると

   4 10 $ D=:|E-T

0.042 0.032 0.042 0.041 0.060 0.117 0.092 0.099 0.095 0.070

0.056 0.031 0.018 0.004 0.021 0.034 0.047 0.048 0.049 0.049

0.074 0.099 0.124 0.137 0.149 0.174 0.124 0.075 0.076 0.101

0.115 0.140 0.130 0.155 0.169 0.003 0.018 0.008 0.010 0.012

のような結果になる。これらの値の最大値は

   >./D

0.174489

である。また、「stat_d」という関数を適用すれば

   stat_d stand WAT

0.174489

のように同じ結果が得られる。

ところで片側統計量の5%点や2.5%点は、「percent_5」や「percent_2_5」という関数を用いて

   percent_5 40

0.210115

   percent_2_5 40

0.18913

のように求められる。片側統計量の 2.5%点がほぼ両側統計量の5%点になるので、Dnの値が0.21より大きければ正規分布の仮定が棄却できる。観測結果は 0.1745であるから、有意水準5%では棄却できない。つまり正規分布でないとは結論できない。また最近の幕内力士の体重データも「WTW」という変数に入力して、

   stat_d stand WTW

0.140572

といった結果で、「WAT」の場合より偏差が小さく、より正規分布に近いことが分かる。

 さらに、一様分布からのデータについても試してみると

   stat_d stand ([:?]$10000"_) 40

0.0959442

   stat_d stand ([:?]$10000"_) 40

0.103789

といった例が示すように、正規分布であるとの仮説棄却されない。つまりデータ数が40個程度では、一様分布や2重指数分布などからの乱数を使って実験してみても、正規分布であるとの仮説はほとんど否定できない!

 

【Jによるプログラム】

mvalue=:0:>.(-@{:*%@{.)+(>:@i.%])@{. NB. m(i)=m(i;n,d)(1≦i≦n)

weight=:(i.!])@[*(<:@[{mvalue@])^[:|.1:+i.@[

        NB. W(k,i)=B(k,i)m(k)^(k-i),(0≦i≦k-1)

qvalue=:>@{.;>@{:,[:-[:+/#@>@{:weight>@{.

    NB. Q(k)=Q(n,k)=-{Q(0)W(k,0)+Q(1)W(k,1)+  +Q(k-1)W(k,k-1)

function=:4 :'+/((1.x.+1)!x.)*,>}.qvalue^:x.(x.,y.);1'

        NB. P(n,c)=B(n,0)Q(n,0)+B(n,1)Q(n,1)+  +B(n,n)Q(n,n)

derivative=:4 :'(x."_function])D.1 y.'

ratio_5=:(function-0.95"_)%derivative

percent_5=:([:([:{:({.,{:-{.ratio_5{:)^:_)],1:)%]

ratio_2_5=:(function-0.975"_)%derivative

percent_2_5=:([:([:{:({.,{:-{.ratio_2_5{:)^:_)],1:)%]

nden=:[:%&(%:o.2)[:^*:%_2: NB. standard normal density]

grid=:*&((>i.499)%500)

ndf=:0.5"_+%&500*([:mean[:nden 0:,])+[:+/[:nden grid

       NB. standard normal distribution function

stat_d=:[:>./[:|([:ndf"0/:〜)-[:((i.+0.5"_)%])# NB. D-stat.

 


   §2 順位和検定(Mann-Whitney test)

 §1では、データ数が少なすぎても多すぎても分布が特定できないことを指摘したが、データが多ければ、ある種の条件を満たす限り『中心極限定理』が成り立つことから、正規分布を仮定して議論を行なうことができる。問題は、データ数が少なすぎる場合である。正規分布であると確信できれば、t−推定やt−検定を用いることができるが、どんな分布に従うのか判然としない場合、強引に正規性を仮定して議論するのは危険である。そこで、分布形を規定しなくても利用できる、いわゆる「ノンパラメトリック法」を用いて検定などを行なえばよいということになる。

 例えば

   X=:9.3 19.9 0.2 10.8 0.1 0.3 2.4 32.5 0.4 1.3

   Y=:8.4 23.5 113 97.3 53.2 28.7 7.6 44.7 85.2 49.1 105.1 30.3

といった2組のデータが与えられていて、平均の差異について検定したいとする。

 まず平均は

   (MX=:mean X),(MY=:mean=:+/%#)Y

7.72 53.8417

であるから、明らかに差がありそうである。

 さらに分散も

   (VX0=:var X),(VY0=:var=:[:mean[:*:]-mean)Y

106.716 1284.64

と算出されるから、等分散も仮定できない。

 そこでまず、正規分布を仮定できるものとして議論してみよう。まず最初に、不偏分散を求めると

   (VX=:var_u X),VY=:(var_u=:[:(+/%<:@#)[:*:]-mean)Y

118.573 1401.43

である。さらに

     ww=[VX/m+VY/n]2/{VX2/m2(m-1)+VY2/n2(n-1)}

という値を超える最小の整数をwとするとき、

という統計量が、等平均の仮定の下では自由度wのt−分布で近似できることを利用して検定できる。これは「ウェルチの検定」と呼ばれているものである。wwや統計量tの値を求める関数は次のように与えられる。

   welch0=:[:(*:@{.%{:)[:+/"1[:>(](],*:@]%<:@#@[)var_u%#)L:0

   welch=:[:(-/@{.%[:%:+/@{:)[:|:[:>(mean,var_u%#)L:0


この例では、

   (welch0 X;Y),welch X;Y

0.556572 _4.06642

である。tの絶対値4.066は、自由度(w=1)のt−分布の片側5%点6.314より小さいから、等平均の仮説は棄却できない。平均がかなり異なっている感じなのに等平均の仮説を棄却できないということは、正規分布の仮定がオカシイということになる。

 そこで、分布形によらない「ノンパラメトリック検定」を試してみよう。平均の差異を検出する最も典型的なものが、「順位和検定」とか「Mann-Whitney検定」と呼ばれている方法である。

 一般に、同じ連続分布から

    

のような2組のデータが与えられているとする。全体を一緒にして並べ直したときの、Xのほうの順位を()として

という「順位和」を考える。そこで、次のような関数を考える。

定義から、次のような漸化関係の成り立つことが分かる。

 したがって、mやnの値がそれほど大きくない場合には、正確な%点の値を直接求めることができる。

 さらに、Wの期待値と分散は

     E{W} = m(m+n+1)/2 = M ,V{W} = mn(m+n+1)/12 = nM/6

で与えられる。また4次のキュミュラントも

     K = -mn(m+n+1)×{(m+n)(m+n+1)-mn}/120 = -nM{M+nM/m-mn/4}/30

のように与えられるから、

と置いて

     F(t;m,n) = Φ(V)

という関係によって近似することができる(Φは標準正規分布関数)。

 さて、「rank」という関数は2組のデータ全体に順位を与えるものである。

   rank X,Y

9 11 2 10 1 3 6 15 4 5 8 12 22 20 18 13 7 16 19 17 21 14

   +/ 10{. rank X,Y

66

   X ranksum Y

66

 つまり「ranksum」は、Xの順位の合計Wである。なお、WにYのほうの順位和183を合計すると253になり、これは(m+n)(m+n+1)/2のように、標本数によって定まる量である。つまり、XかYのどちらか一方の順位和だけに着目すればよい。さらに、m=10,n=12のときの5%点(の近似値)を求める関数「rminimax」を用いて

   rminimax 10,12

89 141

という値が算出できる。この近似式は驚くほど精度がよい(m≦n≦12という範囲では、(3,9)と(5,10)のときだけ1だけ安全なほうに異なっている)。

 さて与えられた順位和の66は、下方の限界値89より小さく、メデタク、平均に差のあることが検出できた。また、力士の体重のデータにも適用してみると

   WTW ranksum WAT

1133

が順位和である。この場合の5%点が

   rminimax 40,43

1498 1862

であるから、平均差は高度に顕著であることが分かる。さらに

   W1=:88 103 110 110 112 124 128 141 148

   W2=:113 132 137 138 139 142 144 148 165 204

は、それぞれ WTW,WATというデータから10個ずつランダムに抜き出したものである。このような小さめのデータにも適用してみると

   W1 ranksum W2

70

が順位和である。この場合の5%点は

   rminimax 10,10

82 128

であるから、やはり平均に差のあることが検出できる。

 


【Jによるプログラム】

NB. TEST STATISTIC

rank=:1:+(i.@#=//:)+/ .*i.@#

ranksum=:[:+/#@[{.[:rank, NB. rank-sum test statistic(dyadic)

NB. Exact Probability

add=:(([:<"1[:(],.[:}.i.)1:+#),.]);(<<1 0)"_

left=:>@{.@{.>@{.

right=:{:@{.@{.,{:@>@{:

middle=:([:+/[*(({:@[$0:),>@{.@]),>@{:@])%+/@[

next=:}.@>@{.;[:<>@{:,[:<left middle right

first=:3 :'}.>{:next^:(#y.)add y.'

second=:([$0:)(,+,〜)-:@]

prob_list=:[:(],[:<>:@#@]second[:([:>{:)first

start=:[:<:[:([:+/>:@i."_1>:@i.

minpt=.start@#+([:+/0.05"_>:+/\)&>

total=:(1:+i."0)([*>:@+)]

minmax=:[:(],.total@#-])minpt

percent=:}."1,([:(minmax;<@,.@])prob_list>@{:"1

NB. Approximation

rm=:[:-:<./**1:++/

ru=:%:@(12"_%*/*1:++/)@]*[-0.5"_+rm@]

rv=:ru*1:+(*:@ru-3:)*(((]*>:)@+/-*/)%40"_*}:*rm)@]

range=:-:@(]*>:)@{.+[:i.*/ NB. a=m(m+1)/n,a+1,---,b=a+mn

rmin=:_1:+[:{:(_1.65"_>:[:,range rv"0 1])#range

rminimax=:(*/+(]*>:)@<./)(],.-)rmin NB. lower & upper 5%-points

 

   §3 置換検定(Pittman test)

 §2では、平均の差異を検出する最も典型的なノンパラメトリック検定として、「順位和検定」を考えた。もし与えられたデータが「順序尺度」であれば、順位だけを利用する方法を考えるのが自然である。しかしデータが「間隔尺度」の場合には、順位だけでなく数値まで利用する考えのほうが合理的である。そこで、ピットマンが提案した「置換検定」という方法を紹介してみよう。

 一般に、同じ平均をもつ母集団から

    

のような2組のデータが与えられているとする。全体を一緒にして並べ直したものを

    

とする。(m+n)個からm個を取り出す組合せの数はN=あり、各組合せのm個の和を小さいほうから順番に並べたものを

(*) 

とする。最小値と最大値は、明らかに

である。そこで、Xのほうのデータの和が(*)のどの辺に位置するかに注目する。もし平均に差がないとすれば、真中辺に位置するし、Xの平均がYの平均より小さい(大きい)場合には、左端(右端)のほうに位置することが期待される。結局、Xのほうの和が左右いずれかの極端に偏っているときに、等平均の仮説が棄却できることになる。ただこの方法は、mやnがほどほどに小さくないと、組合せの数が大きくなって計算量が膨大になるという欠陥もある。例えば、m=10n=12としてみると、N=646646である。

 まず

   A=:3 4 ] B=:2 6 7

といった2組のデータで考えてみる。m=2,n=3の場合の組合せは

   comball 2 3

1 1 0 0 0

1 0 1 0 0

1 0 0 1 0

1 0 0 0 1

0 1 1 0 0

0 1 0 1 0

0 1 0 0 1

0 0 1 1 0

0 0 1 0 1

0 0 0 1 1

のように10通りある。各組合せに対する2個のほう和は

   2 5([:+/"1([:comball#&>@;)#,)3 6 7

7 5 8 9 8 11 12 13 10 9

のように与えられる。

 そこで、これらの値を順番に並べ直したとき、Aという組合せの和がどの位置にあるかは、「position」という関数を用いて

   2 5 position 3 6 7

1 1 0 0 0 0 0 0 0 0

のように小さいほうから2番目にあることが分かる。


   2 5 permute 3 6 7

0.2

 次に、m=6,n=8の場合の例として

   ]XX=:4}./:〜X

1.3 2.4 9.3 10.8 19.9 32.5

   ]YY=:_4}./:〜Y

7.6 8.4 23.5 28.7 30.3 44.7 49.1 53.2

といった2組のデータについて考えてみる。順位和検定を行なってみると

   XX ranksum YY

32

という結果が得られる。しかし、順位和の5%点は

   rminimax 6 8

31 59

であるから、等平均の仮説が棄却されない。ところが置換検定では

   +/ XX position YY

78

となり、XXの和は小さいほうから78番目のところにあることが分かる。N=3003であるから

   78 % 3003

0.025974

また、「permute」という関数を使って

   XX permute YY

0.025974

となる。5%の有意水準で、等平均性が棄却できる。

 

【Jによるプログラム】

zero=:[:<(1:,1:+[:{:[:$[:>{.)$0:

body=:[:(((1:,.>@{.),0:,.>@{:)L:1)2:<\]

more=:]`(],[:-.L:0{:)@.(([:([:(+/<-:@#)@{.@>{:)])

comball=:3 :'>({.y.){(zero,[:more body)^:(_1++/y.)(,.0);,.1'

position=:+/@[>:[:+/"1([:comball#&>@;)#,

permute=:([:+/position)%[:({.!+/)#&>@;

 


   §4 連の個数による検定

 §2のデータを大きさの順に並べて

   X=:0.1 0.2 0.3 0.4 1.3 2.4 9.3 10.8 19.9 32.5

   Y=:7.6 8.4 23.5 28.7 30.3 44.7 49.1 53.2 85.2 97.3 105.1 113

データ全体に順位をつけたとき、Xの方の順位は

     1 2 3 4 5 6 9 10 11 15

のように与えられ、またYのほうの順位は

     7 8 12 13 14 16 17 18 19 20 21 22

である。これらの順番にしたがって、Xのデータには1、Yのデータには0をつけると

     1 1 1 1 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0

といった系列が得られる。このとき、1の連なりや0の連なりのことを「連」といい、各連の大きさを連の「長さ」という。上の例では、長さ6の1の連、長さ2の0の連、長さ3の1の連、長さ3の0の連、長さ1の1の連、長さ7の0の連があり、合計6つの連がある。連の個数をR=R(m,n)と書くと、1と0が左右に完全に分離した場合は連の個数Rが2で、これが最小値である。一般に、mやnの数が大きくなればRの値も大きくなる。また、2組の標本が同一母集団からのものであれば1や0が入り乱れて連の個数は多くなる傾向を示すが、平均などが違っていれば、1と0が分離しやすくなるから連の個数Rは少なくなる。かくて、Rの大小によって、2組の標本が同一母集団からのものであるか否かを識別することができる。

 一般に、Rの確率分布は

で与えられる。ここで

 この確率の計算はWの場合より簡単で、かなり大きなm,nでも正確な値を瞬時に算出できる。この統計量は、m,nがかなり大きいときにだけ正規近似を用いるほうがよい。

 Rの期待値と分散は

のように与えられる。

 さて、「run」という関数を用いれば

   X run Y

1 1 1 1 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0

のように0,1の系列が得られ、さらに「number」という関数により


   X number Y

6

といったように、連の個数が求められる。

 次に、「prob」という関数を用いれば、Rの正確な確率分布が算出できる。

m=10,n=12の場合には

   0.4": 3 7 $ prob 10 12

0.0000 0.0000 0.0003 0.0014 0.0061 0.0163 0.0429

0.0750 0.1286 0.1543 0.1800 0.1500 0.1200 0.0686

0.0367 0.0138 0.0046 0.0010 0.0002 0.0000 0.0000

である。さらに

   0.05(lower,upper)10 12

7 17

が、下側と上側の5%点を与えている。上の例ではR=6であったから、下側5%点より小さく、等平均の仮説を有意水準5%で棄却できる。この結論は、順位和検定の場合と同じである。

 また力士の体重のデータにも適用してみると

   WTW number WAT

22

が連の個数である。またこの場合の5%点は

   0.05(lower,upper)40 43

34 51

と与えられ、連の数は極端に少ないことから、やはり平均体重の間にはかなりの差のあることが識別できる。さらに小さなW1,W2といったデータにも適用してみると

   W1 number W2

8

が連の個数である。またこの場合の5%点は

   0.05(lower,upper)10 10

6 16

であるから、平均が等しいとの仮説は棄却できない。順位和検定では棄却できたから、連の数による検定は、やや感度が鈍い趣のあることが分かる。

 


【Jによるプログラム】

run=:(1:+i.@#@,)e.#@[{.[:rank,

number=:[:(1:+[:+/[: }.-}:)run NB. number of run test

even=:2:*[:*/! NB. 2*C(m,k)*C(n;k) x=k,y=m,n

odd=:([:*/@["0])+[:*/[!"0] NB. C(m,k)*C(n;k-1)+C(m,k-1)*C(n;k)

comb=:(<:@-:@[even<:@])`(([:(<:,])<.@-:)@[odd<:@])@.(2: [)

prob=:((1:+[:}.i.@+/)comb"0 1])%<./!+/

lower=:1:+[:[:+/[>:[:+/\prob@]

upper=:2:+[:[:+/[<:[:+/\[: .prob@]

 

   X run Y

1 1 1 1 1 1 0 0 1 1 1 0 0 0 1 0 0 0 0 0 0 0

   X number Y

6

   0 even 3 2                     1 2 3 even"0 1(3 2)

2                              ┃12 6 0

   ([:(,.|.)],<:)2                2

2 1

1 2

   2(([:(,.|.)],<:)@[);])3 2      2(([:(,.|.)],<:)@[)!"1])3 2

-─┬-─┐                   ┃3 1

2 1│3 2│                   ┃3 2

1 2│  

└─-┴-─┘

   2 odd 3 2                      1 2 3 odd"0 1(3 2)

9                              ┃5 9 1

   (2+i.6)comb"0 1(4 3)

2 5 12 9 6 1

   prob 4 3

0.05714429 0.142857 0.342857 0.357143 0.171429 0.0285714

   +/ prob 4 3                    0.05 lower 4 3

1                              ┃1

   ($,+/) prob 10 12              (lower,upper)prob 10 12

21 1                           ┃7 17


   §5 Mosesの検定

 一般に

のような2組のデータが与えられているとき、全体を大きさの順番に並べたものを

と置いて、X(Y)のデータには1(0)を付ける。そこで、hを適当な自然数として、左から(h+1)番目の1の順位を、右から(h+1)番目の1の順位をとして

という統計量を考え

となるような限界値を求めて、観測結果がこの値より小さいときに、散布度が等しいとする帰無仮説を有意水準αで棄却できる。hの定め方には任意性があるが、h=[m/4]を考えるのが妥当なところだろう。そこで、特にh=[m/4]としたをM(m,n)と表すことにする。

 Moses([2])は、k=m−2hと置いて、次の関係を示した。

【J言語によるプログラム】

   position=:([:/:〜;)e.>@{.

   stat_m=:3 :'-/+/(+/\position y.)<:/(#>{.y)(-,])<.y.%4'

   pair=:[:(,: .)[:i.>:

   start=:](<:@],-)[:<:[:<.[:-:[

   numerator=:[:*/"[:!/"1[: :start(],:+)pair@[

   prob_m=:[:(]%+/){.numerator{:

   prob5=:[:((1:+[:+/0.05"_>:]){.])[:+/\prob_m

   index=:[:+/"1[>:/[:prob5]

   range0=:0:,([:i.>:)+2:*[:<.]%4:

   percent_m=:index{[:range0[:{.]

 

     X=.6 8 9 10 11 12 13 14 15 16 20 24

     Y=.2 3 4 5 7 21 22 23 25 26 30 32

   position X;Y

0 0 0 0 1 0 1 1 1 1 1 1 1 1 1 1 0 0 0 1 0 0 0 0

   3 stat_m X;Y

6

   pair 8                         8([: :start@](],:+)pair@[)10

0 1 2 3 4 5 6 7 8              ┃ 1 6435

8 7 6 5 4 3 2 1 0              ┃ 3 3432

   8 start 10                  ┃ 6 1716

2 7                            ┃10  792

   A=.8(start@](],:+)pair@[)10 ┃15  330

 0  1  2  3  4 5 6  7  8       ┃21  120

 8  7  6  5  4 3 2  1  0       ┃28   36

                               ┃36    8

 2  3  4  5  6  7 8 9 10       ┃45    1

15 14 13 12 11 10 9 8  7

 

   8 numerator 10

6435 10296 10296 7920 4950 2520 1008 288 45

 

   0.3": prob_m 8 10

0.147 0.235 0.235 0.181 0.113 0.058 0.023 0.007 0.001

   prob5 8 10                     prob5 12 12

0.147059                       ┃0.0186335 0.0774763

   prob5 16 16

0.00339884 0.0186256 0.0566926

0.05 を1つ超えたところまでの累積確率の系列を出力する】

 

   (p=.0.01 0.025 0.05)index 12┃   p index 16

0 1 1                          ┃1 2 2

   range0 12

0 6 7 8 9 10 11 12 13 14 15 16 17 18

   range0 16

0 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24

   p percent_m 12 12              p percent_m 16 16

0 6 6                          ┃8 9 9

   p percent_m 20 20              p percent_m 24 24

11 12 13                       ┃14 15 16

 


   p percent_m"1(12,.12+i.5)      p percent_m"1(16,.16+i.5)

0 6 6                          ┃8 9 9

0 6 6                          ┃8 9 9

0 0 6                          ┃8 9 9

0 0 6                          ┃9 9 9

0 0 6                          ┃8 9 9

 

   §6 アンサリ・ブラッドレィの検定

 一般に

   

のような2組のデータが与えられているとする。全体を大きさの順番に並べてたものを

    

と置いて、

という1,0の系列に

というウエィトを掛けて

という統計量を定義する。この統計量の範囲は

     a = {c(c+1)+d(d+1)}/2 , c=[m/2],d=[(m+1)/2]

     b = a+d×[n/2]+c×[(n+1)/2]

である。そこで、Aという統計量がxという値になる、m個の1とn個の0の組合せの数をとすると、Aの確率関数は

のように与えられる。

 Ansari-Bradley([1])は

といった漸化関係の成り立つことを示した。さらに期待値と分散も

のように与えられる。

 次のような2組のデータが与えられたとする。

   A=:47 58 33 40 65 28 39 54 46 60 51 69

   B=:53 44 49 56 43 55 45 57 59 52 49 48

そこで、2組のデータを大きさの順番に並べて、X(Y)のデータには1(0)を付けると

     U = 1 1 1 1 0 0 0 1 1 0 0 0 1 0 0 1 0 0 0 1 0 1 1 1

といった系列が得られ、これに

     R = 1 2 3 4 5 6 7 8 9 10 11 12 12 11 10 9 8 7 6 5 4 3 2 1

というウエィトを掛けて加えると、統計量の値としては

     A = 59

のように算出される。ところで

     Pr{A(12,12)≦59} = 0.0151

という関係が得られるから、有意水準3%で等分散の帰無仮説を棄却することができる。

 

【J言語によるプログラム】

   position=:([:/:〜;)e.>@{.

   weight=:[:;[:(1:+[:i.<.)L:0[:<.[:(;-)-:

   stat_ab=:[:+/[:*/((weight#@;),:position

   first=:((([:>.]%2:)}.(2: >:)$1:)"0

   length=:1:+[:+/ .*[:[:(<.,>.)&>;

   min=:[:+/]{.2:#1:+i.

   range=:([:min[)+[:i.length

   left=:[:+/"1[>:/([:(]%{:)+/\)@]"1

   right=:[:+/"1[>:/([:(]%{.)+/\.)@]"1

   both=:(left, .@[right])`(left,:right)@.(1:=#@[)

   percent_ab=:[: :([both[:{.]){[:}.]

 

   more=:4 :0

r=. :,.f=. .first x.+h=.0

while. h<(#{."1 y.)-1

  do. b=.(a=.x.length h+1)-#s=.((>0:)#](h=.h+1){y.

      r=.r,f=.((b$0),s)+f,(a-#f)$0

end.

)


   table=:4 :0

r=.<q=.((k=.1)range y.),first 1+i.y.

while. k<x.

  do. r=.r,<q=.(k range y.),(k=.k+1)more}.q

end.

)

   position A;B

1 1 1 1 0 0 0 1 1 0 0 0 1 0 0 1 0 0 0 1 0 1 1 1

   weight 24

1 2 3 4 5 6 7 8 9 10 11 12 12 11 10 9 8 7 6 5 4 3 2 1

   */ ((weight#@;),:position)A;B

1 2 3 4 0 0 0 8 9 0 0 0 12 0 0 9 0 0 0 5 0 3 2 1

   +/ */ ((weight#@;),:position)A;B

59

   stat_ab A;B

59

  first 1+i.5                     min"0(1+i.5)

2 0 0                          ┃1 2 4 6 9

2 1 0                             2 length"0(1+i.5)

2 2 0                          ┃2 3 4 5 6

2 2 1                             3 length"0(1+i.5)

2 2 2                          ┃2 4 5 7 8

 

   2 range"0(1+i.5)               3 range"0(1+i.5)

2 3 0 0 0 0                    ┃4 5 0 0 0 0  0  0

2 3 4 0 0 0                    ┃4 5 6 7 0 0  0  0

2 3 4 5 0 0                    ┃4 5 6 7 8 0  0  0

2 3 4 5 6 0                    ┃4 5 6 7 8 9 10  0

2 3 4 5 6 7                    ┃4 5 6 7 8 9 10 11

   2 more first 1+i.5

2 0 0 0 0 0                    ┃f(x 2,1) , x=2,3

1 4 1 0 0 0                    ┃f(x 2,2) , x=2,3,4

1 4 3 2 0 0                    ┃f(x 2,3) , x=2,3,4,5

1 4 5 4 1 0                    ┃f(x 2,4) , x=2,3,4,5,6,7

1 4 5 6 3 2                    ┃f(x 2,5) , x=2,3,4,5,6,7,8


   3 mre 2 more first 1+i.5

2 2  0  0  0  0 0 0            ┃f(x 3,1) , x=4,5

2 3  4  1  0  0 0 0            ┃f(x 3,2) , x=4,5,6,7

2 4  8  4  2  0 0 0            ┃f(x 3,3) , x=4,5,6,7,8

2 4  9  8  7  4 1 0            ┃f(x 3,4) , x=4,5,6,7,8,9,10

2 4 10 12 12 10 4 2            ┃f(x 3,5) , x=4,5,6,7,8,9,10,11

 

   3 table 5

1 2 3

2 0 0

2 1 0

2 2 0

2 2 1

2 2 2

2 3 4 5 6 7

2 0 0 0 0 0

1 4 1 0 0 0

1 4 3 2 0 0

1 4 5 4 1 0 1 4 5 6 3 2

4 5  6  7  8  9 10 11

2 2  0  0  0  0  0  0

2 3  4  1  0  0  0  0

2 4  8  4  2  0  0  0

2 4  9  8  7  4  1  0

2 4 10 12 12 10  4  2

    4 range 8

6 7 8  9 10 11 12 13 14 15 16 17 18 19 20 21 22

   ]T48=.{:4 table 8

6 7 8  9 10 11 12 13 14 15 16 17 18 19 20 21 22

1 2 2  0  0  0  0  0  0  0  0  0  0  0  0  0  0

1 4 5  4  1  0  0  0  0  0  0  0  0  0  0  0  0

1 4 7  8  9  4  2  0  0  0  0  0  0  0  0  0  0

1 4 9 12 18 12  9  4  1  0  0  0  0  0  0  0  0

1 4 9 14 22 22 21 16  9  4  2  0  0  0  0  0  0

1 4 9 16 26 32 34 32 26 16  9  4  1  0  0  0  0

1 4 9 16 28 36 44 46 46 36 29 18 11  4  2  0  0

1 4 9 16 30 40 54 60 67 60 54 40 30 16  9  4  1

   +/"1}.T48

5 15 35 70 126 210 330 495

   4 ! 5+i.8

5 15 35 70 126 210 330 495

   ]A48=.{:T48

1 4 9 16 30 40 54 60 67 60 54 40 30 16 9 4 1

   0.3":100*9{.A48%495

0.202 0.808 1.818 3.232 6.061 8.081 10.909 12.121 13.535

 


   +/\ A48

1 5 14 30 60 100 154 214 341 395 435 465 481 490 494 495

   100*8{.(]%{:)+/\ A48

0.20202 1.0101 2.82828 6.06061 12.1212 20.202 31.1111 43.2323

   100*8{.(]%{.)+/\. A48

43.2323 31.1111 20.202 12.1212 6.06061 2.82828 1.0101 0.20202

   0.005 0.025 ([>:/]%{:@]) +/\ T

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

   0.005 0.025 ([:+/"1[>:/]%{:@]) +/\ A48

1 2

   0.005 0.025 left A48

1 2

   0.005 0.025 right A48

16 16 15

   0.005 0.025 both A48

1 2 15 16

   1 2 15 16{ {.A48

7 8 21 22

   ('',0.005 0025 percent_m T48);T48

0 0  0  0

6 6  9  9

6 6 11 11

6 6 13 13

6 7 14 15

6 7 16 17

7 8 17 18

7 8 19 21

7 8 21 22

6 7 8  9 10 11 12 13 14 15 16 17 18 19 20 21 22

1 2 2  0  0  0  0  0  0  0  0  0  0  0  0  0  0

1 4 5  4  1  0  0  0  0  0  0  0  0  0  0  0  0

1 4 7  8  9  4  2  0  0  0  0  0  0  0  0  0  0

1 4 9 12 18 12  9  4  1  0  0  0  0  0  0  0  0

1 4 9 14 22 22 21 16  9  4  2  0  0  0  0  0  0

1 4 9 16 26 32 34 32 26 16  9  4  1  0  0  0  0

1 4 9 16 28 36 44 46 46 36 29 18 11  4  2  0  0

1 4 9 16 30 40 54 60 67 60 54 40 30 16  9  4  1

 


§7 2×2分割表に対するフィッシャーの精密検定

 一般に

Z =

a

b

c

d

のような2×2分割表のデータが与えられた場合、独立性の検定には「カイ2乗統計量」

を用いればよいとされている。この統計量は、a,b,c,d といった値がある程度大きければ、自由度1のカイ2乗分布で近似できる。ところがデータ数がかなり小さいときには、この近似が必ずしもよくないので、正確な確率の値を求めて検定を行なうべきであると、

フィッシャーが提案した。つまり、Zのような分割表を得る確率は、周辺度数を固定して考えるとき

という「超幾何分布」によって与えられる。例えば

 

B =

1

6

4

2

 

といったデータが与えられていれば、

という確率になる。ところで与えられたデータBと同じ周辺度数をもつ分割表は、

A

 

B

 

C

 

D

 

E

 

F

0

7

1

6

2

5

3

4

4

3

5

2

5

1

4

2

3

3

2

4

1

5

0

6

のように6通り考えられ、それぞれの確率は

である。分類項目の間が独立であれば、Aという結果のほうがBよりさらに極端なケースと考えられるから、AまたはBの起きる確率は

     Pr{A}+Pr{B} = 2/429 + 35/429 = 37/429 = 0.0862471

となる。5%の有意水準では、独立性の仮説を棄却できない。

 


【Jによるプログラム】

change=:(2 2$1 _1 _1 1)"_+]

translate=.4 :'change`(change^:_1)@.({.^:2=[:<./,)^:x. y.'

extreme=:([:i.1:+[:<./,)translate"0 2]

probc=:([:*/[:!+/,+/"1)%[:*/[:![:(],+/),

probe=:[:+/[:probc"_1 extreme

NB. lower 5% point

box=:[:<"_1](-,.])[:}:i.

matrix=:[:>L:1[:{L:2box@{.;[:<[: .L:0 box@{:

matrix_5=:[:((0.05"_>:L:0[:probe L:0])#L:0)matrix

eliminate=:[:(((<0 2$'')"_〜:[:{."1])#])matrix

lower_5=:[:>./"_1[:>[:{."1 L:0 eliminate

index_a=:[: .(1:+4:>:])}.[:>:i.

lower_table=:<,.](];"_1[:<&lower_5"_1,.)index_a

 

   ]B=:2 2$1 6 4 2                change B

1 6                            ┃2 5

4 2                            ┃3 3

   ]A=:change^:_1 B               1 tranlate B

0 7                            ┃0 7

5 1                            ┃5 1

   extreme B                      probc B

1 6                            ┃0.0815851

2 4                               probc 1 tranlate B

                               ┃0.004662

0 7                               prob B

5 1                            ┃0.0862471

 

   box 3                          box 4

┌─-┬─-┐                   ┃┌─-┬─-┬─-┐

3 0│2 1│                   ┃│4 0│3 1│2 2│

└─-┴─-┘                   ┃└─-┴─-┴─-┘


   matrix 4 3                     matrix 4 4

┌─-┬─-┐                    ┃┌─-┬─-┬─-┐

4 0│4 0│                   ┃│4 0│4 0│4 0│

0 3│1 2│                   ┃│0 4│1 3│2 2│

├─-┼─-┤                   ┃├─-┼─-┼─-┤

3 1│3 1│                   ┃│3 1│3 1│3 1│

0 3│1 2│                   ┃│0 4│1 3│2 2│

├─-┼─-┤                   ┃├─-┼─-┼─-┤

2 2│2 0│                   ┃│2 2│2 2│2 2│

0 3│1 2│                   ┃│0 4│1 3│2 2│

└─-┴─-┘                   ┃└─-┴─-┴─-┘

   matrix_5 4 3                   matrix_5 4 4

┌─-┬─┐                       ┃┌─-┬─┬─┐

4 0│                         ┃│4 0│   

0 3│                         ┃│0 4│   

├─-┼─┤                       ┃├-−┼─┼−┤

├─-┼─┤                       ┃├-−┼─┼−┤

└─-┴─┘                       ┃└─-┴─┴─┘

   $ L:0 matrix_5 4 3           ┃   $ L:0 matrix_5 4 4

┌─-┬─-┐                      ┃┌─-┬─-┬─-┐

2 2│2 0│                      ┃│2 2│2 0│2 0│

├─-┼─-┤                      ┃├─-┼─-┼─-┤

2 0│2 0│                      ┃│2 0│2 0│2 0│

├─-┼─-┤                      ┃├─-┼─-┼─-┤

2 0│2 0│                      ┃│2 0│2 0│2 0│

└─-┴─-┘                      ┃└─-┴─-┴─-┘

   eliminate 4 3               ┃   eliminate 4 4

┌─-┬─┐                       ┃┌─-┬─┬─┐

4 0│                         ┃│4 0│   

0 3│                         ┃│0 4│   

└─-┴─┘                       ┃└─-┴─┴─┘


   eliminate 5 4               ┃   eliminate 5 5

┌─-┬─-┬─┐                  ┃┌─-┬─-┬─┬─┐

5 0│5 0│                    ┃│5 0│5 0│   

0 4│0 3│                    ┃│0 5│1 4│   

├─-┼─-┼─┤                  ┃├─-┼─-┼─┼─┤

4 1│                        ┃│4 1│      

0 4│                        ┃│0 5│      

└─-┴─-┴─┘                   ┃└─-┴─-┴─┴─┘

   lower_5 4 3                      lower_5 4

4 0                              ┃4 0

   lower_5 5 4                   ┃   lower_5 5 5

5 1                              ┃5 1

4 0                              ┃4 0

 

   (index_a 3);(index_a 4); index_a

-┬─-┬───-┐

3│4 3│5 4 3 2│

-┴─-┴───-┘

   lower_table 3                    lower_table 4

-┬-┬─-┐                      ┃┌-┬-┬─-┐

3│3│3 0│                      ┃│4│4│4 2│

-┴-┴─-┘                      ┃├-┼-┼─-┤

                                 ┃│4│3│4 2│

                                 ┃└-┴-┴─-┘

 

  , lower_table 5

-┬-┬─-┬-┬-┬─-┬-┬-┬─-┬-┬-┬─-┐

5│5│5 1│5│4│5 1│5│3│5 0│5│2│5 0│

│ │ │4 0│ │ │4 0│ │ │   │ │ │  

-┴-┴─-┴-┴-┴─-┴-┴-┴─-┴-┴-┴─-┘

   , lower_table 6

-┬-┬─-┬-┬-┬─-┬-┬-┬─-┬-┬-┬─-┬-┬-┬─-┐

6│6│6 2│6│5│6 1│6│4│6 1│6│3│6 0│6│2│6 0│

│ │ │5 1│ │ │5 0│ │ │5 0│ │ │5 0│ │ │  

│ │ │4 0│ │ │4 0│ │ │   │ │ │   │ │ │  

-┴-┴─-┴-┴-┴─-┴-┴-┴─-┴-┴-┴─-┴-┴-┴─-┘


§8 メディアン(中央値)検定

 §7では、2×2分割表に対する正確な確率の値を求めて独立性の検定を行なう方法を示した。この検定法を応用して、平均など位置母数の差異を検出する「メディアン検定」が考えられる。

 §5で考えたデータ

   X=:0.1 0.2 0.2 0.4 1.3 2.4 9.3 10.8 19.9 32.5

   Y=:7.6 8.4 23.5 28.7 30.3 44.7 49.1 53.2 85.2 97.3 105.1 113

に対して、全体を大きさの順に並べたときのメディアン(中央値)は、11、12番目のデータの平均値で「21.7」のように与えられる。この値より大きいXのデータは1個だけで、残り9個は小さい。またYのほうのデータについては、メディアンより小さいデータが2個で、残り10個はメディアンより大きな値になっている。そこで

 

メディアンより大きなデータ数

1

10

メディアンより小さいデータ数

9

2

といった2×2分割表が得られる。2組の標本が同一母集団からのものであれば、各マスに入るデータの個数はほぼ均等になる傾向を示すが、平均などがずれていれば、左右いずれかの対角部分に多くのデータが入るため、より偏った分割表が得られることになる。したがって、上記のような分割表の偏りの程度によって、2組の標本が同一母集団からのものであるか否かが識別できることになる。

 

   ([:(>.,<.)[:-:[:<:#)X,Y        median X,Y

11 10                          ┃21.7

   21.7 divide X                  21.7 divide Y

1 9                            ┃10 2

   X med_table Y                  probe X med_table Y

1 10                           ┃0.000592608

9  2

であるから、X,Yに基づく分割表はかなり偏っていて、等平均の仮説は有意水準1%でも棄却できる。さらに、力士の体重に対するデータに対しても

   WTW med_table WAT              WTW med_test WAT

 6 35                          ┃6.4878e_6

37  5

   W1 med_table W2                W1 med_test W2

2 8                            ┃0.0115071

8 2

といった結果が得られる。10個という少数データの場合でも、ほぼ有意水準1%で等平均の仮説が棄却できる。「メディアン検定」は、存外、検出力の高いことが分かる。

 

【Jによるプログラム】

median=:[:(+/%#)([:(>.,<.)[:-:[:<:#)

divide=:[:({:,-/)#@],[:+/<

med_table=:[: :([:median,)divide&>;

med_test=:[:probe med_table NB. median test

 

   §9 各種検定法の性能(ロバストネス)の比較

 

 これまで紹介してきたノンパラメトリック検定法が、どの程度の識別力があるかを検討するために、乱数を使ったシミュレーション実験を行なってみよう。

 そのためにまず、分散が同じで平均だけが異なる10個の乱数の対と、さらに平均と分散も異なるデータの対を生成する関数を次のように定義する。

   unif10=.(_49.5"_+[:?10"_$100"_)%28.866"1_

      NB. uniform random number

   normal10=:_6:+[:+/([:?(12 10$100)"_)%100"_

     NB. normal random number

   dexp10=:[:-/0.7071"_*[:^.(1:+[:?(2 10$100)"_)%100"_

     NB. double-exponential random number

   make_data=:[:>([:<[:".[),.[:<"1(1:,{:@])*/{.@]+[:".[

     NB. pair of diff. means and pair of diff. means &variences

 

rand10」が一様分布、「normal10」が正規分布、そして「exp10」が両側指数分布から10個のデータを取り出す関数で、いずれも平均が0、分散が1である。

 例えば、「unif10」に「make_data」という関数を実行してみると

   'unif10 1'make_data 1 1.5

0.087 0.745 _1.299 1.230 _0.052 _0.918  0.156  0.918 _1.438 0.121

1.121 2.126  1.814 2.576 _0.230  0.706 _0.680 _0.611  0.013 2.403

 

0.087 0.745 _1.299 1.230 _0.052 _0.918  0.156  0.918 _1.438 0.121

1.682 3.241  2.765 3.864 _0.345  1.058 _1.020 _0.916  0.019 3.605

2つのテーブルの最初の行が平均0、分散1の分布からのデータで、最初のテーブルの最後の行が平均を1(右引数の先頭の要素)にしたときのデータで、さらにこのデータを1.5(右引数の最後の要素)倍したものが、最後のテーブルの最後の行に与えられる。

 さらに、右引数にテーブル(ランク2)を入力すると

   $ U=:'unif10 1'make_data"1(10 2$1 1.5)

10 2 2 10

からも分かるように、上記のような対の2組のデータが10セット(ランク4)生成される。 さらに、このようなUに順位和を算出する関数を適用してみると

   (82>:TEST);TEST=:({.ranksum{:)"2 U

┌─-┬───-┐

0 1│ 85  79│

1 1│ 70  62│

0 0│ 83  83│

1 1│ 74  75│

0 0│102 102│

1 1│ 61  63│

1 1│ 67  66│

1 1│ 81  71│

1 1│ 78  75│

1 1│ 82  68│

└─-┴───-┘

そこで、上で与えたような10組の実験で、何回、等平均の仮説が棄却されたかは

   rankt=:[:+/82>:({.ranksum{:)"2

といった関数を定義してやれば

   rankt U

7 8

のようにして算出できる。さらに、連の個数による検定、メディアン検定、さらにt−検定等に対しても

   runt=:[:+/6:>:({.nuber{:)"2

   medt=:[:+/0.05"_>:({.med_test{:)"2

   var_u=:[:(+/%<:@)[:*:]-mean=:+/%#

   tstat=:2.236"_*[:-/[:mean&>;)%[:%:var_u

   tet=:[:+/_1.734"_>:({.tstat{:)"2

のように定義して、Uに適用してみると


   ([:>rankt;runt;medt;tet)U

7 8

3 5

2 5

9 9

といった結果が得られる。

 このような実験を1000回行なう関数を次のように定義する

   simulation=:[:([:>rankt;runt;medt;tet)[make_data"1(1000 2)"_$]

 

   ]U1=:'unif10 1'simulation 1 1.5

614 720

179 281

186 303

653 792

以下同じような実験を9回行なって、結果をU2〜U10という変数に代入すると

   U1;U2;U3;U4;U5

┌───-┬───-┬───-┬───-┬───-┐

614 720│623 737│618 730│648 738│601 701│

179 281│189 299│175 287│178 254│179 254│

186 303│204 295│183 290│171 268│171 268│

653 792│659 792│672 784│672 790│649 767│

└───-┴───-┴───-┴───-┴───-┘

   U6;U7;U8;U9;U10

┌───-┬───-┬───-┬───-┬───-┐

626 737│619 741│611 715│637 745│589 713│

188 289│168 271│181 301│191 281│173 257│

178 288│184 272│185 278│202 299│177 267│

660 786│685 811│650 781│687 814│627 773│

└───-┴───-┴───-┴───-┴───-┘

のようになる。10回分を合計して

   +/U1,U2,U3,U4,U5,U6,U7,U8,U9,:U10

6186 7277

1801 2820

1863 2861

6614 7890

となる。

 さらに、正規分布や両側指数分布のデータに対しても同じ実験を行なうと

   +/N1,N2,N3,N4,N5,N6,N7,N8,N9,:N10

6460 7517

1637 2546

3008 4130

6554 7672

   +/E1,E2,E3,E4,E5,E6,E7,E8,E9,:E10

7657 8650

2987 4454

5273 6695

7249 8330

各結果は、上の行から順位和検定、連の個数による検定、メディアン検定、そしてt−検定に対するものである。総じて連の個数による検定が最も悪く、次いでメディアン検定も良くない。ところが順位和検定は、t−検定と比較しても遜色ないことが分かる。特に、両側指数分布からのデータに対しては、t−検定よりも良い結果になっている。

 さらに

   'unif10 1'simulation 1.5 1.5

907 964

441 647

473 656

940 987

   'normal10 1'simulation 1.5 1.5

930 976

449 615

647 788

939 979

   'dexp10 1'simulation 1.5 1.5

956 983

604 768

815 912

953 983

のような結果が得られる。平均で1.5位離れていれば、ほぼ帰無仮説が棄却される。何故か、分散が大きく異なっている場合のほうが検出力が高い。

 なお、2組の分布が異なる場合の実験結果も示してみると

   (rankt,runt,medt,tet)(unifl10,:1:+dexp10)"0 i.1000

731 256 461 692

   (rankt,runt,medt,tet)(normal10,:1:+dexp10)"0 i.1000

731 259 461 723

といった結果が得られる。順位和検定がややt−検定を上まわるかなという感じである。

 さらに平均が1.5離れている場合の実験結果では、順位和検定やt−検定では、ほぼ帰無仮説が棄却される。何故か、分散が大きく異なっている場合のほうが検出力が高い。また、2組の分布が異なる場合についても実験してみると、順位和検定がややt−検定を上まわる感じである。

 最後に、置換検定と順位和検定の比較をしてみるために、5個のデータを生成するための関数を

   unif5=:(_49.5"_+[:?(5$100"_)%28.866"1_

   normal5=:_6:+[:+/([:?(12 5$100)"_)%100"_

   dexp5=:[:-/0.7071"_*[:^.(1:+[:?(2 5$100)"_)%100"_

と定義し、さらに順位和検定と置換検定の実験結果を与える関数も次のように与える。

   rankt55=:[:+/19"_>:({.ramksum{:)"2

   C55=:comball 5 5

   permt55=:[]+/12.6"_>({.([:+/"1+/@[>:+/"1 C55"_#,){:)"2

 1000回のシミュレーションの結果は以下のようになる。

   (permt55,:rankt55)'unif5 1'make_data"1(1000 2$1 1.5)

368 471

357 459

   (permt55,:rankt55)'unif5 1'make_data"1(1000 2$1.5 1.5)

686 814

643 776

   (permt55,:rankt55)'normal5 1'make_data"1(1000 2$1 1.5)

432 508

422 487

   (permt55,:rankt55)'normal5 1'make_data"1(1000 2$1.5 1.5)

658 774

649 736

   (permt55,:rankt55)'dexp5 1'make_data"1(1000 2$1 1.5)

510 611

512 601

   (permt55,:rankt55)'dexp5 1'make_data"1(1000 2$1.5 1.5)

762 853

735 828

 総じて、置換検定のほうが順位和検定よりわずかながらよい。


以上の実験結果から、次のような結論が得られる。

 

分布が比較的かたまった状態で分離している場合にはt−検定が良い。

分布が正規分布より散らばっているケースでは順位和検定が良い。

また、分布形が異なる感じで分離しているときも順位和検定が良い。

データが少なく、間隔尺度なら置換検定が良いが、順位和検定でも構わない。

連の個数やメディアン検定を積極的に使う理由はなさそうである。

 

§10 逐次差の符号によるランダムネスの検定

 ある講座での出席者数を1回目から10回目まで数えてみたら

 

講義回数

1  2  3  4  5  6  7  8  9  10

出席者数

128 113 110 100 108 101 104  98  92  82

 

のような結果が得られたとする。これより、出席者数は減少傾向にあるといえるだろうか? 隣あった回数の間で、出席者数が減少していれば「0」、増加していれば「1」をつけてみると、(4,5)と(6,7)のところだけが1で、残りの7か所では0である。

 一般に、

という観測系列が与えられたときに

といった統計量を定義すると、増加傾向にある個数を示す値である。

 Y(n)の期待値と分散は

     E{Y(n)} = (n-1)/2 , V{Y(n)} = (n+1)/12

のように与えられる。さらに、を{1,2, ………, n}の置換の中でY(n)=kとなる組合せの数とすると、

といった漸化関係の成り立つことが分かる。n個の観測系列がランダムであるという仮説の下では、確率の値が


 


のように与えられる。したがって

     Pr{Y(n)≦y(n,α)} ≦ α

となるy(n,α)が表1のように与えられる。

例題の場合には、X(10)=2であるから、有意水準 2.5%でランダムネスの仮説が棄却でき、出席者数は減少傾向にあると判断できる。

1 の表

n

1% 2.5% 5%

6

7

8

9

10

11

12

13

14

15

  0   0   0

 0   1   1

1   1   1

1   1   2

1   2   2

2   2   2

2   3   3

3   3   3

3   3   4

3   4   4

【Jによるプログラム】

   stat_line=:[:+/{:<{.

   more=:[:+/([:(,.|.)1:+i.@>:@#)*[:>2:<;.0\0:,]0:

   point_line=:4 :'<:+/"1 x.>:/+/\(more^:(y.-1)1)%!y.'

 

   X=:128 113 110 100 108 101 104 98 92 82

   (}:<}.)X

0 0 0 1 0 1 0 0 0

   stat_line X

2

   (1:+i.@>:@#)1 1                (0:,]0:)1 1

1 2 3                          ┃0 1 1 0

   ([:(,.|.)1:+i.@>:@#)1 1        ([:>2:<;.0\0:,]0:)1 1

1 3                            ┃1 0

2 2                            ┃1 1

3 1                            ┃0 1

   (([:(,.|.)1:+i.@>:@#)*[:>2:<;.0\0:,]0:)1 1

1 0

2 2

0 1

 


   more 1 1                       more^:2(1 1)

1 4 1                          ┃1 11 11 1

   more^:3(1 1)                   more^:4(1 1)

1 26 66 26 1                   ┃1 57 302 302 57 1

   more^:5(1)                     +/more^:5(1)

1 57 302 302 57 1              ┃720

 

   (i.10),:more^:9(1)

0    1     2      3       4       5      6     7    8 9

1 1013 47840 455192 1310354 1310354 455192 47840 1013 1

   +/more^:9(1)                   !10

3628800                        ┃3.6288e6

 

   0.01 0.025 0.05 point_line 10

1 2 2

 

 §11 全ての対の符号によるランダムネスの検定

 観測系列のi<jというペア(Xi,Xj)に対して、Xi<Xjなら「1」、Xi≧Xjなら「0」を与えて、

    

といった統計量が定義できて、これまた増加傾向にある個数を示す値で、0から nC2=

(n−1)/2の範囲の値をとり、Z(n)の期待値と分散は

     E{Z(n)} = n(n-1)/4 , V{Z(n)} = n(n-1)(2n+5)/72

のように与えられる。さらに、を{1,2,………,n}の置換の中でZ(n)=kとなる組合せの数とすると、

    

 例えばn=4の場合、からまでの値が{1,3,5,6,5,3,1}だから

     Pr{Z(4)=0}=Pr{Z(4)=6}= 1/24 ,Pr{Z(4)=1}=Pr{Z(4)=5}= 3/24

     Pr{Z(4)=2}=Pr{Z(4)=4}= 5/24 ,Pr{Z(4)=3}= 6/24

といった具合に算出される。

 一般に

     Pr{Z(n)≦p(n,α)} ≦ α

となるようなp(n,α)が表2のように与えられる。

 さて与えられた例題の場合には、表2に示したように、1の個数が4であるから、有意水準1%でもランダムネスの仮説が棄却できる。

表2 ランダムネスの検定のパーセント点

 n 1% 2.5% 5%

n 1% 2.5% 5%

n 1% 2.5% 5%

5  0   0   1

6  1   1   2

7  2   3   4

8  4   5   6

9  6   8   9

10  9   11  12

11  12  14  16

12  15  18  20

13  19  22  25

14  24  27  29

15  28  32  35

16  34  37  41

17  39  43  47

18  45  50  54

19  52  57  61

20  59  64  69

21  66  72  78

22  74  80  85

23  82  89  94

24  91  98  104

25 100  107  114

26 109  117  124

27 119  128  135

28 130  139  146

29 140  150  158

30 152  162  170

31 164  174  183

32 176  187  196

33 188  200  210

34 202  214  224

35 215  228  239

36 229  242  254

37 244  257  269

38 259  273  285

39 274  289  301

40 290  305  318

45 376  394  410

50 473  495  513

60 702  731  755

 

【Jによるプログラム】

   stat_pair=:[:+/([:+/{.<}.)\.

   distribution=:3 :'[:+/"1[:=[:stst_pair"1([:i.!)A.i.

   next=:1:+[:-:1:+[:%:(8:*#)-7:

   index=:[:>.[:-:1:+[:-:[:(*<:)next

   take=:[:+/\(index{.]

   modify=:3 :'(t{.u),-/((-,])(#u)-t=.(#<.next)y.){."0 1 u=.take y.'

   more_pair=:3 :'h,|.((1+-:(*<:)next y.)-#h){.h=.modify y.'

   point_pair=:4 :'<:+/"1 x.>:/(+/\more_pair^:(y.-1)1)%!y.'

   var_pair=:((5:+2:*])*]*<:)%72"_

   forth=:3 :'(y.*(y.-1)*_372 _997 _127 328 100&p.y.)%43200'"0

   range_pair=:([:(]--:@{:)[:i.1:+[:-:]*<:)%[:%:var_pair

   nden=:([:^[:-*:%2:)%[:%:[:o.2:

   middle=:[:}.]*(0.002*i.500)"_

   ndfs=:]*([:+/([:-:[:nden 0:,]),[:nden middle)%500"_

   ndf=:ndfs+0.5"_

   orgf=:nden*0 _3 0 1&p.

   approximation=:ndf@[-orgf@[*(_3:+forth%[:*:var_pair)@]%24"_

   point_ap=:[:<:[:+/"1[>:/[:(range_pair approximation"0])]

 

   ]Y=./:^:2 X.

9 8 7 3 6 4 5 2 1 0

   ]\.Y                          ({.<}.)"1]\.Y

9 8 7 3 6 4 5 2 1 0            ┃0 0 0 0 0 0 0 0 0

8 7 3 6 4 5 2 1 0 0            ┃0 0 0 0 0 0 0 0 0

7 3 6 4 5 2 1 0 0 0            ┃0 0 0 0 0 0 0 0 0

3 6 4 5 2 1 0 0 0 0            ┃1 1 1 0 0 0 0 0 0

6 4 5 2 1 0 0 0 0 0            ┃0 0 0 0 0 0 0 0 0

4 5 2 1 0 0 0 0 0 0            ┃1 0 0 0 0 0 0 0 0

5 2 1 0 0 0 0 0 0 0            ┃0 0 0 0 0 0 0 0 0

2 1 0 0 0 0 0 0 0 0            ┃0 0 0 0 0 0 0 0 0

1 0 0 0 0 0 0 0 0 0            ┃0 0 0 0 0 0 0 0 0

0 0 0 0 0 0 0 0 0 0            ┃0 0 0 0 0 0 0 0 0

 

   ([:+/{.<}.)\.) X              stat_pair X

0 0 0 3 0 1 0 0 0 0            ┃4

 

   (([:i.!)A.i.)3

0 1 2

0 2 1

1 0 2

1 2 0

2 0 1

2 1 0

   ([:stat_pair"1([:i.!)A.i.)3

3 2 2 1 1 0

   ([:=[:stat_pair"1([:i.!)A.i.)3

1 0 0 0 0 0

0 1 1 0 0 0

0 0 0 1 1 0

0 0 0 0 0 1

 

   ]D3=.distribution 3            ]D4=.distribution 4

1 2 2 1                        ┃1 3 5 6 5 3 1

   ]D5=.distribution 5

1 4 5 9 15 20 22 20 15 9 4 1

 

   6!:2'distribution 6'           6!:2'distribution 7'

1.15                           ┃12.13

   6!:2'distribution 8'           6!:2'distribution 9'

98.71                          ┃out of memory

【直接計算するのは時間がかかりすぎる!】

 

   next L:0 D3;D4;D5              index L:0 D3;D4;D5

-┬-┬-┐                    ┃┌-┬-┬-┐

4│5│6│                    ┃│4│6│8│

-┴-┴-┘                    ┃└-┴-┴-┘

   take L:0 D3;D4;D5

1 3 5 6

1 4 9 15 20 23

1 5 14 29 49 71 91 106

   (3 :'((-,])(#u)-(#<.next)y.){."0 1 u=.take y.')D4

23

 1

   (3 :'((-,])(#u)-(#<.next)y.){."0 1 u=.take y.')D5

91 106

 1   5

   modify L:0 D3;D4;D5

1 3 5 6

1 4 9 15 20 22

1 5 14 29 49 71 90 101

   D4 ; more_pair D3

1 3 5 6 5 3 1

1 3 5 6 5 3 1

   D5 ; more_pair D4

1 4 5 9 15 20 22 20 15 9 4 1

1 4 5 9 15 20 22 20 15 9 4 1

   ]D6=.more_pair D5

1 5 14 29 49 71 90 101 101 90 71 49 29 14 5 1

   ($D6);(1+-:6*5);(+/D6);!6

16

16

720

720

   D4 ; more_pair^:3(1)

1 3 5 6 5 3 1

1 3 5 6 5 3 1

   D5 ; more_pair^:4(1)

1 4 5 9 15 20 22 20 15 9 4 1

1 4 5 9 15 20 22 20 15 9 4 1

 

   0.01 0.025 0.05 point_pair 5┃   0.01 0.025 0.05 point_pair 6

0 0 1                          ┃1 1 2

   0.01 0.025 0.05 point_pair 7┃   0.01 0.025 0.05 point_pair 8

2 3 4                          ┃4 5 6

   0.01 0.025 0.05 point_pair 9

6 8 9

 

   §12 順位相関係数

 一般に、対で観測されたn個の個体の、各変量ごとの順位を(),()として、これらを通常の計量値のように考えて相関係数を求めてみると

のようになる。これは、「スピアマンの順位相関係数」と呼ばれているもので、2変量間に相関がない場合には

     E{R(n)} = 0 , V{R(n)} = 1/(n-1)

となる。さらにnが30くらい大きければ、近似的に正規分布に従う。

 2組のデータが全く同順位の場合には、順位の差の平方和は

であるから、R(n)=1となり、反対に全く逆順のときには

であるから、R(n)=−1となる。

 さて、表1のデータに対する順位の差の平方和は

と算出されるから、D≦14となる順列の総数は表2に示したように167で、また1から7までの数字の順列の総数は7!=5040であるから

のように与えられる。したがって、無相関であるという仮説を有意水準5%で棄却でき、国語と数学の成績の間には「正の相関がある」と結論できる。

 


   順位相関係数の計算

 

   rank=:i.@#(=/+/ .*[)/:

   stst_sqd=:[:+/\[:*:i.@#@]-rank@]/:[

   stst_cor=:1:-6:*stat_sqd%[:(]**:@]-1:)#@]

 

   seqall=:([:i.[:!#)A.] NB. all permutations

   sqdall=:[:/:〜[:+/"1[:*:]-"1 seqall

       NB. sum of quadratic difference

   sqdtable=:[:(〜.,:[:+/"1=)sqdall

   index=:_1:+[:+/([*[:!])>:[:+/\]){:sqdtable@]

   lower=:([*[:!])index[:sqdtable[:i.]

   cvalue=:1:-6:*lower%[:(]**:@]-1:

       NB. critical value of rank correlation

 

   ]RX=.rank X=:50 60 65 75 35 80 70

1 2 3 5 0 6 4

   ]RY=.rank Y=:40 35 80 55 25 70 45

2 1 6 4 0 5 3

   *:RX-RY

1 1 9 1 0 1 0

   +/*:RX-RY                      X stat_sqd Y

14                             ┃14

   X stat_cor Y

0.75

 

   seqall i.3                     sqdall i.3

0 1 2                          ┃0 2 2 6 6 8

0 2 1                             ]R3=.sqdtable i.3

1 0 2                          ┃0 2 6 8

1 2 0                          ┃1 2 2 1

2 0 1

2 1 0

 

   sqdall i.4

0 2 2 6 6 8 2 4 6 12 10 14 6 10 8 14 16 18 12 14 14 18 18 20

   ]D4=.sqdtable i.4

0 2 4 6 8 10 12 14 16 18 20

1 3 1 4 2  2  2  4  1  3  1

   ]D5=.sqdtable i.5

0 2 4 6 8 10 12 14 16 18 20 22 24 26 28 30 32 34 36 38 40

1 4 3 6 7  6  4 10  6 10  6 10  6 10  4  6  7  6  3  4  1

   +/"1 D5

420 120

   16{."1 D6=.sqdtable 6

0 2 4 6  8 10 12 14 16 18 20 22 24 26 28 30

1 5 6 9 16 12 14 24 20 21 23 28 24 34 20 32

   $ D6

2 36

   +/"1 D6

1260 720

   29{. {: D7=.sqdtable 7

1 6 10 14 29 26 35 46 55 54 74 70 84 90 78 90 123 134 147 98 168 130 175 144 168 144 168 144 184

   +/"1 D7

3192 5040

   43{. {: D8=.sqdtable 8

1 7 15 22 47 54 70 94 126 124 178 183 237 238 276 264 379 349 380 400 517 394 542 492 640 557 666 595 776 684 786 718 922 745 917 781 982 826 950 844 1066 845 936

   $(([:i.!)A."0 1 i.) 9

out of memory

 

   (0.05*!3)index sqdtable i.3 ┃   (0.05*!4)index sqdtable i.4

0                              ┃0

   (0.05*!5)index sqdtable i.5 ┃   (0.05*!6)index sqdtable i.6

2                              ┃6

 

   0.05 lower 5                   0.01 lower 5

2                              ┃0

   0.05 lower 6                   0.01 lower 6

6                              ┃2

 

   0.05 cvalue 5                  0.05 cvalue 6

0.9                            ┃0.828571

   0.05 cvalue 7                  0.01 cvalue 7

0.714286                       ┃0.892857

 

   §13 順位相関係数の標本分布とパーセント点の導出

 

   seqsum=:4 :0

'L R'=.((i.!#x.)A.x.);(i.!#y.)A.y.

r=.((k=.0){L),"1 R

while. k<<:#L

  do. r=.r,((k=.k+1){L),"1 R

end.

)

   sqd11=:4 :0

'L R'=.(t;-.t=.((<.-:y)=+/"1 t)#t=.#:i.-2^y=.#y.)#L:0 y.

A=.〜.r [ B=.+/"1=r=./:〜+/"1*:x.-"1(k{L)seqsum(k=.0){R

while. k<<:#L

  do. r=./:〜+/"1*:x.-"1(k{L)seqsum(k=.k+1){R

      B=.+/"1((A=/〜.r)#+/"1=r),.((A=./:〜〜.A,r)=/A)#B

end.

A,:B

)

   sqd00=:4 :0

'a b'=.(d=.y.-x.)({.;}.)i.y.

'A B'=.(s;-.s=.(d=+/"1 t)#t=.#:i.-2^y.)#L:0 i.y.

C=.a sqd11 k{A [ D=.b sqd11(k=.0){B

F=.+/,.((E=./:〜〜.,s)=/"1 s=.({.C)+/{.D)#"_1({:C)*/{:D

while. k<<:#A

  do. C=.a sqd11 k{A [ D=.b sqd11(k=.k+1){B

      F=.((E=./:〜〜.E,t)=/E)#F [ t=./:〜〜.,s=.({.C)+/{.D

      F=.+/"1 F,.(E=/t)#+/,.(t=/"1 s)#"_1({:C)*/{:D

end.

E,:F

)

sqd_p=:4 :'(<:+/t<{.x.){"1({.y.),:t=.+/\({:y.)%!{:x.'

 

【演算時間】

 

   6!:2'(i.k)sqd11 1;i.k=7'       6!:2'(i.k)sqd11 0;i.k=8'

0.05(0.49)                     ┃0.33(3.4)

   6!:2'(i.k)sqd11 1;i.k=9'       6!:2'(i.k)sqd11 2;i.k=10'

2.69(33.09)                    ┃140.55(274.52)

   6!:2'(i.k)sqd11 2;i.k=11'      6!:2'(i.k)sqd11 2;i.k=12'

1605.42(27323.5)              

 

   6!:2'D9=.5 sqd00 9'            6!:2'D10=.5 sqd00 10'

2.14                           ┃7.3

   6!:2'D11=.6 sqd00 11'          6!:2'D12=.6 sqd00 12'

27.02                          ┃102.1

   6!:2'D12=.7 sqd00 12'          6!:2'D12=.8 sqd00 12'

94.96                          ┃199.49

   6!:2'D12=.7 sqd00 13'          6!:2'D13=.8 sqd00 13'

335.92                         ┃1377.48

   6!:2'D14=.7 sqd00 14'          6!:2'D15=.8 sqd00 15'

1157.5                         ┃5270.49

   6!:2'D16=.8 sqd00 16'          6!:2'D17=.8 sqd00 17'

18347                          ┃123318

 

   0.05 7 sqd_p D7                0.05 8 sqd_p D8

16 0.0440476                   ┃30 0.0480903

   0.025 7 sqd_p D7               0.025 8 sqd_p D8

12 0.0240079                   ┃22 0.0229167

   0.01 7 sqd_p D7                0.01 8 sqd_p D8

 

6 0.00615079                   ┃14 0.00768849

   0.05 9 sqd_p D9                0.05 10 sqd_p D10

48 0.0480903                   ┃72 0.0481393

   0.025 9 sqd_p D9               0.025 10 sqd_p D10

36 0.021627                    ┃58 0.0244894

   0.01 9 sqd_p D9                0.01 10 sqd_p D10

26 0.00861166                  ┃42 0.00870288

 

   0.05 11 sqd_p D11              0.05 12 sqd_p D12

102 0.0469582                  ┃142 0.0494539

   0.025 11 sqd_p D11             0.025 12 sqd_p D12

84 0.0238787                   ┃118 0.0244408

   0.01 11 sqd_p D11              0.01 12 sqd_p D12

64 0.0090973                   ┃92 0.0092612

 

   0.05 13 sqd_p D13              0.05 14 sqd_p D14

188 0.0485266                  ┃244 0.0486266

   0.025 13 sqd_p D13             0.025 14 sqd_p D14

160 0.0249319                  ┃210 0.0249824

   0.01 13 sqd_p D13              0.01 14 sqd_p D14

128 0.00970532                 ┃170 0.00953338

 

   0.05 15 sqd_p D15              0.05 16 sqd_p D16

310 0.00486229                 ┃388 0.0492784

   0.025 15 sqd_p D15             0.025 14 sqd_p D16

268 0.0244026                  ┃338 0.0246591

   0.01 15 sqd_p D15              0.01 16 sqd_p D16

222 0.00972995                 ┃284 0.00999043

 

   0.05 17 sqd_p D17              0.05 18 sqd_p D18

478 0.0049834                  ┃580 0.0499117

   0.025 17 sqd_p D17             0.025 18 sqd_p D18

418 0.0245113                  ┃512 0.0249773

   0.01 17 sqd_p D17              0.01 18 sqd_p D18

354 0.00983329                 ┃436 0.00985633

 

【「D(n)」のパーセント点】

 n    5 6  7  8  9 10  11  12  13  14  15  16  17  18

5%    2 6 16 30 48 72 102 142 188 244 310 388 478 580

2.5%  0 4 12 22 36 58  84 118 160 210 268 338 418 512

1%    0 2  6 14 26 42  64  92 128 170 222 284 354 436

 


 さらに高次の項を用いて連続補正を行なった場合はより改良される。

     F(x) = Φ(x)-φ(x){c4H3(x)/24+c6/720H5(x)+c422H7(x)/1152

                     +c8H7(x)/40324+c4c6H9(x)/17280

                     +c43H11(x)/82944

のように表現できる。ここで

     c4 = 6(36-5n-19n2)/25n(n+1)(n-1)

     c6 = 48N6/245n3(n+1)3(n-1)

     c8 = 144N8/875n5(n+1)5(n-1)3

     N6 = -1800+2760n+4054n2-2637n3-2603n4+723n5+583n6

     N8 = 846720-1080567n-1616688n2+2358048n3+1800776n4

        -1690125n5-1012323n6+578442n7+304254n8-83709n9-41939n10

     H3(x) = -3x+x3

     H5(x) = 15x-10x3+x5

     H7(x) = -105x+105x3-21x5+x7

     H9(x) = 945x-1260x3+378x5-36x7+x9

     H11(x) = -10395x+17325x3-6930x5+990x7-55x9+x11

である。

 

この近似式を用いて、パーセント点を求める関数を次のように定義する。

   cum4=:_0.24"_*_36 5 19&p.%]*<:@*:

   numerator6=:_1800 2760 4054 _2637 723 583&p.

   cum6=:48r245"_*numerator6%*:@<:*(]*>:)^3:

   numerator8=:846720 1080567 _1616688 2358048 1800776 _1690125

               _1012323 578442 304254 _83709 _41939&p.

   cum8=.144r875*numerator8%(<:^3:)*(]*>:)^5:

   orgf3=:nden*0 _3 0 1&p.

   orgf5=.nden*0 15 0 _10 0 1&p.

   orgf7=.nden*0 _105 0 105 0 _21 0 1&p.

   orgf9=.nden*0 945 0 _1260 0 378 0 _36 0 1&p.

   orgf11=.nden*0 _10395 0 17325 0 _6930 0 990 0 _55 0 1&p.

   dev1=:orgf3@[*cum4@]%24"_

   dev2=:(orgf5@[*cum6@]%720"_)+orgf7@[**:@cum4@]%1152"_

   dev3=:(orgf7@[*cum8@]%40324"_)+orgf9@[*(cum6*cum4)@]%17280"_)         +orgf11@[*(cum4@[^3:)%82944"_

   index1=:_1:+[:+/[>range0@](ndf@[-dev1+dev2+dev3)"0]

   cvalue_sqd=:index1"0{range_sqd@]

 

   0.05 0.25 0.01 cvalue_sqd 7

16 10 6

   0.05 0.25 0.01 cvalue_sqd 8

30 22 14

   0.05 0.25 0.01 cvalue_sqd 9

48 36 26

   0.05 0.25 0.01 cvalue_sqd 10

72 58 42

   0.05 0.25 0.01 cvalue_sqd 11

102 84 64

   0.05 0.25 0.01 cvalue_sqd 12

142 118 92

   0.05 0.25 0.01 cvalue_sqd 13

188 160 128

   0.05 0.25 0.01 cvalue_sqd 14

244 208 170

   0.05 0.25 0.01 cvalue_sqd 15

310 268 222

   0.05 0.25 0.01 cvalue_sqd 16

388 338 282

   0.05 0.25 0.01 cvalue_sqd 17

478 418 354

   0.05 0.25 0.01 cvalue_sqd 18

580 512 436


 

PERCENTAGE POINT OF D(n) for n=4(1)30

n  5%  2.5%  1%

 n  5%  2.5%  1%

n  5%  2.5%  1%

4   0    -    -

5   2    0    0

6   6    4    2

7  16   12    6

8  30   22   14

9  48   36   26

10  72   58   42

11 102   84   64

12 142  118   92

13 188  160  128

14 244  210  170

15 310  268  222

16 388  338  284

17 478  418  354

18 580  512  436

19 694  616  530

20 824  736  636

21 970  868  756

22 1132 1018  890

23 1310 1182 1040

24 1508 1364 1206

25 1724 1566 1388

26 1958 1784 1588

27 2214 2022 1806

28 2492 2282 2044

29 2794 2564 2304

30 3118 2866 2584

 

      【参考文献】

[1]Ansari,A.R.-Bradley,R.A.(1961).Rank-sum test for dispersions.                            Ann.Math.Stat.31

[1]David,S.T.,Kendall,M.G.and Stuart,A.(1951).Some questions of distribution

in the theory of rank correlation.Biometrika 38. 131-140.

[2]Kendall,M.G.,Kendall,S.F.H. and Smith,B.B.(1938).The distribution of Speaman's coefficient of rank correlation in a universe in which all rankings occur an equal number of times.Biometrika 30. 251-273.

[3]Mann,H.B.(1945).Nonparametric tests against trend.Econometrica 13. 245-259.

[2]Moses,L.E.(1952).A two sample test. Psychometrika 17

[3]Pitman,E.J.G.(1937).Significance tests which may be applied to samples from any populations. .The correlation coefficient test. J.R.S.S.Supp.4 225-232.

[4]清水良一(1976)『中心極限定理』(教育出版)