せいびし

専門性を持たない自分の正体を求めて

【各種検定の理論】順序統計量概論(ノンパラメトリック検定)

導入

ノンパラメトリック検定は、測定したデータの性質が質的変数であったり、正規分布に従わないデータに使われます。ノンパラメトリック検定の根底に存在する仕組みはデータの順番です。僕自身、ノンパラメトリック検定の存在は知っていたのですが、各検定がどういう仕組みで動いているのかを全く知らなかったので、これを機会に勉強しています。

各検定の仕組みに入る前に重要な話として順序統計量が存在しています。順序統計量は手に入れたデータを単に昇順に並べただけのものですが、仕組みが簡単な割に数式が爆発していてとても面白い事になっています。座学だけなのもつまらないので、手を動かしてみましょう。

1.順序統計量の概要

順序統計量の定義

データの順番の考え方は、特定の確率分布に依存しません。 \displaystyle X_i(i=1..n)を任意の連続な分布関数 \displaystyle F(x)から得られる標本とします。ここで、 \displaystyle F(x)は次の式で定義されます。

 \displaystyle F(x) = Pr(X_i \leq x)

ここから、順序統計量 \displaystyle X_{(i)}を次のように定義します。
 \displaystyle X_{(i)}:標本を昇順に並べた時の \displaystyle i番目のデータ

この定義から、最小値は \displaystyle X_{(1)}、最大値は \displaystyle X_{(n)}となります。

順序統計量の確率分布

導出をすると長くなるので、すっ飛ばして式だけ書いておきます(考え方は注釈に回します*1*2*3)。

 \displaystyle F_{(i)}(x) =  \sum_{k=i}^{n}{}_n C_k F(x)^k(1-F(x))^{n-k}

 \displaystyle f_{(i)}(x) = \frac{n!}{(i-1)!1!(n-i)!}F(x)^{i-1}(1-F(x) )^{n-i}f(x) ( -\infty \leq x \leq \infty )

1つ1つの量に確率分布が存在するのはなんだか不思議な感じがしますね。

順序統計量ではよく中央値を求める事があります。中央値はデータの個数が奇数の時はちょうど真ん中の値であり、偶数の時は最小値と最大値から辿ってちょうど交差する2つの値で平均を取るものでした。中央値を \displaystyle X_{(M)}とした時、データの個数が偶数、奇数の時の中央値の確率密度関数は上の公式から次のように導出する事が出来ます。

  • データの個数が偶数の時

