確率・統計

(10) 順序統計量

順序尺度では、平均値や標準偏差の推定や検定を行うことはできなくなります。例えば、同学年の各クラスごとの学力に偏りがないかをしらべるとき、クラスごとの得点の分布状態を使うのではなく学年内の順位から調べなければならないとしたとき、今まで見てきた推定・検定方法は利用できなくなります。この章では、順序統計量の分布からはじめて、順序尺度を使った検定法などについて紹介したいと思います。

(注) 数式などの記法について (ドキュメントの中で使用している数式の表現方法に関する注意点です)

1) 順序統計量(Order Statistcs)

要約統計量の中で、順序尺度に対して利用できるパラメータに「中央値」などの「分位数」や「最頻値」がありました。これらは、母集団から抽出した確率変数を大きさの順で並べた上で求めることができるので、標本の並べ替え処理が必要になります。抽出した独立な確率変数 Xi ( i = 1, 2, ... N ) を大きさの順で並べ替えたものを

{X(j)} = X(1) ≤ X(2) ≤ ... ≤ X(N)

としたとき、これを「順序統計量(Order Statistcs)」といいます。ここで、母集団は確率密度関数 p(x) に従う離散分布として、その確率変数を xi ( i = 1, 2, ... ) と小文字で表すことにします。また、i ≤ k ならば xi ≤ xk が常に成り立っているとします。つまり、xi は昇順で並べ替えられているとするわけです。各確率変数 xi の確率密度は pi = p( xi ) で表し、その累積密度を Pi = Σk{1→i}( pk ) とします。

ある i に対して Xk ≤ xi となるならば Yk = 1、そうでなければ Yk = 0 と定めます。i = 1 ならば、Xk = x1 となる Xk がない限り常に Yk = 0 が成り立ち、i が大きくなるほど 0 の要素は減少していくことになります。Xk ≤ xi となる確率は、Xk が xi と等しくなるまでの累積密度に等しく、それは Pi を意味するので

P( Yk = 1 ) = P( Xk ≤ xi ) = Pi

になります。Yk は互いに独立であり、確率 Pi に従うベルヌーイ分布と考えることができるので( i は固定であり、Pi も一つに定まることに注意してください )、その和を

R = Σk{1→N}( Yk )

とすれば、R は二項分布 BN,Pi(r) = NCr Pir( 1 - Pi )N-r に従うことになります。ここで、ある確率変数 xi に対して、順序統計量 X(j) が xi 以下となる確率は、X(1) から X(j) までが xi 以下となる確率を意味し、このとき R ≥ j が成り立つことから

P( X(j) ≤ xi ) = P( R ≥ j ) = Σr{j→N}( BN,Pi(r) ) = Σr{j→N}( NCr Pir( 1 - Pi )N-r )

という結果になります。また、

p( X(j) = xi )=P( X(j) ≤ xi ) - P( X(j) ≤ xi-1 )
=Σr{j→N}( NCr Pir( 1 - Pi )N-r ) - Σr{j→N}( NCr Pi-1r( 1 - Pi-1 )N-r )
=Σr{j→N}( NCr { Pir( 1 - Pi )N-r - Pi-1r( 1 - Pi-1 )N-r } )

が成り立ちます。ただし、i = 1 の場合を考えて、P0 = 0 とします。

次に、母集団が連続分布 p(x) ( x ∈ R ) に従うとして、その累積密度を P(x) = ∫{-∞→x} p(t) dt と表します。

ある x に対して Xk ≤ x となるならば Yk = 1、そうでなければ Yk = 0 と定めます。x が限りなく小さければ Yk = 0 となり、x が大きくなるほど 0 の要素は減少していくことになります。Xk ≤ x となる確率は、Xk が x と等しくなるまでの累積密度に等しく、それは P(x) を意味するので

P( Yk = 1 ) = P( Xk ≤ x ) = P(x)

になります。Yk は互いに独立であり、確率 P(x) に従うベルヌーイ分布と考えることができるので、その和を

R = Σk{1→N}( Yk )

とすれば、R は二項分布 BN,P(x)(r) = NCr P(x)r( 1 - P(x) )N-r に従うことになります。ここで、ある確率変数 x に対して、順序統計量 X(j) が x 以下となる確率は、X(1) から X(j) までが x 以下となる確率を意味し、このとき R ≥ j が成り立つことから

P( X(j) ≤ x ) = P( R ≥ j ) = Σr{j→N}( BN,P(x)(r) ) = Σr{j→N}( NCr P(x)r( 1 - P(x) )N-r )

という結果になります。P( X(j) ≤ x ) は X(j) に対する累積密度関数を表しているので、これを x で微分すると X(j) の確率密度関数を得ることができます。これを p( X(j) = x ) と表すと、

p( X(j) = x )=(d/dx)( P( X(j) ≤ x ) )
=(d/dx)Σr{j→N}( NCr P(x)r( 1 - P(x) )N-r )
=Σr{j→N}( NCr (d/dx){ P(x)r( 1 - P(x) )N-r } )
=Σr{j→N}( NCr { r・P(x)r-1・P'(x)・( 1 - P(x) )N-r - ( N - r )・P(x)r・( 1 - P(x) )N-r-1・P'(x) } )
=Σr{j→N}( NCr・r・P(x)r-1・( 1 - P(x) )N-r - NCr・( N - r )・P(x)r・( 1 - P(x) )N-r-1 )・P'(x)
={ Σr{j→N}( { N! / ( r - 1 )!( N - r )! }・P(x)r-1( 1 - P(x) )N-r ) - Σr{j→N-1}( { N! / r!( N - r - 1 )! }・P(x)r( 1 - P(x) )N-r-1 ) }・P'(x)
={ Σr{j→N}( { N! / ( r - 1 )!( N - r )! }・P(x)r-1( 1 - P(x) )N-r ) - Σr{j+1→N}( { N! / ( r - 1 )!( N - r )! }・P(x)r-1( 1 - P(x) )N-r ) }・P'(x)
={ N! / ( j - 1 )!( N - j )! }・P(x)j-1( 1 - P(x) )N-j P'(x)
={ N! / ( j - 1 )!( N - j )! }・{ ∫{-∞→x} p(t) dt }j-1・{ ∫{x→∞} p(t) dt }N-j・p(x)

となります。連続分布も、離散分布の場合とほとんど同じ方法で確率密度関数を求めることができました。

また、抽出した N 個の標本において二個の要素 ( X(k), X(l) ) ( k < l ) がそれぞれ ( x, y )だったときの同時確率密度を p( x, y ) とすると、X(j) < x となる標本 X(j) が k - 1 回、x < X(j) < y となる標本 X(j) が l - k - 1 回、y < X(j) となる標本 X(j) が N - l 回、そして X(k), X(l) がそれぞれちょうど x, y となるときの多項分布に従うと考えれば、

p( x, y )={ N! / ( k - 1 )!( l - k - 1 )!( N - l )! }{ P(x) }k-1p(x){ P(y) - P(x) }l-k-1p(y){ 1 - P(y) }N-l
={ N! / ( k - 1 )!( l - k - 1 )!( N - l )! }{ ∫{-∞→x} p(t) dt }k-1{ ∫{x→y} p(t) dt }l-k-1{ ∫{y→∞} p(t) dt }N-lp(x)p(y)

と表すことができます。これはさらに一般化できて、K 個の要素 ( X(j1), X(j2) ... X(jK) ) が x = ( x1, x2, ... xK ) だったときの同時確率密度を p(x) とすれば

p(x) = { N! / ( j1 - 1 )!( j2 - j1 - 1 )!...( jK - jK-1 - 1 )!( N - jK )! }{ ∫{-∞→x1} p(t) dt }j1-1{ ∫{x1→x2} p(t) dt }j2-j1-1...{ ∫{xK→∞} p(t) dt }N-jKp(x1)...p(xK)

となります。


例として、次のような一様分布に対する順序統計量の確率密度関数を考えてみます。

p(x)=1 / θ( 0 ≤ x ≤ θ )
=0( その他 )

このとき、順序統計量の j 番目 X(j) に対する確率密度関数は

p( X(j) = x ) = { N! / ( j - 1 )!( N - j )! }・{ ∫{-∞→x} p(t) dt }j-1・{ ∫{x→∞} p(t) dt }N-j・p(x)

と表されますが、x < 0 かつ θ < x のときは明らかにゼロになるので、0 ≤ x ≤ θ の場合のみを考えれば

p( X(j) = x )={ N! / ( j - 1 )!( N - j )! }・{ ∫{0→x} 1 / θ dt }j-1・{ ∫{x→θ} 1 / θ dt }N-j・( 1 / θ )
={ N! / ( j - 1 )!( N - j )! }・( x / θ )j-1・( 1 - x / θ )N-j・( 1 / θ )

t = x / θ とすると、dt / dx = 1 / θ となって、変数を t に変換した確率密度関数を q(t) とすれば p(x)dx = q(t)dt なので、

q(t) = { N! / ( j - 1 )!( N - j )! }・tj-1・( 1 - t )N-j ( 0 ≤ t ≤ 1 )

になります。この期待値 E[t] は

E[t] = { N! / ( j - 1 )!( N - j )! }∫{0→1} tj・( 1 - t )N-j dt

で計算できますが、積分値は

∫{0→1} tj・( 1 - t )N-j dt=[ jtj-1・( 1 - t )N-j ]{0→1} + { j / ( N - j + 1 ) }∫{0→1} tj-1・( 1 - t )N-j+1 dt
={ j / ( N - j + 1 ) }{ [ ( j - 1 )tj-2・( 1 - t )N-j+1 ]{0→1} + { ( j - 1 ) / ( N - j + 2 ) }∫{0→1} tj-2・( 1 - t )N-j+2 dt
:
=-{ j! / ( N - j + 1 )( N - j + 2 )...N・( N + 1 ) }[ ( 1 - t )N+1 ]{0→1}
=j!( N - j )! / ( N + 1 )!

なので、

E[t] = { N! / ( j - 1 )!( N - j )! }{ j!( N - j )! / ( N + 1 )! } = j / ( N + 1 )

になります。また、分散は V[t] = E[t2] - ( E[t] )2 より

V[t] = { N! / ( j - 1 )!( N - j )! }∫{0→1} tj+1・( 1 - t )N-j dt - { j / ( N + 1 ) }2

で計算できますが、積分値は ∫{0→1} tj+1・( 1 - t )(N+1)-(j+1) dt と表せることから ( j + 1 )!( N - j )! / ( N + 2 )! と求めることができるので、

V[t]={ j( j + 1 ) / ( N + 1 )( N + 2 ) } - { j / ( N + 1 ) }2
={ j( j + 1 )( N + 1 ) - j2( N + 2 ) } / ( N + 1 )2( N + 2 )
=j( jN + j + N + 1 - jN - 2j ) / ( N + 1 )2( N + 2 )
=j( N - j + 1 ) / ( N + 1 )2( N + 2 )

になります。

一様分布に対する順序統計量の確率密度関数

上のグラフは、一様分布における順序統計量の確率密度関数を描画したもので、100 個の標本を抽出したときの最小値・中央値・最大値に対する分布を示しています。左側にある赤い線が最小値 (j = 1)、中央の緑色の分布が中央値 (j = 50)、右側の青い線が最大値 (j = 100) を表します。最小値は限りなくゼロに近い値、最大値は限りなく 1 に近い値になる傾向になり、中央値は t = 0.5 付近の確率が最も高くなっています。これは、直感的に考えた結果とも一致しています。

X(1) と X(N) の同時確率密度は

p( x, y )={ N! / ( 1 - 1 )!( N - 2 )!( N - N )! }{ P(x) }1-1{ ∫{x→y} p(t) dt }N-1-1{ 1 - P(y) }N-Np(x)p(y)
=N( N - 1 ){ ∫{x→y} p(t) dt }N-2p(x)p(y)

で表され、ここでも 0 ≤ x, y ≤ θ 以外の範囲ではゼロになります。また、x が X(1) を、y が X(N) をそれぞれ表すとすれば必ず x < y となることから、0 ≤ x < y ≤ θ の範囲のみを考えて

p( x, y )=N( N - 1 ){ ∫{x→y} 1 / θ dt }N-2( 1 / θ2 )
=N( N - 1 ){ ( y - x ) / θ }N-2( 1 / θ2 )
={ N( N - 1 ) / θN }( y - x )N-2 ( 0 ≤ x < y ≤ θ )

となります。u = x, v = y - x とすれば、x = u, y = u + v となって、ヤコビアンは det( J( u, v ) ) = ∂x/∂u・∂y/∂v - ∂x/∂v・∂y/∂u = 1 となるので、変数変換した結果を q( u, v ) とすれば

q( u, v ) = { N( N - 1 ) / θN }vN-2 ( 0 ≤ v ≤ θ )

となって、v に対する周辺分布 qv(v) は、u が 0 から θ - v までの値を取り得ることから

qv(v) = { N( N - 1 ) / θN }vN-2∫{0→θ-v}du = { N( N - 1 ) / θN }vN-2( θ - v )

と求めることができます。これは、X(N) - X(1) すなわち「範囲(Range)」の確率密度関数を表しています。下図は、標本数 100 における範囲の確率密度関数を示したもので、かなり 1 に近い値を取り得ることがわかります。

一様分布に対する範囲の確率密度関数

また、u = x, v = ( x + y ) / 2 とすれば、x = u, y = 2v - u となって、ヤコビアンは det( J( u, v ) ) = ∂x/∂u・∂y/∂v - ∂x/∂v・∂y/∂u = 2、 y - x = ( 2 v - u ) - u = 2 ( v - u ) となるので、変数変換した結果を r( u, v ) とすれば

r( u, v )={ N( N - 1 ) / θN }{ 2( v - u ) }N-2・2
={ 2N-1N( N - 1 ) / θN }( v - u )N-2

0 ≤ x ≤ θ, 0 ≤ y ≤ θ より、0 ≤ u ≤ θ, 2v - θ ≤ u ≤ 2v であり、x < y という制限から u < v でもあるので、u が正値を持つ範囲は

0 ≤ v < θ / 2 のとき 0 ≤ u < v

θ / 2 ≤ v ≤ θ のとき 2v - θ ≤ u < v

と場合分けされます。それぞれの v に対する周辺分布 rv(v) は

rv(v)={ 2N-1N( N - 1 ) / θN }∫{0→v} ( v - u )N-2 du
={ 2N-1N( N - 1 ) / θN }( -1 / N - 1 )[( v - u )N-1]{0→v}
=( 2v )N-1N / θN( 0 ≤ v < θ / 2 )
={ 2N-1N( N - 1 ) / θN }∫{2v-θ→v} ( v - u )N-2 du
={ 2N-1N( N - 1 ) / θN }( -1 / N - 1 )[( v - u )N-1]{2v-θ→v}
={ 2( θ - v ) }N-1N / θN( θ / 2 ≤ v < θ )

と求められます。これは、ミッドレンジの確率密度関数を表します。

一様分布に対するミッドレンジの確率密度関数

上の図は、θ = 1, N = 100 の場合のミッドレンジの確率密度関数を示したものです。0.5 で確率密度は最大となり、そこから離れるに従って急激に小さくなって、0.46 から 0.54 の範囲より外側はほとんどゼロになります。

任意の確率密度関数に対して順序統計量の確率密度を求めるためのサンプル・プログラムを以下に示します。

/*
  orderedPDF : j 番目の順序統計量 x(j) における確率密度を返す

  double x : 確率変数
  unsigned int n : 標本数
  unsigned int j : 順序統計量の番号
  ContDist& dist : 確率変数の確率密度関数

  戻り値 : 確率密度
*/
double orderedPDF( double x, unsigned int n, unsigned int j, ContDist& dist )
{
  if ( j == 0 || n < j ) return( NAN );
  if ( &dist == 0 ) return( NAN );

  double k = fact( n ) / ( fact( j - 1 ) * fact( n - j ) ); // N! / ( j - 1 )!( N - j )!
  double lowerP = dist.lower_p( x );                        // distの下側確率

  return( k * pow( lowerP, j - 1 ) * pow( 1.0 - lowerP, n - j ) * dist[x] );
}

サンプル・プログラムを使って、標準正規分布 N(0,1) と指数分布 P1(x) = e-x に対し、100 個の標本を抽出した場合の X(j) を j = 1, 50, 100 の三通りで求めた結果を以下に示します。

正規分布指数分布
正規分布に対する順序統計量の確率密度関数 指数分布に対する順序統計量の確率密度関数

左右対称な分布である正規分布では、中央値が平均値 0 で最大の確率密度となるので、標本中央値と標本平均値はほぼ一致すると考えられることを意味します。それに対し、指数分布の場合は対象な分布ではないので、中央値は平均値 1 よりも低い側へ分布しています。指数分布の場合、最小値はほとんどゼロに近い値になる傾向があり、逆に最大値はかなり広い領域に分布しているため、標本を抽出する操作を繰り返して最大値を求めると、そのバラツキは最小値の場合と比べてかなり大きくなります。


2) U-分布

確率密度関数 p(x), q(y) に従う二つの分布からそれぞれ m 個の x と n 個の y を抽出して、それらを昇順に並べることを考えます。もし、q(y) の分布が p(x) よりも値の大きい方へシフトしていれば、並べた結果の列において、前半は x が多く見られ、後半に進むに従って y が増えていくことが予想できます。もし、両方の分布が平均も分散もほぼ等しい状態であるならば、両方のデータはほぼ均等に存在することになるでしょう。この度合いを示す指標として、「H.B.マン(Henry Berthold Mann)」と「D.R.ホイットニー(Donald Ransom Whitney)」が U という値を 1947 年に提唱しました。

U を、y が x よりも前側になる場合の数を表す値とします。例えば、"yxx" の場合は二つの x の前に y があるので、U = 2 になります。m 個の x と n 個の y からなるデータ列に対して、U に対する並べ方の数を Nm,n(U) で表します。例えば、( m, n ) = ( 2, 2 ) の場合、考えられるデータ列の並び方は

xxyy, xyxy, xyyx, yxxy, yxyx, yyxx

の 6 通りで、それぞれの U の値を求めると

データ列U
xxyy0
xyxy1
xyyx2
yxxy2
yxyx3
yyxx4

になります。それぞれの y に対してその後ろ側にある x の個数を数えた和が U になることに注意してください。よって、

N2,2(0) = N2,2(1) = N2,2(3) = N2,2(4) = 1
N2,2(2) = 2

と求めることができます。以下に、他の例を示しておきます。

mn数列U
10x0
01y0
11xy0
yx1
21xxy0
xyx1
yxx2
12xyy0
yxy1
yyx2
22xxyy0
xyxy1
yxxy2
xyyx2
yxyx3
yyxx4

上表において ( m, n ) = ( 2, 2 ) の場合と ( m, n ) = ( 2, 1 ), ( m, n ) = ( 1, 2 ) の場合を見比べると、( m, n ) = ( 2, 2 ) での上半分 xxyy, xyxy, yxxy は ( m, n ) = ( 2, 1 ) のデータ列 xxy, xyx, yxx の末尾に y を付加したもの、上半分の xyyx, yxyx, yyxx は ( m, n ) = ( 1, 2 ) のデータ列 xyy, yxy, yyx の末尾に x を付加したものになっています。末尾に y を付加しても U の値に変化はありませんが、末尾に x を付加した場合、元のデータ列に含まれる y の数の分だけ U の値は増加します。この例でも、上半分での U の値は変化していないのに対し、下半分では U の値が n = 2 だけ増加しています。このことを考慮すると、以下の漸化式が成り立つことがわかります。

Nm,n(U) = Nm-1,n( U - n ) + Nm,n-1(U)

末尾が x であるときは、それを除外した ( m - 1, n ) 個のデータ列において U = U - n であるものが、x を末尾に付加することで U となり、末尾が y であるときは、それを除外した ( m, n - 1 ) 個のデータ列において U であるものがそのまま適用されるので、上記漸化式が成り立つことになります。

x を m 個、y を n 個並べるときの並べ方は ( m + n )! / m!・n! 通りあります。全ての並べ方が等しい確率で発生すると仮定したとき、その確率を Pm,n(U) とすれば、

Pm,n(U) = { m!・n! / ( m + n )! }・Nm,n(U)

と表すことができます。上式に先ほど求めた漸化式を代入すると、

Pm,n(U)={ m!・n! / ( m + n )! }・{ Nm-1,n( U - n ) + Nm,n-1(U) }
={ m / ( m + n ) }{ ( m - 1 )!・n! / ( m + n - 1 )! }・Nm-1,n( U - n ) + { n / ( m + n ) }{ m!・( n - 1 )! / ( m + n - 1 )! }・Nm,n-1(U)
={ m / ( m + n ) }・Pm-1,n( U - n ) + { n / ( m + n ) }・Pm,n-1(U)