データの個数 \displaystyle nが偶数の時は2つの順序統計量を用いるため、順序統計量の同時確率密度関数から求めます。式と考え方は注釈*4に載せておきます。

 \displaystyle f_M(x) = \frac{n!}{\{(\frac{n}{2}-1)!\}^2} (F(x_1)(1-F(x_2))^{\frac{n}{2}-1}f(x_1)f(x_2) (-\infty < x_1 < x_2 < \infty)

  • データの個数が奇数の時

データの個数 \displaystyle nが奇数の時はちょうど真ん中の値になるため、1つの順序統計量の確率密度関数から求めます。

 \displaystyle f_M(x) = \frac{n!}{\{(\frac{n-1}{2})!\}^2}(F(x)(1-F(x)))^{\frac{n-1}{2}}f(x) (-\infty < x < \infty)

順序統計量のモーメント母関数

モーメント母関数の導出は割と簡単です。例えば、連続型確率変数 Xの原点周りのr次モーメントは次の式で表されました。

 \displaystyle E[X^r] = \int_{-\infty}^{\infty}x^rf(x)dx

この定義に順序統計量を丸々当てはめると、
 \displaystyle E[X_{(i)}^r] = \frac{n!}{(i-1)!1!(n-1)!} \int_{-\infty}^{\infty} x^r \cdot F(x)^{i-1}(1-F(x) )^{n-i}f(x) dx

あっさり終わってしまいましたね

2.シミュレーション

やること

個々の順序統計量に分布が存在する事が分かりました。モーメントも導出したので、期待値と分散も定義する事が出来ます。本記事のシミュレーションでは、順序統計量として、最小値、中央値、最大値の分布を可視化してみたいと思います。シミュレーションする確率密度関数は連続一様分布 \displaystyle U(0, 1)とします。実験の手順を図で示します。

f:id:seibibibi:20191105163820p:plain

一様分布から \displaystyle n個(奇数個)のデータをサンプリングし、最小値、中央値、最大値を求めます。この手順を複数回行い、それぞれの順序統計量に関するヒストグラムを作成します。また、この分布では \displaystyle X_{(i)}の期待値と分散は次の式で表されます。ここで、 \displaystyle M=\frac{n+1}{2}です。*5

 \displaystyle E[ X_{(i)} ] = \frac{i}{n+1}, V[ X_{(i)} ]=\frac{i(n-i+1)}{ (n+1)^2 (n+2) }

この事から、最小値 \displaystyle X_{(1)}、中央値 \displaystyle X_{(M)}、最大値 \displaystyle X_{(n)}の期待値と分散は次の式で表すことができます。
 \displaystyle E[ X_{(1)} ] = \frac{1}{n+1}, V[ X_{(1)} ]=\frac{n}{ (n+1)^2 (n+2) }

 \displaystyle E[ X_{(M)} ] = \frac{1}{2}, V[ X_{(M)} ]=\frac{1}{4(n+2)}

 \displaystyle E[ X_{(n)} ] = \frac{n}{n+1}, V[ X_{(n)} ]=\frac{n}{ (n+1)^2 (n+2) }

実験結果

実験1 とりあえず回してみる

まずは適当なパラメータを設定して回してみます。 \displaystyle n=9とし、サンプリング回数を100回にしました。ヒストグラムと密度関数を以下に示します。

f:id:seibibibi:20191105161449p:plain

見た感じ中央値のヒストグラムが収束していないような気がしますね。ちなみにこの状況で各統計量の平均は \displaystyle 0.09, 0.51, 0.90、分散は \displaystyle 0.0077, 0.023, 0.0094となりました。期待値と分散の式にそのまま代入すると、期待値は \displaystyle 0.10, 0.50, 0.90、分散は \displaystyle 0.0082,0.023,0.0082なので、期待値はほとんど一致している事が確認できました。今の所分散は多少の差が見られますが、多分シミュレーション回数を増やせば一致してくるのでしょう。でもこれで終わるのもつまらないので、もう少し遊んでみます。

実験2 1回当たりのサンプリング数 \displaystyle nを変える

次に、1回当たりのサンプリング数を適当に \displaystyle n=5, 25, 45にした時のヒストグラムを見てみましょう。サンプリング回数は100回です。

f:id:seibibibi:20191105162431p:plain

サンプリング数が多いと分散が小さくなるので、期待値周辺に値が集中するようになっていますね。面白いなあ。

実験3 サンプリング回数を変える

最後は、サンプリング回数を変えてみましょう。 \displaystyle n \displaystyle 15とし、サンプリング回数を10回、100回、1000回と変えていきます。ヒストグラムは以下です。

f:id:seibibibi:20191105163102p:plain

サンプリング数は変えていないので、分散は変わらず、単に山が滑らかになりつつあるというかなんというか…

考察

そりゃそう

まとめ

今回は順序統計量で遊んでみました。当たり前のことですが、シミュレーションをしてもきちんと理論通りの値を吐き出してくれます。逆に考えれば、こうやっていちいちシミュレーションをしなくても解析的に解が求められる事がいかに面白いかという事を再発見しました。現実世界ではそう何回もデータを取っては推定量出して、またデータを取っては推定量を出して…なんてできないので、何か理論上の方針というか目安というか…そういう物があると考え方として面白いですよね。以上でした!


※今回作成したプログラムはこちら
seibibibi.hatenadiary.jp

*1: \displaystyle F_{(i)}(x) = Pr(X_{(i)} \leq x)から始まります。 \displaystyle X_{(i)} \leq xという事は、 \displaystyle X_{(i)}さえx以下であれば良い訳で、 \displaystyle X_{(i+1)}より先は大きかろうが小さかろうが関係ない訳です。言い換えれば、大きかろうが小さかろうがの全ての事象を和の事象として足し合わせなければなりません。任意の順序統計量がx以下である確率は \displaystyle F(x)であり、xより大きい確率は \displaystyle 1-F(x)となります。後は何個中何個がx以下で…と考えれば、この考え方は二項分布になるため、結局二項分布の足し合わせになります。

*2:確率密度関数を求める時は、 \displaystyle F(x)=pとして、 \displaystyle \frac{df_{(i)}}{dx} = \frac{df_{(i)}}{dp} \cdot \frac{dp}{dx}で導出します。後はシグマの中で打ち消し合いが起きて要らないものが勝手に消えてくれます。

*3:確率密度関数の導出は多項分布から導出する事も出来ます。多項分布は次の式で表されます。

 \displaystyle f(x_1, ..., x_k|n, p_1, ..., p_k) = \frac{n!}{x_1!...x_k!}p_1^{x_1}...p_k^{x_k}
多項分布の簡単なイメージは多面サイコロです。各面で出る確率 \displaystyle p_i (i=1..k)が割り振られており、 \displaystyle n回振ると各面が何回出てきますか?の分布です。これを順序統計量に応用します。順序統計量の個数を \displaystyle n \displaystyle -\infty < x < \infty上の \displaystyle x_i < x < x_i+\delta x区間 \displaystyle X_{(i)}が1つだけ存在しているものとし、この時の \displaystyle Pr(x < X_{(i)} < x+\delta x)について考えてみます。多面サイコロで考えれば、確率 \displaystyle F(x)で生起する \displaystyle X_{(i)}以下の面、確率 \displaystyle F(x+\delta x)-F(x)で生起する \displaystyle X_{(i)}の面、確率 \displaystyle 1-F(x+\delta x)で生起する \displaystyle X_{(i)}より大きい順序統計量の面の3面に分ける事が出来ます。 \displaystyle n個の順序統計量の内、 \displaystyle X_{(i)}以下となる順序統計量の個数は \displaystyle i-1個、 \displaystyle X_{(i)}より大きい順序統計量の個数は \displaystyle n-i個になる事は、多面サイコロで言い換えると、 \displaystyle n回振った時に、それぞれ \displaystyle i-1回、1回、 \displaystyle n-i回出てくることになります。従って、 \displaystyle Pr(x < X_{(i)} < x+\delta x)は多項分布の式に当てはめると次のように表されます。
 \displaystyle \frac{n!}{(i-1)!1!(n-i)!}F(x)^{i-1} ( F(x+\delta x)-F(x) ) ( 1-F(x+\delta x) )^{n-i}
この式を \displaystyle \delta x→0導関数を求めてあげれば、
 \displaystyle \lim_{\delta x \to 0} \frac{n!}{(i-1)!1!(n-i)!}F(x)^{i-1} ( 1-F(x+\delta x) )^{n-i} \frac{F(x+\delta x)-F(x)}{\delta x}
 \displaystyle = \frac{n!}{(i-1)!1!(n-i)!}F(x)^{i-1}(1-F(x) )^{n-i}f(x)
となります。個人的には、こっちの導出方法の方が順序統計量の考え方も網羅されていて好きです。順序統計量の同時確率密度関数を導出する時も多項分布を利用すると全く同じ考え方で導出する事が出来ます。同時確率密度関数の場合は多面サイコロで言えば5つの面が現れます。

*4:という事で、多項分布を用いた導出を簡単に載せておきます。同時確率密度関数の場合は

 \displaystyle Pr( x_i < X_{(i)} \leq x_i+\delta x, x_j < X_{(j)} \leq x_j+\delta j)
で考えます。多面サイコロで言う所の5つの面が見えましたでしょうか。あとは5つの区間に所属する順序統計量の個数と確率を求めて、多項分布に代入した物を \displaystyle \delta x, \delta y \to 0導関数を取るだけです。最終的に得られる同時確率密度関数は以下になります。
 \displaystyle C_{i, j} F(x_i)^{i-1}  ( F(x_j) - F(x_i ) )^{j-i-1} ( 1-F(x_j) )^{n-j} f(x_i) f(x_j)
 \displaystyle C_{i,j} = \frac{n!}{ (i-1)!(j-i-1)!(n-j)! }

*5:注釈欄が記事の半分くらいになってしまって笑ったので、導出は省略です。確率密度関数と確率分布関数を定義の式に突っ込むとベータ関数が現れるので、ベータ関数をガンマ関数で表してごにょごにょすると導出できます。