という漸化式が求められ、これを利用して U の確率分布を求めていくことができます。但し、U < 0 ならば Pm,n(U) = 0 とし、Pm,0(U) と P0,n(U) は U = 0 の時 1 でその他は 0 になります。また、U は全ての y がデータ列の前側に集まったときに最大値となるので、n 個の y それぞれについて m 個の x が後ろ側になることから U の最大値は mn と求められます。


U の確率分布を求めるためのサンプル・プログラムを以下に示します。

// UMap用のキー
struct UKey {
  unsigned int m;
  unsigned int n;
  // コンストラクタ
  UKey( unsigned int _m, unsigned int _n )
  : m( _m ), n( _n ) {}
  // 比較用演算子の多重定義
  bool operator<( const UKey& dest ) const
  {
    if ( m == dest.m )
      return( n < dest.n );
    else
      return( m < dest.m );
  }
};

// U値を保持する配列用Mapの型
typedef map< UKey, vector<double> > UMap;

/*
  uProb : U の確率分布 Pm,n(U) を求める

  int m, n : データ列の要素数
  int u : U値

  戻り値 : Pm,n(U)
*/
double uProb( int m, int n, int u )
{
  if ( m < 0 || n < 0 || u < 0 || u > m * n )
    return( 0 );

  // Pm,0(U) = P0,n(U) = 1 ( u = 0 ) or 0 ( u ≠ 0 )
  if ( m == 0 || n == 0 )
    return( ( u == 0 ) ? 1 : 0 );

  return(
         (double)m * uProb( m - 1, n, u - n ) / (double)( m + n ) +
         (double)n * uProb( m, n - 1, u ) / (double)( m + n )
         );
}

/*
  uGet : UMapから U 値に対する確率を取得する

  const UMap& uMap : U値に対する確率密度を保持したUMap
  int m, n : 確率密度を取得する対象のキー(データ列の大きさ)
  int u : U値

  戻り値 : Pm,n(U) (未登録の場合は負数を返す)
*/
double uGet( const UMap& uMap, unsigned int m, unsigned int n, unsigned int u )
{
  UKey key( m, n ); // UMap内の対象キー
  UMap::const_iterator cit = uMap.find( key ); // 対象キーに対するuMapの反復子
  if ( cit == uMap.end() )
    return( -1 );

  const vector<double>& uVec = cit->second; // 値を取得する対象配列
  if ( u >= uVec.size() )
    return( 0 );

  return( uVec[u] );
}

/*
  uSet : UMapに U 値に対する確率をセットする

  UMap& uMap : U値に対する確率密度を保持するUMap
  int m, n : 確率密度をセットする対象のキー(データ列の大きさ)
  int u : セットする対象の U 値
  double d : セットする値(確率)

  戻り値 : True ... セット成功 ; False ... u が不正
*/
bool uSet( UMap& uMap, unsigned int m, unsigned int n, unsigned int u, double d )
{
  UKey key( m, n ); // UMap内の対象キー

  // 対象キーが見つからない場合は新規追加
  if ( uMap.find( key ) == uMap.end() ) {
    // あらかじめ負数で初期化しておく
    uMap.insert( pair< UKey, vector<double> >( key, vector<double>( m * n + 1, -1 ) ) );
  }

  vector<double>& uVec = uMap[key]; // 値をセットする対象配列
  if ( u >= uVec.size() )
    return( false );
  uVec[u] = d;

  return( true );
}

/*
  uProb : U の確率分布 Pm,n(U) を求める

  一度再帰的に求めた確率は、UMapにセットして再利用できるようにする
  UMap 内の値はあらかじめ負数で初期化しておく必要がある

  int m, n : データ列の要素数
  int u : U値
  UMap& uMap : 値を保持するUMap

  戻り値 : Pm,n(U)
*/
double uProb( int m, int n, int u, UMap& uMap )
{
  if ( m < 0 || n < 0 || u < 0 )
    return( 0 );

  double d = uGet( uMap, m, n, u ); // 最初に UMap から確率を取得

  // 負数ならばまた求められていないので計算を行う
  if ( d < 0 ) {
    d = uProb( m, n, u );
    uSet( uMap, m, n, u, d ); // 求めた値を UMap にセットしておく
  }

  return( d );
}

/*
  uProb : U の確率分布 Pm,n(U) を求めて配列の形で返す

  一度再帰的に求めた確率は、UMapにセットして再利用できるようにする

  int m, n : データ列の要素数
  vector<double>& uVec : 結果を保持する配列

  戻り値 : true ... 正常終了 ; false ... 引数の異常
*/
bool uProb( int m, int n, vector<double>& uVec )
{
  if ( m < 0 || n < 0 )
    return( false );

  // 確率分布を保持する UMap
  UMap uMap;

  // 確率を計算して代入
  uVec.resize( m * n + 1 );
  for ( unsigned int u = 0 ; u < uVec.size() ; ++u )
    uVec[u] = uProb( m, n, u, uMap );

  return( true );
}

/*
  uProb : 0〜mMax, 0〜nMax までの U の確率分布 Pm,n(U) を求めて UMap の形で返す

  int mMax, nMax : データ列の要素数最大値
  UMap& uMap : 結果を保持するMap

  戻り値 : true ... 正常終了 ; false ... 引数の異常
*/
bool uProb( int mMax, int nMax, UMap& uMap )
{
  if ( mMax < 0 || nMax < 0 )
    return( false );

  // 確率分布を保持する UMap の初期化
  uMap.clear();

  // 確率を計算して代入
  for ( unsigned int n = 0 ; n <= (unsigned int)nMax ; ++n ) {
    for ( unsigned int m = 0 ; m <= (unsigned int)mMax ; ++m ) {
      for ( unsigned int u = 0 ; u <= m * n ; ++u )
        uProb( m, n, u, uMap );
    }
  }

  return( true );
}

uProb という名の関数は 4 つありますが、最も上側の関数が、指定した m, n および U から確率密度を求めるための関数で、二番目の関数は上側の関数を利用して確率密度を求めた上で、その値を再利用できるように UMap ( 実体は map< UKey, vector<double> > ) 型の配列に登録する処理を加えたものです。確率密度は再帰的に求められるので、再帰処理の回数を減らすための対処法になります。三番目の関数は、指定した m, n に対して確率密度を求めるためのもの、最後の関数は、指定した mMax, nMax より小さい場合に対して全ての確率密度を求めるために作成したものです。UMap は連想配列であり、そのキーは UKey という構造体で表現されています。UKey は m, n を保持して、この二つの値の大小関係によってキーが判別されます。UMap の値は配列(vector)であり、この中に各 U 値に対する確率密度を計算します。


サンプル・プログラムを利用して、m = 3 の場合と m = 8 の場合について、n を 0 から 8 まで変化させたときの、各 U 値に対する確率密度を求めた結果を以下に示します。

m = 3

un
012345678
01.000e+002.500e-011.000e-015.000e-022.857e-021.786e-021.190e-028.333e-036.061e-03
12.500e-011.000e-015.000e-022.857e-021.786e-021.190e-028.333e-036.061e-03
22.500e-012.000e-011.000e-015.714e-023.571e-022.381e-021.667e-021.212e-02
32.500e-012.000e-011.500e-018.571e-025.357e-023.571e-022.500e-021.818e-02
42.000e-011.500e-011.143e-017.143e-024.762e-023.333e-022.424e-02
51.000e-011.500e-011.143e-018.929e-025.952e-024.167e-023.030e-02
61.000e-011.500e-011.429e-011.071e-018.333e-025.833e-024.242e-02
71.000e-011.143e-011.071e-018.333e-026.667e-024.848e-02
85.000e-021.143e-011.071e-019.524e-027.500e-026.061e-02
95.000e-028.571e-021.071e-019.524e-028.333e-026.667e-02
105.714e-028.929e-029.524e-028.333e-027.273e-02
112.857e-027.143e-028.333e-028.333e-027.273e-02
122.857e-025.357e-028.333e-028.333e-027.879e-02
133.571e-025.952e-027.500e-027.273e-02
141.786e-024.762e-026.667e-027.273e-02
151.786e-023.571e-025.833e-026.667e-02
162.381e-024.167e-026.061e-02
171.190e-023.333e-024.848e-02
181.190e-022.500e-024.242e-02
191.667e-023.030e-02
208.333e-032.424e-02
218.333e-031.818e-02
221.212e-02
236.061e-03
246.061e-03

m = 8

un
012345678
01.000e+001.111e-012.222e-026.061e-032.020e-037.770e-043.330e-041.554e-047.770e-05
11.111e-012.222e-026.061e-032.020e-037.770e-043.330e-041.554e-047.770e-05
21.111e-014.444e-021.212e-024.040e-031.554e-036.660e-043.108e-041.554e-04
31.111e-014.444e-021.818e-026.061e-032.331e-039.990e-044.662e-042.331e-04
41.111e-016.667e-022.424e-021.010e-023.885e-031.665e-037.770e-043.885e-04
51.111e-016.667e-023.030e-021.212e-025.439e-032.331e-031.088e-035.439e-04
61.111e-018.889e-024.242e-021.818e-027.770e-033.663e-031.709e-038.547e-04
71.111e-018.889e-024.848e-022.222e-021.010e-024.662e-032.331e-031.166e-03
81.111e-011.111e-016.061e-023.030e-021.399e-026.660e-033.263e-031.709e-03
98.889e-026.667e-023.434e-021.709e-028.325e-034.196e-032.176e-03
108.889e-027.273e-024.242e-022.176e-021.099e-025.594e-032.953e-03
116.667e-027.273e-024.646e-022.564e-021.332e-026.993e-033.730e-03
126.667e-027.879e-025.455e-023.108e-021.698e-029.013e-034.895e-03
134.444e-027.273e-025.657e-023.497e-021.965e-021.088e-025.983e-03
144.444e-027.273e-026.263e-024.040e-022.364e-021.336e-027.537e-03
152.222e-026.667e-026.263e-024.429e-022.697e-021.570e-029.013e-03
162.222e-026.061e-026.667e-024.895e-023.130e-021.865e-021.096e-02
174.848e-026.263e-025.128e-023.430e-022.129e-021.274e-02
184.242e-026.263e-025.439e-023.863e-022.455e-021.507e-02
193.030e-025.657e-025.517e-024.096e-022.735e-021.717e-02
202.424e-025.455e-025.672e-024.462e-023.061e-021.981e-02
211.818e-024.646e-025.517e-024.629e-023.326e-022.207e-02
221.212e-024.242e-025.439e-024.862e-023.621e-022.479e-02
236.061e-033.434e-025.128e-024.895e-023.838e-022.704e-02
246.061e-033.030e-024.895e-025.028e-024.087e-022.976e-02
252.222e-024.429e-024.895e-024.227e-023.178e-02
261.818e-024.040e-024.862e-024.382e-023.419e-02
271.212e-023.497e-024.629e-024.429e-023.582e-02
281.010e-023.108e-024.462e-024.491e-023.776e-02
296.061e-032.564e-024.096e-024.429e-023.877e-02
304.040e-032.176e-023.863e-024.382e-024.002e-02
312.020e-031.709e-023.430e-024.227e-024.033e-02
322.020e-031.399e-023.130e-024.087e-024.087e-02
331.010e-022.697e-023.838e-024.033e-02
347.770e-032.364e-023.621e-024.002e-02
355.439e-031.965e-023.326e-023.877e-02
363.885e-031.698e-023.061e-023.776e-02
372.331e-031.332e-022.735e-023.582e-02
381.554e-031.099e-022.455e-023.419e-02
397.770e-048.325e-032.129e-023.178e-02
407.770e-046.660e-031.865e-022.976e-02
414.662e-031.570e-022.704e-02
423.663e-031.336e-022.479e-02
432.331e-031.088e-022.207e-02
441.665e-039.013e-031.981e-02
459.990e-046.993e-031.717e-02
466.660e-045.594e-031.507e-02
473.330e-044.196e-031.274e-02
483.330e-043.263e-031.096e-02
492.331e-039.013e-03
501.709e-037.537e-03
511.088e-035.983e-03
527.770e-044.895e-03
534.662e-043.730e-03
543.108e-042.953e-03
551.554e-042.176e-03
561.554e-041.709e-03
571.166e-03
588.547e-04
595.439e-04
603.885e-04
612.331e-04
621.554e-04
637.770e-05
647.770e-05

結果を見てまず気付くのが、m, n を固定したとき、U に対する分布は左右対称となって、中央が最も確率密度が高くなるということです。また、( m, n ) = ( 3, 8 ) と ( m, n ) = ( 8, 3 ) の分布を比較すると、どちらも同じ分布であることがわかります。これらの性質は、次のように考えることで理解することができます。U' を、今度は x が y よりも前側になる場合の数とします。( m, n ) = ( 2, 2 ) の場合で U と U' を比較すると、次のようになります。

データ列UU'
xxyy04
xyxy13
xyyx22
yxxy22
yxyx31
yyxx40

このように、U と U' の和はどのデータ列に対しても mn と等しくなります。データ列の中で、x, y それぞれの要素を選ぶ場合の数は mn になります。それぞれに対して y が x より前側にあれば U が、逆に x が y より前側にあれば U' が加算されることになるので、両者の和は mn になります。従って、U' = mn - U であり、U に対するデータ列の個数と U' = mn - U に対するデータ列の個数は等しくなることから、Pm,n(mn - U) = Pm,n(U) が成り立つことになります。また、U' の分布は Pn,m(U') であることを意味していることから、Pn,m(U) = Pm,n(U) であることもわかります。この結果を考慮すると、U 値に対する確率分布を計算するサンプル・プログラムに以下の処理を加えることで高速化させることができます。

// 左右対称のため u の値は mn / 2 以下にできる
if ( u > m * n / 2 )
  u = m * n - u;

分布が左右対称で中央部分が最も大きく、確率変数 U は 0 から mn までの値をとることから、Pm,n(U) の期待値を Em,n[U] とすると、

Em,n[U] = mn / 2

になります。これは、mn 個ある x, y の選び方それぞれに対して、y が x の前側にある確率が 1 / 2 であることからも導くことができます。

υ2 = ( U - mn / 2 )2 として、分散 Em,n2] を求めてみます。

Em,n2] = ΣU{0→mn}( ( U - mn / 2 )2Pm,n(U) )

より、Pm,n(U) の漸化式を代入すると

Em,n2]=ΣU{0→mn}( ( U - mn / 2 )2[ { m / ( m + n ) }・Pm-1,n( U - n ) + { n / ( m + n ) }・Pm,n-1(U) ] )
={ m / ( m + n ) }ΣU{0→mn}( ( U - mn / 2 )2・Pm-1,n( U - n ) ) + { n / ( m + n ) }ΣU{0→mn}( ( U - mn / 2 )2・Pm,n-1(U) )

前半の項の和を計算すると、Pm-1,n( U - n ) は n ≤ U ≤ ( m - 1 )n + n = mn で正値となるので 0 ≤ U ≤ n - 1 までの範囲は無視できて、U' = U - n とすれば 0 ≤ U' ≤ ( m - 1 )n となって

{ m / ( m + n ) }ΣU{0→mn}( ( U - mn / 2 )2・Pm-1,n( U - n ) )
={ m / ( m + n ) }ΣU'{0→(m-1)n}( ( U' + n - mn / 2 )2・Pm-1,n( U' ) )
={ m / ( m + n ) }ΣU'{0→(m-1)n}( { U' + ( n / 2 ) - ( m - 1 )n / 2 }2・Pm-1,n( U' ) )
={ m / ( m + n ) }ΣU'{0→(m-1)n}( [ { U' - ( m - 1 )n / 2 }2 + n{ U' - ( m - 1 )n / 2 } + n2 / 4 ]・Pm-1,n( U' ) )
={ m / ( m + n ) }ΣU'{0→(m-1)n}( [ { U' - ( m - 1 )n / 2 }2 + nU' - ( 2m - 3 )n2 / 4 ]・Pm-1,n( U' ) )
={ m / ( m + n ) }{ Em-1,n2] + nEm-1,n[U] - ( 2m - 3 )n2 / 4 }
={ m / ( m + n ) }{ Em-1,n2] + ( m - 1 )n2 / 2 - ( 2m - 3 )n2 / 4 }
={ m / ( m + n ) }( Em-1,n2] + n2 / 4 )

後半の項の和も同様に計算すると、Pm,n-1(U) は 0 ≤ U ≤ m( n - 1 ) で正値となるので m( n - 1 ) + 1 ≤ U ≤ mn までの範囲は無視できて、

{ n / ( m + n ) }ΣU{0→mn}( ( U - mn / 2 )2・Pm,n-1(U) )
={ n / ( m + n ) }ΣU{0→m(n-1)}( ( U - mn / 2 )2・Pm,n-1(U) )
={ n / ( m + n ) }ΣU{0→m(n-1)}( { U - ( m / 2 ) - m( n - 1 ) / 2 }2・Pm,n-1(U) )
={ n / ( m + n ) }ΣU{0→m(n-1)}( [ { U - m( n - 1 ) / 2 }2 - m{ U - m( n - 1 ) / 2 } + m2 / 4 ]・Pm,n-1(U) )
={ n / ( m + n ) }ΣU{0→m(n-1)}( [ { U - m( n - 1 ) / 2 }2 - mU + m2( 2n - 1 ) / 4 ]・Pm,n-1(U) )
={ n / ( m + n ) }{ Em,n-12] - mEm,n-1[U] + m2( 2n - 1 ) / 4 }
={ n / ( m + n ) }{ Em,n-12] - m2( n - 1 ) / 2 + m2( 2n - 1 ) / 4 }
={ n / ( m + n ) }( Em,n-12] + m2 / 4 )

よって、

Em,n2]={ m / ( m + n ) }( Em-1,n2] + n2 / 4 ) + { n / ( m + n ) }( Em,n-12] + m2 / 4 )
={ m / ( m + n ) }Em-1,n2] + { n / ( m + n ) }Em,n-12] + mn2 / 4( m + n ) + m2n / 4( m + n )
={ m / ( m + n ) }Em-1,n2] + { n / ( m + n ) }Em,n-12] + mn / 4

となって、これが分散の漸化式となります。m = 0 または n = 0 の場合、直接計算することで E0,n2] = Em,02] = 0 と求められるので、これを利用して E1,12] を求めると

E1,12] = E0,12] / 2 + E1,02] / 2 + 1 / 4 = 1 / 4

また、E2,12], E1,22] は

E2,12] = 2E1,12] / 3 + E2,02] / 3 + 2 / 4 = 2 / 3

E1,22] = E0,22] / 3 + 2E1,12] / 3 + 2 / 4 = 2 / 3

よって、E2,22] は

E2,22] = 2E1,22] / 4 + 2E2,12] / 4 + 4 / 4 = 5 / 3

これだけでは規則性がよくわかりませんが、m を 1 に固定した上で n の値を大きくしていくと Em,n2] は次のようになります。

n012345678
Em,n2]0 / 123 / 128 / 1215 / 1224 / 1235 / 1248 / 1263 / 1280 / 12
Δ-3 / 125 / 127 / 129 / 1211 / 1213 / 1215 / 1217 / 12

値は全て分母を 12 に揃えています。また、下側の値 Δ は、前の値との差を表しています。この結果を見ると、分子の値の増分は 2 ずつ大きくなっていることがわかります。このような場合、漸化式の解として例えば n2 が含まれているとすれば、( n + 1 )2 - n2 = 2n + 1 なので、増分は n が 1 増えるに従い 2 ずつ大きくなります。漸化式の対称性から、n を 1 に固定して m の値を変化させた場合も全く同じ結果になるので、m2 も含まれていると考えられます。少なくとも、m, n それぞれ次数が 2 の項は存在することになります。

もうひとつ、m, n それぞれを 1 ずつ増加した場合の Em,n2] の値を計算すると、次のようになります。

m,n012345678
Em,n2]0 / 123 / 1220 / 1263 / 12144 / 12275 / 12468 / 12735 / 121088 / 12
Δ-3 / 1217 / 1243 / 1281 / 12131 / 12193 / 12267 / 12353 / 12
ΔΔ--14 / 1226 / 1238 / 1250 / 1262 / 1274 / 1286 / 12

ΔΔ は、前の値からの増分 Δ に対してさらに増分を計算した結果です。ΔΔ が 1 ずつ増加していることから、Em,n2] は m = n のとき 3 次項であることになります。なぜなら、n3 が含まれているとすれば、( n + 1 )3 - n3 = 3n2 + 3n + 1 で、{ 3( n + 1 )2 + 3( n + 1 ) + 1 } - ( 3n2 + 3n + 1 ) = 6( n + 1 ) なので、n が 1 増えれば ΔΔ は定数 6 だけ増えることになるからです。

E0,n2] = Em,02] = 0 なので、漸化式の解は mn を含むと考えれば、m = n のとき二次式となります。E1,12] = 3 / 12 となるような、m と n を含む項として ( m + n + 1 ) /12 を加えて、mn( m + n + 1 ) / 12 とすれば、m = n のとき三次式となり、m, n の一方を固定したときの次数は 2 となります。実際、この式が Em,n2] を求める一般解となります。


U-分布を求めるには漸化式を使って再帰的に処理する必要があり、これは非常に時間がかかります。しかし、m, n が大きい場合、U 値に対する分布は正規分布に近づくため、これを利用して近似的に確率密度を求めることができます。下図は ( m, n ) = ( 12, 12 ) の場合の U 値の分布と、平均 mn / 2 = 72, 分散 mn( m + n + 1 ) / 12 = 300 の正規分布 N( 72, 300 ) のグラフを重ねあわせたもので、棒で示されたほうが U 値の分布、曲線が正規分布をそれぞれ表しています。これを見ると、両者が非常に近い分布形状を持っていることがわかります。

U分布(m=n=12)

U 値は、データ列の中から mn 通りある x と y のペアを抽出して、順序が xy なら 0、yx なら 1 としたときの和を表しているのでした。xy になるか、yx になるかは 1 / 2 の確率であり、それを何回も繰り返した時の和が U 値であるとすれば、中心極限定理から、その分布が正規分布に従うことは直感的に理解できると思います。また、その平均が mn / 2、分散が mn( m + n + 1 ) / 12 になることから、正規分布は N( mn / 2, mn( m + n + 1 ) / 12 ) で表されることになります。


3) マン・ホイットニーのU検定(Mann-Whitney U Test)

前述の通り、二つの分布の平均と分散が等しければ、U 値の確率密度は中心部が最も高く、0 や mn に近づくに従って小さくなっていきます。分布が完全に離れていれば、U 値は逆に 0 や mn に非常に近い値になるので、この性質を利用することで二つの分布が等しいものかを判断することができます。このような検定法を「マン・ホイットニーの U 検定(Mann-Whitney U Test)」といいます。ここで、前回紹介した例題を U 検定を使って検定してみたいと思います。

装置にある改良を行った後の製品の歩留まりを調べたところ、改良前と後で次のような結果が得られた。

改良前85.382.686.482.188.978.490.388.881.480.6
改良後79.481.687.284.186.891.190.190.482.683.0

改良前後において、歩留まりに差は生じたか。

改良前を x 群、改良後を y 群としてデータをソートすると、次のような結果になります。

データ列78.479.480.681.481.682.182.682.683.084.185.386.486.887.288.888.990.190.390.491.1
群名xyxxyxxyyyxxyyxxyxyy

y 群の各要素に対して、x 群の中でより大きい要素の数 u を数えます。

データ列78.479.480.681.481.682.182.682.683.084.185.386.486.887.288.888.990.190.390.491.1
群名xyxxyxxyyyxxyyxxyxyy
u-9--7--5.555--33--1-00

ここで、x 群と y 群の値が等しい場合は、0.5 を加算することになります。よって、82.6 という値を持つ要素に対しては、本来整数値であるはずの u 値が 5.5 という半端な数になっています。

上記結果から、U = Σu = 9 + 7 + 5.5 + 5 + 5 + 3 + 3 + 1 + 0 + 0 = 38.5 になります。要素数の積は 10 * 10 = 100 なので中央値よりは小さな値であり、y 群の方が分布は大きい側に偏っていると考えることができます。同じ分布だとした場合に確率としてどの程度発生しうるのかを調べるため、まずは ( m, n ) = ( 10, 10 ) の場合の U 値の確率密度を計算すると、

U0123456789
確率密度0.0000.0000.0000.0000.0000.0000.0000.0000.0000.000
累積確率0.0000.0000.0000.0000.0000.0000.0000.0000.0000.001
U10111213141516171819
確率密度0.0000.0000.0000.0010.0010.0010.0010.0010.0020.002
累積確率0.0010.0010.0010.0020.0030.0030.0040.0060.0070.009
U20212223242526272829
確率密度0.0020.0030.0030.0040.0050.0050.0060.0070.0080.009
累積確率0.0120.0140.0180.0220.0260.0320.0380.0450.0530.062
U30313233343536373839
確率密度0.0100.0110.0120.0140.0150.0160.0180.0190.0200.021
累積確率0.0720.0830.0950.1090.1240.1400.1570.1760.1970.218
U40414243444546474849
確率密度0.0230.0240.0250.0260.0270.0280.0280.0290.0290.029
累積確率0.2410.2640.2890.3150.3420.3700.3980.4270.4560.485

のようになります(U ≥ 50 の範囲は除外しています)。累積確率から、U ≤ 38.5 となる確率はおよそ 20% で、帰無仮説を「分布は等しい」としたとき、危険率 10% (片側検定) としても仮説は棄却できないので、差があったとは認められないことになります(どちらとも判別できないというのが結論になります)。y 群の分布の方が高い側にシフトしているとみなせるのは、危険率 10% の場合 U ≤ 32、危険率 5% なら U ≤ 27 程度になる必要があります。

正規分布で近似した場合、その平均は 10 * 10 / 2 = 50、分散は 10 * 10 * ( 10 + 10 + 1 ) / 12 = 175 になります。U ≤ 38.5 のときの p 値は 0.192 であり、危険率 10% (片側検定) での棄却域は U ≤ 33、危険率 5% なら U ≤ 28 程度となって、実際の分布をよく近似していることがわかります。両分布の標本数が 10 程度以上であれば、正規分布で近似しても精度的にはそれほど問題ないと思います。

ところで、x 群の各要素に対して y 群の中でより大きい要素の数として U を求めてみた場合の結果がどのようになるかを調べてみると、次のようになります。

データ列78.479.480.681.481.682.182.682.683.084.185.386.486.887.288.888.990.190.390.491.1
群名xyxxyxxyyyxxyyxxyxyy
u10-99-87.5---55--33-2--

より、U = Σu = 10 + 9 + 9 + 8 + 7.5 + 5 + 5 + 3 + 3 + 2 = 61.5 なので、先程の結果 38.5 との和は 100 で、両要素数の積と等しくなります。これは、U-分布のところでも考察した通り、mn 個ある x, y のペアのなかで順番が yx の場合に 1 を加算するか、xy の場合に 1 を加算するかによって両方の U 値が計算できるため、両者の和は mn と等しくなることから得られる結果です。分布の対称性から Pm,n( U ) = Pm,n( mn - U ) であり、さらに Pm,n( U ) = Pn,m( U ) でもあることから、結局どちらの方法で計算しても結果は同じことになります。しかし、U 値と mn / 2 の大小関係から、x, y のどちらの分布がより大きい側へ偏っているかが変化することは注意が必要です。

U 値の別の求め方として、各要素の順位から計算する手法があります。

データ列78.479.480.681.481.682.182.682.683.084.185.386.486.887.288.888.990.190.390.491.1
群名xyxxyxxyyyxxyyxxyxyy
u-9--7--5.555--33--1-00
順位1234567.57.591011121314151617181920

上の表は、先ほど求めた u 値 ( y に対して x が後ろ側にある場合の数 ) に順位を付加したものです。2 番目にある y の u 値は 9 となっていますが、この y より後ろ側には x, y 合わせて 20 - 2 = 18 個の要素があり、そのうち y は自分自身を除く 9 個なので、それを引いた数の 9 が u 値と等しくなります。同様な考え方によって u 値を計算すると、

20 - 2   - 9 = 9
20 - 5   - 8 = 7
20 - 7.5 - 7 = 5.5
20 - 9   - 6 = 5
20 - 10  - 5 = 5
20 - 13  - 4 = 3
20 - 14  - 3 = 3
20 - 17  - 2 = 1
20 - 19  - 1 = 0
20 - 20  - 0 = 0

と求められることがわかります。各式の最初の値は m + n を、2 番目は順位を、そして 3 番目は n - 1 で表されるのでその和は n( n - 1 ) / 2 で計算できます。従って、y の各順位の和を Sy とすれば、

U = n( m + n ) - Sy - n( n - 1 ) / 2 = mn + n( n + 1 ) / 2 - Sy

で計算することもできます。順位を求めるためには x, y 両方のデータを並べ替える必要があり、コンピュータを使って処理する場合は前の求め方に比べて早くなるというわけではないですが、手計算で行う場合は順位を求める手法のほうが簡単です。


マン・ホイットニーのU検定を行うためのサンプル・プログラムを以下に示します。

/*
  calcUValue : U 値を求めて返す

  配列 x はソートをするため、元のデータ列が変換されることに注意

  vector<double>& x : データ列 x
  const vector<double>& y : データ列 y

  戻り値 : U 値
*/
double calcUValue( vector<double>& x, const vector<double>& y )
{
  sort( x.begin(), x.end() );

  // y の各要素より大きな x の要素数を求め、その和を U 値として計算
  double u = 0; // U値
  for ( unsigned int i = 0 ; i < y.size() ; ++i ) {
    vector<double>::const_iterator cit = lower_bound( x.begin(), x.end(), y[i] );
    u += x.end() - cit;
    // x,y の要素が等しい場合は 0.5 と換算
    while ( cit != x.end() ) {
      if ( *cit == y[i] ) u -= 0.5;
      ++cit;
    }
  }

  return( u );
}

/*
  UTest : マン・ホイットニーの U-検定

  配列 x はソートをするため、元のデータ列が変換されることに注意

  vector<double>& x : データ列 x
  const vector<double>& y : データ列 y
  double a : 危険率

  戻り値 : true ... 帰無仮説(分布が等しい)が棄却された ; false ... 帰無仮説が保留された
*/
bool UTest( vector<double>& x, const vector<double>& y, double a )
{
  cout << "***** U-Test *****" << endl << endl;

  // 不偏分散
  double ux = sampleVariance( x, true );
  double uy = sampleVariance( y, true );

  // 不偏分散比の確認
  if ( ux != 0 && uy != 0 ) {
    cout << "The ratio of unbiased variance ux^2 / uy^2 = " << ux / uy << endl << endl;
  } else {
    cout << "Failed to get the ratio of unbiased variance ux^2 / uy^2." << endl;
    cout << "Maybe all x data ( or all y data ) are same, or data size is only one." << endl << endl;
  }

  // データ列の要素数
  unsigned int m = x.size();
  unsigned int n = y.size();

  double u = calcUValue( x, y ); // U値
  cout << "U = " << u << endl;

  // U 値から分布の大小を判定
  if ( u >= (double)( m * n ) / 2.0 ) {
    cout << "Distribution X seems to be more than distribution Y." << endl << endl;
    u = (double)( m * n ) - u;
  } else {
    cout << "Distribution Y seems to be more than distribution X." << endl << endl;
  }

  // 下側確率の計算
  double p = 0; // 下側確率
  if ( m <= 10 && n <= 10 ) {
    vector<double> uVec;
    uProb( m, n, uVec );
    unsigned int i = 0;
    for ( ; (double)i < u ; ++i )
      p += uVec[i];
    p += ( u - (double)( i - 1 ) ) * uVec[i];
  } else {
    NormalDistribution nd( m * n / 2.0, sqrt( m * n * ( m + n + 1 ) / 12.0 ) );
    p = nd.lower_p( u );
  }

  // p-値と検定結果の出力
  cout << "p-value = " << p << endl;
  if ( p < a / 2 ) {
    cout << "Null hypothesis is rejected." << endl;
    return( true );
  } else {
    cout << "Null hypothesis is accepted." << endl;
    return( false );
  }
}

関数 calcValue は、与えられたデータ列から U 値を求めるためのもので、y の要素に対して x の要素の方が大きい場合の数から求めるため、x の要素が y よりも大きい傾向にあるほど U 値は大きくなります。ここでは、順位から求めるのではなく yx の要素を直接比較しています。x をソートした上で、y の各要素が x の各値に対して大きくなる(または等しくなる)最初の位置を求め、末尾までにある要素数を求める操作を繰り返し、和を求めています。この時に利用しているのが STL(Standard Template Library) にある lower_bound 関数で、ソート済みのコンテナクラスから指定した値と等しくなる最初の要素の位置を返します。等しい値が見つからない場合、より大きな要素の中で最初の位置を返すので、今回の処理には最適な関数となります。

関数 UTest は、U-検定を行うためのプログラムです。プログラムの中で、U 値を求めて検定を行う前に、両方のデータから不偏分散を求め、その比率を計算しています。二つの標本が同じ分布から抽出されたものであるかを検定する手法を前回「二標本の検定(Two Sample Test)」として紹介しました。分布が等しいかを検定する前に分散が等しいと見なせるかを確認する必要があるため、不偏分散比が F-分布に従うことを利用したのですが、それは標本が正規分布に従うと仮定した場合であり、今回はそのまま利用することはできません。従って、不偏分散比を計算するまでを行って、それが極端に 1 から外れていないかを確認することができるようにだけしてあります。
U-検定の場合も分布の分散の比によって検定結果は影響を受けます。例えば、分散が極端に小さな分布と極端に大きな分布について U-検定を行った場合、それぞれの分布の平均や中央値が外れていても、分散の大きな分布の中に分散の小さな分布が含まれてしまっていれば、同じ分布である( U 値に偏りがない)という結果になってしまいます。


マン・ホイットニーの U 検定は、「ウィルコクソンの順位和検定(Wilcoxon Rank-sum Test)」とも呼ばれています。マン・ホイットニーの U 値に対する論文は 1947 年に発表されましたが、それより 2 年前の 1945 年に「フランク・ウィルコクソン(Frank Wilcoxon)」が同様の手法を発表しています。しかし、ウィルコクソンは二つの標本の数が同じ場合のみを想定しており、U 値の概念を導入し、異なる標本数の場合まで拡張したのがマンとホイットニーです。このことから、「マン・ホイットニー・ウィルコクソン検定(Mann-Whitney-Wilcoxon Test ; MWW Test)」と呼ばれることもあるようです。

前回紹介した平均値の検定( t-検定 )分散の検定は、母集団が正規分布に従うという強い前提条件があるため、どのような母集団に対しても利用できるというわけではありませんでした(その代わり、かなり定量的にデータ解析することができます)。このように、特定の分布に従うことを前提条件に行う統計手法は「パラメトリック手法(Parametric Procedure)」と呼ばれます。それに対し、U-検定のように分布に依存しない手法は「ノンパラメトリック手法(Nonparametric Procedure)」と呼ばれ、次のような特徴があります。

  1. 結果に対して、母集団の分布に対する影響が小さい ... 頑健性(Robustness)
  2. 精度は悪いが計算が簡単 ... 簡便性(Simple and Quick)

コンピュータを利用する場合、簡便性はあまり重要ではないので、ノンパラメトリック手法を利用する主な理由は「頑健性」にあると思います。特に分布に依存しないことから利用できる場合は非常に多くなります。その代わり、精度的にはパラメトリック手法よりも悪くなります。

また U 検定は、実際の値がわからず順序だけしか情報がないような場合にも適用できるという利点もあります。但し、この場合は分散を求めることができず、分布のバラツキを確認することは不可能になります。順序に関する情報のみから U 値を求めるサンプル・プログラムを以下に示します。

/*
  calcUValue : U 値を求めて返す

  x, y をそれぞれ true, false で表現した bool型配列から計算する

  const vector<bool>& data : データ列

  戻り値 : U 値
*/
double calcUValue( const vector<bool>& data )
{
  // x, y 各要素数を先に計算
  unsigned int xCnt = 0;
  unsigned int yCnt = 0;
  for ( unsigned int i = 0 ; i < data.size() ; ++i ) {
    if ( data[i] )
      ++xCnt;
    else
      ++yCnt;
  }

  if ( xCnt == 0 || yCnt == 0 ) return( 0 );

  // y の各要素より大きな x の要素数を求め、その和を計算
  double u = 0; // U 値
  for ( unsigned int i = 0 ; i < data.size() ; ++i ) {
    if ( data[i] ) {
      --xCnt;
    } else {
      u += xCnt;
      if ( --yCnt == 0 ) break;
    }
  }

  return( u );
}

配列の中にある x, y 各要素の数を先に求めておき、あとは y が見つかったら残りの x の要素数を加算する処理を繰り返すだけの簡単な処理です。なお、この場合は要素が等しいかどうかを判断することができないため、U 値は必ず整数値となります。


4) ウィルコクソンの符号順位検定(Wilcoxon Signed-rank Test)

二つの分布に対応がある場合、一標本の検定に置き換えることが可能であることを、前回の中で少しだけ紹介しました。例えば、被験者に対して投薬前後の分布に差があったかを調べる場合などがそれに該当し、各被験者に対して投薬前後での値から差や比率などを求めた上で、その分布状態を調べることになります。差を利用するのであれば、それがゼロであるという仮説から検定を行うことになります。「ウィルコクソンの符号順位検定(Wilcoxon Signed-rank Test)」は、同様の処理を分布に依存しない(ノンパラメトリックな)場合について行うための手法になります。次の例題を使って説明したいと思います。

10人の被験者に対して投薬を実施して、ある年とその一年後の体重を測定した結果を調べたところ、次のような表になった。

投薬前73.369.779.887.578.367.251.983.188.681.5
投薬後70.072.680.386.575.570.753.381.786.781.5

投薬前後において、体重に差は生じたか。

まずは、投薬前のデータを ai, 投薬後のデータを bi としてその差 di = ai - bi を求めます。

投薬前(ai)73.369.779.887.578.367.251.983.188.681.5
投薬後(bi)70.072.680.386.575.570.753.381.786.781.5
差分(di)3.3-2.9-0.51.02.8-3.5-1.41.41.90.0

次に、差分の絶対値 |di| が小さい順に順位付けを行います。

投薬前(ai)73.369.779.887.578.367.251.983.188.681.5
投薬後(bi)70.072.680.386.575.570.753.381.786.781.5
差分(di)3.3-2.9-0.51.02.8-3.5-1.41.41.90.0
順位8712693.53.55
差分の符号+--++--++

順位付けのとき、差分がゼロの要素は無視して、以降の処理においては対象外とします。また、同値の要素があった場合、それぞれを適当に順位付けした上でその順位の平均を求め、すべての要素に適用します。上記の例では差分の絶対値が 1.4 の要素が二つあるので、それぞれに順位付けをした場合 3, 4 になります。よって、( 3 + 4 ) / 2 = 3.5 が求める順位となります。

最後に、差分が正値・負値であるものに分けて順位和を計算します。

正値 : 8 + 2 + 6 + 3.5 + 5 = 24.5
負値 : 7 + 1 + 9 + 3.5     = 20.5

両方の順位和を合計すると、差分がゼロのものを除いたデータ数 N に対する 1 から N までの和 Σi{1→N}( i ) = N( N + 1 ) / 2 と等しくなります。よって、どちらかが小さな値になるほど、片側は大きくなる傾向にあります。両分布に対して差異がなければ、正値および負値になる度合いは半々であると考えられ、両者の差の期待値はゼロに等しくなります。しかし、分布に偏りがあれば差は大きくなり、その絶対値の最大は N( N + 1 ) / 2 に近づきます。この例では、正値の方が多少大きいため投薬前の方が分布が大きい方へシフトしているものの、両者の差は期待値に非常に近いため分布に差異はないと判断することができます。

順位和の差分を Δs を、( 正値の順位和 - 負値の順位和 ) で定義したとき、分布の要素数が N のときの、Δs となる場合の数を CN(Δs) で表します。差分が正値となる要素を [+], 負値となる要素を [-] で表し、差分の絶対値から昇順で並べる場合を考えると、 正値の個数を m としたとき、N = 1 では m = 0 での並べ方 [-] と m = 1 での並べ方 [+] しか存在しないので、

C1(-1) = 1
C1(1) = 1

になります。N = 2 では、m = 0 のとき [-][-], m = 1 では [+][-] と [-][+], m = 2 は [+][+] で、Δs は順番に -3, -1, 1, 3 となるので、

C2(-3) = 1
C2(-1) = 1
C2(1) = 1
C2(3) = 1

という結果が得られます。N = 1 の場合の [-] という結果に対し、最も差分の大きい要素として末尾に [+] を追加すると、Δs は -1 から 2 増加して 1 へ、また [-] を追加すると、-1 から 2 減少して -3 へ、それぞれ変化することに注意してください。N = 3 の場合も同様に計算すると、

C3(-6) = 1 ([-][-][-])
C3(-4) = 1 ([+][-][-])
C3(-2) = 1 ([-][+][-])
C3(0) = 2 ([-][-][+],[+][+][-])
C3(2) = 1 ([+][-][+])
C3(4) = 1 ([-][+][+])
C3(6) = 1 ([+][+][+])

と求めることができます。ここで、C3(0) は N = 2 のとき Δs = -3 となるデータ列 [-][-] の末尾に [+] を追加して 3 増加した結果と、Δs = 3 となるデータ列 [+][+] の末尾に [-] を追加して 3 減少した結果が含まれています。これらのことから、U 分布で検討した場合と同様に、漸化式として次の結果が得られます。

CN(Δs) = CN-1(Δs-N) + CN-1(Δs+N)

この漸化式を利用して、例えば C4(2) を求めると、

C4(2) = C3(-2) + C3(6) = 2

になります。但し、Δs は -N( N + 1 ) / 2 から N( N + 1 ) / 2 までの範囲のみ値を持ち、それ以外はゼロになります。また、範囲内でも必ず値を持つとは限らないことも上記の結果から明らかです。

[-] と [+] を並べる場合の数は 2N です。分布に差がなく、これらの場合がすべて同等の確率で起こるとしたとき、差分が Δs になる確率 PN(Δs) は

PN(Δs)=CN(Δs) / 2N
={ CN-1(Δs-N) + CN-1(Δs+N) } / 2N
={ 2N-1PN-1(Δs-N) + 2N-1PN-1(Δs+N) } / 2N
={ PN-1(Δs-N) + PN-1(Δs+N) } / 2

となって、これで確率密度を求めることができます。正値と負値を入れ替えた場合、Δs は符号が逆転するだけなので、PN( Δs ) = PN( -Δs )、つまり分布は左右対称であり、期待値はゼロになります。分散を VN とすると、

VN=ΣΔs{-N(N+1)/2→N(N+1)/2}( ( Δs - 0 )2PN(Δs) )
=ΣΔs{-N(N+1)/2→N(N+1)/2}( Δs2{ PN-1(Δs-N) + PN-1(Δs+N) } / 2 )

前半の項は、PN-1(Δs-N) が | Δs - N | ≤ ( N - 1 ){ ( N - 1 ) + 1 } / 2 = N( N - 1 ) / 2 で正値となるのでそれ以外の範囲は無視できて、Δs - N = k とすれば

ΣΔs{-N(N+1)/2→N(N+1)/2}( Δs2PN-1(Δs-N) / 2 )
=Σk{-N(N-1)/2→N(N-1)/2}( ( k + N )2PN-1(k) / 2 )
=Σk{-N(N-1)/2→N(N-1)/2}( ( k2 + 2Nk + N2 )PN-1(k) / 2 )
=( VN-1 + N2 ) / 2

ここで、Σk{-N(N-1)/2→N(N-1)/2}( k・PN-1(k) ) は期待値でゼロになるためこの項は消去できることに注意してください。後半の項は、PN-1(Δs+N) が | Δs + N | ≤ N( N - 1 ) / 2 で正値となるのでそれ以外の範囲は無視できて、Δs + N = k とすれば

ΣΔs{-N(N+1)/2→N(N+1)/2}( Δs2PN-1(Δs+N) / 2 )
=Σk{-N(N-1)/2→N(N-1)/2}( ( k - N )2PN-1(k) / 2 )
=Σk{-N(N-1)/2→N(N-1)/2}( ( k2 - 2Nk + N2 )PN-1(k) / 2 )
=( VN-1 + N2 ) / 2

従って、

VN = VN-1 + N2

という非常に単純な漸化式が得られます。これは、VN が自然数の二乗和を使って表されることを意味しており、その初項は V1 = 1 なので、結局自然数の二乗和そのものになって、( 2N3 + 3N2 + N ) / 6 になります。二乗和の一般式の解き方は、以前「数値積分法 -1-」で紹介した「ベルヌーイ数」の章の中で説明してありますので、そちらをご覧ください。

上記は、正値と負値の差に対する分布の漸化式となっていますが、一般的には正値と負値のうち小さい方を選び、その確率密度を使って検定を行います。正値の順位和を sp で表すと、負値の順位和は N( N + 1 ) / 2 - sp となり、この差が Δs であることから

Δs = sp - { N( N + 1 ) / 2 - sp }

という関係式が得られます。これを sp について解くと

sp = Δs / 2 + N( N + 1 ) / 4

となるので、Δs の分布に対して N( N + 1 ) / 4 シフトした上で幅を半分にした状態が sp の分布になります。左右対称の形は変わらないので、正値の順位和の方が大きくなった場合、分布の半分(順位和の小さい側)だけを使うことによって、負値の順位和を使っても推定値や検定の結果は変わりません。しかし、シフトしたことによって期待値は N( N + 1 ) / 4 に、分布の幅が半分になったことから分散も 1 / 4 の ( 2N3 + 3N2 + N ) / 24 になります。

U 分布の場合と同様に、順位和の分布も N が大きくなると正規分布に近似することができます。N( N + 1 ) / 2 通りある正値と負値の組み合わせが全て同等に発生するとして、その順位の和として(そして負値の場合は符号を負に変えて) Δs を求めているので、中心極限定理から直感的にイメージできると思います。平均はゼロ、分散は ( 2N3 + 3N2 + N ) / 6 なので、近似式は N( 0, ( 2N3 + 3N2 + N ) / 6 ) になります。


順位和とその確率密度を求めるためのサンプル・プログラムを以下に示します。

// 絶対値での比較用関数オブジェクト
class CmpByAbs {
public:
  // 二つのデータの比較
  bool operator()( const double& x, const double& y )
  { return( fabs( x ) < fabs( y ) ); }
};

// 符号付き順位和に対する確率密度を保持する配列用Mapの型
typedef map< unsigned int, map<int,double> > SRankMap;

/*
  calcSRankValue : 符号付き順位和を求めて返す

  const vector<double> &x, &y : データ列
  double& sRank : 求めた順位和
  unsigned int& cnt : 差がゼロ以外の要素数

  戻り値 : True ... 正常終了 ; False ... x と y のサイズが不一致、要素数がゼロ
*/
bool calcSRankValue( const vector<double>& x, const vector<double>& y, double& sRank, unsigned int& cnt )
{
  const double epsilon = 1e-12; // 絶対値を比較するときの精度

  if ( x.size() != y.size() ) return( false );
  if ( x.size() == 0 ) return( false );

  // x, y の差を計算
  vector<double> d; // 差分 x - y
  for ( unsigned int i = 0 ; i < x.size() ; ++i ) {
    double diff = x[i] - y[i];
    if ( diff != 0 ) d.push_back( diff );
  }
  cnt = d.size();

  // 差の絶対値から昇順で並べる
  CmpByAbs cmp;
  std::sort( d.begin(), d.end(), cmp );

  // 各要素の順位を計算
  vector<double> rank( d.size() ); // 差の順位
  for ( unsigned int i = 0 ; i < rank.size() ; ) {
    // 絶対値が等しいグループの末尾を探索
    unsigned int j = i + 1;
    for ( ; j < rank.size() ; ++j ) {
      if ( fabs( fabs( d[i] ) - fabs( d[j] ) ) > epsilon * fabs( d[i] ) ) break;
    }
    // 順位和の計算(値が等しい要素は等分配)
    double rankSum = 0; // 順位和
    for ( unsigned int k = i ; k < j ; ++k )
      rankSum += k + 1;
    rankSum /= j - i;
    // 順位和を配列に代入
    for ( ; i < j ; ++i )
      rank[i] = rankSum;
  }

  // 正値に対する順位和の計算
  double sp = 0;
  for ( unsigned int i = 0 ; i < d.size() ; ++i )
    if ( d[i] > 0 ) sp += rank[i];

  // Δs = 正値の順位和(sp) - 負値の順位和(sm) を計算
  // sp + sm = cnt( cnt + 1 ) / 2 より
  // Δs = sp - sm = sp - { cnt( cnt + 1 ) / 2 - sp }
  //    = 2sp - cnt( cnt + 1 ) / 2
  sRank = sp * 2 - (double)( d.size() * ( d.size() + 1 ) ) / 2.0;

  return( true );
}

/*
  sRankProb : 符号付き順位和 Δs の確率分布 Pn(Δs) を求める

  unsigned int n : データ列の要素数
  int sRank : 順位和 Δs

  戻り値 : Pn(Δs)
*/
double sRankProb( unsigned int n, int sRank )
{
  if ( n == 0 ) return( 0 );
  if ( n == 1 ) {
    if ( abs( sRank ) == 1 )
      return( 0.5 );
    else
      return( 0 );
  }

  if ( abs( sRank ) > n * ( n + 1 ) / 2 )
    return( 0 );

  return( ( sRankProb( n - 1, sRank - n ) + sRankProb( n - 1, sRank + n ) ) / 2.0 );
}

/*
  sRankGet : SRankMapから符号付き順位和に対する確率を取得する

  const SRankMap& sRankMap : 符号付き順位和に対する確率を保持したSRankMap
  unsigned int n : 確率を取得する対象のキー(データ列の大きさ)
  int sRank : 符号付き順位和 Δs

  戻り値 : Pn(Δs) (未登録ならば負数を返す)
*/
double sRankGet( const SRankMap& sRankMap, unsigned int n, int sRank )
{
  int sRankMax = n * ( n + 1 ) / 2;
  if ( sRank > abs( sRankMax ) ) return( 0 );

  SRankMap::const_iterator cit = sRankMap.find( n ); // 対象キーに対する sRankMap の反復子
  if ( cit == sRankMap.end() )
    return( -1 );

  const map<int,double>& sMap = cit->second; // 値を取得する対象配列
  map<int,double>::const_iterator sPair = sMap.find( sRank );
  if ( sPair == sMap.end() )
    return( -1 );

  return( sPair->second );
}

/*
  sRankSet : SRankMapに符号付き順位和に対する確率をセットする

  SRankMap& sRankMap : 符号付き順位和に対する確率を保持するSRankMap
  unsigned int n : 確率をセットする対象のキー(データ列の大きさ)
  int sRank : セットする対象の符号付き順位和 Δs
  double d : セットする値(確率)

  戻り値 : True ... セット成功 ; False ... sRankが範囲外
*/
bool sRankSet( SRankMap& sRankMap, unsigned int n, int sRank, double d )
{
  int sRankMax = n * ( n + 1 ) / 2;
  if ( sRank > abs( sRankMax ) ) return( false );

  map<int,double>& sRankDist = sRankMap[n];
  sRankDist[sRank] = d;

  return( true );
}

/*
  sRankProb : 符号付き順位和 Δs の確率分布 Pn(Δs) を求める

  一度再帰的に求めた確率は、SRankMapにセットして再利用できるようにする
  SRankMap 内の値はあらかじめ負数で初期化しておく必要がある

  int n : データ列の要素数
  int sRank : 符号付き順位和
  SRankMap& sRankMap : 値を保持するSRankMap

  戻り値 : Pm,n(U)
*/
double sRankProb( unsigned int n, int sRank, SRankMap& sRankMap )
{
  double d = sRankGet( sRankMap, n, sRank ); // 最初に SRankMap から確率を取得

  // 負数ならばまた求められていないので計算を行う
  if ( d < 0 ) {
    d = sRankProb( n, sRank );
    sRankSet( sRankMap, n, sRank, d ); // 求めた値を SRankMap にセットしておく
  }

  return( d );
}

/*
  sRankProb : 符号付き順位和 Δs の確率分布 Pn(Δs) を求めて Map の形で返す

  一度再帰的に求めた確率は、SRankMapにセットして再利用できるようにする

  unsigned int n : データ列の要素数
  map<int,double>& sRankDist : 結果を保持する Map (確率密度)
*/
void sRankProb( unsigned int n, map<int,double>& sRankDist )
{
  // 確率分布を保持する SRankMap
  SRankMap sRankMap;

  // 確率を計算して代入
  int sRankMax = n * ( n + 1 ) / 2;
  for ( int sRank = -sRankMax ; sRank <= sRankMax ; ++sRank )
    sRankDist[sRank] = sRankProb( n, sRank, sRankMap );
}

/*
  sRankProb : 0〜nMax までの順位和 Δs の確率分布 Pn(Δs) を求めて Map の形で返す

  unsigned int nMax : データ列の要素数最大値
  SRankMap& sRankMap : 結果を保持するMap
*/
void sRankProb( unsigned int nMax, SRankMap& sRankMap )
{
  // 確率分布を保持する SRankMap の初期化
  sRankMap.clear();

  // 確率を計算して代入
  for ( unsigned int n = 0 ; n <= nMax ; ++n ) {
    int sRankMax = n * ( n + 1 ) / 2;
    for ( int sRank = -sRankMax ; sRank <= sRankMax ; ++sRank )
      sRankProb( n, sRank, sRankMap );
  }
}

calcSRankValue は、与えられた二つのデータ列から符号付き順位和を求めるための関数です。求めた順位和の他に、差分がゼロではない要素の数も返すことができるように、これらはそれぞれ変数 sRankcnt への参照が引数として渡されています。コメントを見ながら内容を確認すれば、処理の流れは簡単に理解できると思います。
CmpByAbs は、絶対値を比較するための関数オブジェクトです。差の絶対値でソートする必要があるため、STL の中の std::sort を使うときに比較用関数オブジェクトとして利用しています。
絶対値が等しいかどうかの判定は、単純に等号で比較するのではなく、しきい値 epsilon に片側の値を掛けたものを使って差がその値より小さい場合に「等しい」としています。浮動小数点数の場合は演算結果にほぼ必ず誤差が生じるため、絶対値が等しくてもプログラム上では「等しくない」と判断される場合があります。そのため、十分に近い値ならば等しいとすることでこの問題を回避しています。

U 値の分布を求める関数 uProb のときと同様に、符号付き順位和を求める関数 sRankProb も四種類あり、最初の関数はデータ列の要素数 n と符号付き順位和 sRank からその確率密度を返すためのもので、二番目は最初の関数を利用して確率密度を求め、一度求めた値は SRankMap 型 ( 実体は map< unsigned int, map<int,double> > 型 ) の配列に保持して再利用できるようにした関数です。この二番目の関数を使い、三番目の関数は要素数が n のときの確率密度を全て求めるため、また最後の関数は要素数が 0 から nMax までの全ての確率密度を求めるために利用します。

下図は、要素数が 12 のときの符号付き順位和の確率分布を示したもので、緑色の曲線が正規分布曲線を表しています。

符号付き順位和の確率分布(n=12)

実際には、正規分布は実際の値を二倍しています。符号付き順位和に対する確率密度が値を持つのは必ず一つおきであり、連続分布に置き換える場合は全体の密度を半分にする必要があります。確率密度がゼロになる順位和に半分だけ確率密度を分け与えることをイメージすれば分かりやすいと思います。グラフ上では、そのような処理をする代わりに正規分布の方を二倍に大きくしてつじつまを合わせているわけです。なお、正規分布で近似する場合、この現象は特に意識する必要はありません(すでに順位和が連続値であるとして計算されるので、二倍するようなことは不要です。そんなことをしたら、全体の確率が 2 になってしまいます)。


ウィルコクソンの符号順位検定を行うためのサンプル・プログラムを以下に示します。

/*
  SRankTest : ウィルコクソンの符号順位検定

  const vector<double> &x, &y : データ列
  double a : 危険率

  戻り値 : true ... 帰無仮説(分布が等しい)が棄却された ; false ... 帰無仮説が保留された
*/
bool SRankTest( const vector<double>& x, const vector<double>& y, double a )
{
  cout << "***** Wilcoxon Signed-rank Test *****" << endl << endl;

  double sRank;   // 符号付き順位和
  unsigned int n; // 有効データ数
  if ( ! calcSRankValue( x, y, sRank, n ) )
    return( NAN );

  // 符号付き順位和と分布の大小関係を出力
  cout << "Signed-Rank = " << sRank << endl;
  if ( sRank >= 0 ) {
    sRank *= -1;
    cout << "Distribution X seems to be more than distribution Y." << endl << endl;
  } else {
    cout << "Distribution Y seems to be more than distribution X." << endl << endl;
  }

  double p = 0; // 下側確率
  if ( n <= 10 ) {
    map<int,double> sRankDist;
    sRankProb( n, sRankDist );
    int i = -( n * ( n + 1 ) / 2 );
    for ( ; (double)i < sRank ; ++i )
      p += sRankDist[i];
    p += ( sRank - (double)( i - 1 ) ) * sRankDist[i];
    p *= 2.0;
  } else {
    NormalDistribution nd( 0, sqrt( n * ( n + 1 ) * ( 2 * n + 1 ) / 6.0 ) );
    p = nd.lower_p( sRank ) * 2.0;
  }

  // p-値と検定結果の出力
  cout << "p-value = " << p << endl;
  if ( p < a ) {
    cout << "Null hypothesis is rejected." << endl;
    return( true );
  } else {
    cout << "Null hypothesis is accepted." << endl;
    return( false );
  }
}

U 検定の場合とほとんど同じような処理を行っています。符号付き順位和を求めた上で、その正負からどちらの分布がより大きいかを判定し、またその絶対値を使って差に有意性があるかを検定します。有効な(差がゼロ以外の)要素数が 10 以下であれば直接確率密度を計算し、10 より大きい場合は正規分布を使って近似を行います。


U 分布と符号付き順位和の分布には密接な関係があります。N 個の差分に対して正値が m 個あった場合、負値は N - m 個で表され、それらの順番によって U 値と符号付き順位和が求められます。正値と負値の数を固定したとき U 値の分布が決まり、正値の個数がゼロから N まで変化したもの全体を使って符号付き順位和の分布が表されるので、符号付き順位和の分布の中には複数の U 分布が含まれることになります。簡単な例として、N = 3 の場合を考えると、

データ列正値の数負値の数U値順位和
[-][-][-]030-6
[+][-][-]120-4
[-][+][-]121-2
[-][-][+]1220
[+][+][-]2100
[+][-][+]2112
[-][+][+]2124
[+][+][+]3006

となって、P0,3(U), P1,2(U), P2,1(U), P3,0(U) の 4 つの U 分布から符号付き順位和の分布 P3(Δs) が構成されていることがわかります。正値・負値の要素数を m, n (このとき、N = m + n が成り立ちます)、順位和をそれぞれ sp, sm としたとき、U と順位和の関係から

U = mn + n( n + 1 ) / 2 - sm

が成り立つことを前節で説明しました。また、符号付き順位和 Δs は

Δs = sp - sm

で表され、sp と sm の和が N( N + 1 ) / 2 と等しいことから

N( N + 1 ) / 2 = sp + sm

でもあるので、これら二つの式から

Δs = N( N + 1 ) / 2 - 2sm

になります。この式と、U に関する式を、sm が消去されるように辺々足し合わせると

2U - Δs=2mn + n( n + 1 ) - N( N + 1 ) / 2
=2m( N - m ) + ( N - m )( N - m + 1 ) - N( N + 1 ) / 2
=2mN - 2m2 + N2 - 2mN + m2 + N - m - N2 / 2 - N / 2
=N2 / 2 + N / 2 - m2 - m
=N( N + 1 ) / 2 - m( m + 1 )

という関係式が得られるので、正値(または負値)の数が決まれば、U 値から符号付き順位和が、またその逆が求められることになります。


今回は、順序統計量をテーマとして取り上げてみたのですが、予想以上に奥が深く、分量も想定していたものより多くなりました。調べればさらに面白いテーマが見つかりそうな気もします。
データが正規分布に従うと仮定した場合、今回の検定手法を用いる必要はないわけで、前回紹介した手法のほうがより精度よく定量的に結果を得ることができます。しかし、データ数が少ないなどの理由で正規分布に従うかどうかわからない、また正規分布に明らかに従わないような場合は、今回紹介した検定手法によって結果を得ることは可能です。結局はデータを見て判断することになるので、やはりヒストグラムなどのグラフを確認した上で検定を行うようにして、ただやみくもに手法だけを利用するようなことは避けるようにしましょう。


<参考文献>
1. 「統計数学入門」 本間鶴千代著 (森北出版)
2. 「日本女子大学 理学部 数物科学科」-「今野良彦のホーム」-「順序統計量
上の二つは、順序統計量の確率分布に関する部分について参考としました。
3. On a Test of Whether One of Two Random Variables is Stochastically Larger than the Other by H.B.Mann and D.R.Whitney
マンとホイットニーによる 1947年の論文です。
4. 「Vassar College」-「Concepts & Applications of Inferential Statistics」-「The Wilcoxon Signed-Rank Test
今回参考にしたのは「符号付き順位和」に関する箇所のみですが、他にも有用な情報がいろいろとあります。
5. 我楽多頓陳館 --- 統計学入門
6. Wikipedia

◆◇◆更新履歴◆◇◆


[Go Back]前に戻る [Back to HOME]タイトルに戻る
inserted by FC2 system