せいびし

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

サイトマップ

はじめまして

本記事はサイトマップとなっています。記事を更新次第、本記事も随時編集していきます。

せいびし(Twitter : @seibibibi)

僕は博士課程の学生です。僕は自分が一体何者なのかを分かっていません。何者かが分からないというのは、研究の領域以外で自分の興味の赴くままに色々と勉強していた所、手を広げすぎて収拾がつかなくなったという所です。興味の勉強と研究がきちんと分離されているならこのような事にはならなかったのですが、最近は入り混じる事が多くなり、本格的に自分が何をしているのか何も分からなくなってしまいました。一応の研究領域は生体医工学で、片麻痺患者を対象とした新しい手指運動機能評価スケールの提案を行っています。

自分が一体何者なのかを明らかにする事が本ブログの主な目的です。何に対しても全く詳しくない自分が、その時々で勉強した内容を寄せ集めた物が本ブログです。恐らく専門の方が内容を見た時には間違いが多々ある事が予想されます。その際はぜひ間違いを指摘してくださいますと、僕の励みにもなり、本当に嬉しく思います。ぜひ宜しくお願いします。(今までの記事をすべて削除して1から作り直しています。トピックはその時その時で小出しにしていきます。1つずつ工事中を開けていくのを楽しみながら頑張っていきたいと思います。(2019年11月1日現在))

トピック1.医療統計の手法としくみ

主に医療統計のトピックです。医療統計で用いられる手法の理論と実例の2つから構成されます。

理論では、あまり難しい数学の話はせず、医療統計で種々に存在する手法の数理統計的な性質のシミュレーションを取り扱いたいと思います。読まれる対象は、今まで医療統計を使ってきたけど数理的な話は何だか良く分からない人数理的な部分に興味が出てきた人など、ある程度医療統計を用いてきた人が対象になるのかなと思います。

実例では、主にRを用いた解析を行います。数理的な所はひとまずおいといて分析と解釈を重視したい人Rで医療統計をいじってみたい人はここから始まると見通しが立ちやすいのかなと思います。

理論はもちろんですが、実例でも基本的な数理統計学の話は出てきます。今までに統計の理論的な話に触れたことが無かった人向けに、本トピックでは数理統計学の基本という内容を用意しておきました。数理統計はかなり魔境の分野ですが、医療統計では用いられる手法や考え方が割と限られているので、医療統計において最低限必要な内容をまとめておきました。ここを読まれると、恐らくこのトピックの全てがすすっと読めると思います。

僕も医療統計はまだまだ勉強中です。一緒に頑張っていきましょう!

数理統計学の基本

工事中…

各種検定

各種検定の理論

主に、1標本の平均値の差の検定、2標本の平均値の差の検定、ノンパラメトリック検定、その他の検定を扱います。
seibibibi.hatenadiary.jp

Rで各種検定

工事中…

分散分析の理論

3標本以上の平均値の差の検定について取り扱います。各種検定の理論に入れても良かったのですが、かなり長くなりそうなので分けました。
seibibibi.hatenadiary.jp

Rで分散分析

工事中…

回帰分析

回帰分析の理論

工事中…

Rで回帰分析

工事中…

書籍

書籍はこちらです。
seibibibi.hatenadiary.jp

カーネルロジスティック回帰で遊ぶ(FUN Advent Calendar #15)

自己紹介

初めましての方も多いと思いますので、自己紹介から始めます。整備士(Twitter:@seibibibi)と申します。博士(後期)課程の2年生です。


学内の一部の人からは指の人と呼ばれている(いた)ようです*1。院棟の奥深くに案内され、謎の手のロボットをいじった事があるぞという方はご存知だと思います。


ヒトの身体に興味があります。特に、ヒトはどのように病気に罹り回復していくのかをモデリングする事に強い興味を持っています。個人的には予後予測リハビリの計画などにほんの少しだけ役に立てたら嬉しいなと思っています。


…という研究を行うためには中々難しい所があります。なので、この目標の似たような試みとして、現在は片麻痺を対象とした手指リハビリ機器による手指運動機能の評価スケールの提案を本学で行っています。専門領域はありませんが、研究では統計と機械学習でごにょごにょしています。


本記事の内容

自己紹介もほどほどにして、本題に入っていきましょう。研究で機械学習を使っていると述べましたが、高尚なモデルは使っていません。主にロジスティック回帰モデルという物を使っています。そこで、そろそろ自分の番だし何書こうかなあ…ロジスティック回帰でなんか書けるかなあと思いつつ調べ物をしていたら、こんな面白いものを見つけました。

ー sklearnライブラリを使用してロジスティック回帰モデルでカーネルを使用する方法を教えてください。
とても素敵な質問ですが、現在scikit-learnはカーネルロジスティック回帰もANOVAカーネルもサポートしていません。

機械学習 – ロジスティック回帰モデルのカーネルLogisticRegression scikit-learn sklearn - コードログ


英語をそのまま機械翻訳したかのようなとても素敵な回答を見つけました。ロジスティック回帰は非常に素朴なモデルであり、真っ直ぐに識別面を引っ張ってくれます。しかし、世の中の多くのデータは簡単に線形分離できるものばかりではありません。こういう時にはロジスティック回帰はあまり有効な手段とは言えません。


カーネルロジスティック回帰は、ロジスティック回帰がベースですが、データに対して後から述べるカーネルなるものを掛ける事でぐにゃっと曲がった識別面を引っ張ってくれるとても素敵なモデルとなっています。そんなものが実装されていないなんて!なんてことだ!


というわけで、とても素敵なモデルなのにどうして!どうして!どうして無いの!なノリで、即席で弾丸実装をしようというのが本記事の内容です。ロジスティック回帰の細かい話は、専門書やネットを引っ張れば大量に出てくるので、そちらに譲る事にします。また、様々な学年の方が読まれると思いますので、数式は1つも出さないようにしました*2。しかし、ロジスティック回帰をカーネル化するにあたってどうしてもここだけは…!という所は、ごめんなさいですが式を載せさせて頂きました。ごめんなさい。

1.ロジスティック回帰の雰囲気

シグモイド関数

まずはシグモイド関数なるものを示します。シグモイド関数はロジスティック回帰を形作る非常に重要な関数です。システムへの入力を \displaystyle xとする時、シグモイド関数は次のような形になります。

f:id:seibibibi:20191208222805p:plain

横軸は入力、縦軸は確率になっています。確率は0から1の間を取るので、入力に対応する縦軸の値もまた0から1の間に収まっている事が分かります。

シグモイド関数は確率を返しますが、識別の時に僕たちが欲しいのは0か1、真か偽、陽性か陰性か…などの2値です。このような場合にはどう考えたら良いのでしょうか。この仕組みは非常に簡単です。次の図を見てみましょう。

f:id:seibibibi:20191210161804p:plain

縦軸を"1"の出力が得られる確率と設定すれば、0.5を境にしてどちらを出力するかを考える事が出来ます。


まずはロジスティック回帰で線を引っ張ってみよう!

ダミーデータを用意し、ロジスティック回帰で識別してみました。今度は入力を \displaystyle x_1, x_2の2つにしてみました。以下がその様子です。

f:id:seibibibi:20191208230304p:plain

真っ直ぐに識別面を引っ張っており、赤い集団と青い集団を上手く分離する事が出来ています。もう1つ見てみましょう。

f:id:seibibibi:20191208232200p:plain

この場合は、ロジスティック回帰は全く役に立たない事が分かります。ロジスティック回帰はその素朴さがゆえに、非線形の識別面を要する場面ではあまり有効ではありません。でも、そんなピンチを救ってくれるテクニックがあります。それがカーネルです*3

2.ロジスティック回帰をカーネル化する

30秒で知った気になるカーネルの雰囲気

いよいよカーネルの話に入っていきます。先ほどの、ロジスティック回帰では太刀打ちできなかったデータセットについて考えてみます。次の図を見てみましょう。

f:id:seibibibi:20191208234146p:plain

こんなことができたらとても嬉しいですよね。入力から更に何かしらを作りだし、それらが張る空間の中でなら素朴なロジスティック回帰が使えるかもしれない…!という発想です。

カーネルを使った回帰

カーネルを使う事でどんなことが起きるのか、基本的な線形回帰と比較する事にします。そうですね、重回帰モデルだとイメージがしやすいのではないでしょうか。

重回帰モデル

重回帰モデルは、例えば筋力を予測するモデルを考えた時に筋力=重み1×身長+重み2×体重+…と、予測する変数に対して特徴量(ここでは身長、体重、…)の1次式で結合していくモデルです。ここでの特徴量は1次なので、変数が1つなら直線、2つなら平面…と、いずれにしてもまっ平らな何かしらが空間に引っ張られます。ちょっと可視化してみましょうか。変数が2つの時のモデルはこんな感じです。

f:id:seibibibi:20191209000005p:plain

カーネルを使った回帰

一方、カーネルを使った回帰では、特徴量に少し小細工を加えます。どのような小細工を加えるかというと、特徴量を高次元の特徴空間へ写像します。特徴空間…?写像…?となんもわからんになってしまった場合は、30秒で知った気になるカーネルの雰囲気で載せた図をイメージしてみて下さい。ここで、高次元の特徴空間上では平面でスパッと切る事が出来ましたが、この識別面を元の特徴空間で見てみるとどうなるでしょうか。ぐにゃぐにゃ曲がって見えるのではないでしょうか。識別は実装で行うとして、重回帰モデルで説明した特徴量に、とあるカーネルを掛けた時の結果を示してみます。

f:id:seibibibi:20191209020344p:plain

うおおなんだこれはああああ…という感じですよね。こんなへんてこな形でも、ちゃんと回帰になっている訳です。

データを与えてフィットしてみよう!

今までは重回帰モデルとカーネルを用いた回帰の一般的な話をし、可視化をしてきました。しかし、これだけでは何だか雰囲気が掴めません。実際にデータを与えてフィットしてみましょう。

sinカーブを立体にしたデータセットに対して、まっ平らな平面を引っ張るモデルとぐにゃぐにゃを引っ張るモデルを当てはめます。今までまっ平らな平面のモデルは重回帰モデルで説明してきましたが、今回はリッジ回帰*4という物を用います。図を見てみましょう。

f:id:seibibibi:20191210125845p:plain

まっ平らな平面を引っ張るモデルでは表現力が足りず、いまいちデータの分布の特徴を掴み切れていない感じがします。一方で、ぐにゃぐにゃを引っ張るモデルでは、見事に波打った面を引っ張ってくれています。このように、カーネル化は引っ張る面の表現力を上げる事が出来ます。これは嬉しい!

ロジスティック回帰のカーネル

カーネルの導入

ここまでで、なんとなくカーネルとは何ぞや??という事の雰囲気が掴めたと思います。ここからは、ロジスティック回帰のカーネル化について説明していきます。少しだけ式が出てきます。

今まで説明してきたものは特徴量が2次元でしたが、ここでは一般化してn次元の特徴量 \displaystyle {\bf x} = (x_1, x_2, ..., x_n) を考えます。この時、出力を特徴量 \displaystyle {\bf x}と重みベクトル \displaystyle {\bf \omega} を用いて、何らかの式を当てはめる事を考えます。例えば、重回帰モデルの場合は \displaystyle {\bf x} \displaystyle {\bf \omega} を用いて、実数を取る出力 \displaystyle yを次の式で表現します。

 \displaystyle y = \omega_0 + \omega_1 x_1 + ... + \omega_n x_n

一方、ロジスティック回帰であれば、出力は \displaystyle yから \displaystyle p (0 \leq p \leq 1) となり、次の式で表現されます。

 \displaystyle \ln{ \frac{p}{1-p} } =  \omega_0 + \omega_1 x_1 + ... + \omega_n x_n

どちらも共通して \displaystyle \omega_0 + \omega_1 x_1 + ... + \omega_n x_nが現れている事が分かります。カーネル化はこの部分に小細工を加えます。カーネルは上で述べた通り、既にある特徴量を用いて特徴空間の次元をふわっと高次元に上げるイメージでした。この事は次の式で表現されます。

 \displaystyle \alpha_1 k( {\bf x^{ (1) }, x} ) + \alpha_2 k( {\bf x^{ (2) }, x} ) + ... + \alpha_N k( {\bf x^{ (N) }, x} )

…なんぞ??

何だか良く分からない…もっと簡単に…

それでは特徴量を \displaystyle xの1次元にしましょう。この時の重回帰モデル(この場合は単回帰ですが)とカーネル化した後の回帰を比較しましょう。

 \displaystyle y = \omega_0 + \omega_1 x

 \displaystyle y = \alpha_1 k(x^{ (1) }, x) + \alpha_2 k(x^{ (2) }, x) + ... + \alpha_N k(x^{ (N) }, x)

さて、どうでしょうか。単回帰の場合は、ある1つの \displaystyle yを求めるのに \displaystyle xしか使っていないという事が出来ます。しかし、カーネル化した場合は、1つの \displaystyle yを求めるために \displaystyle xから \displaystyle N個の \displaystyle k(x^{ (i) }, x)なる量を求めて、それらの重み付き総和を求めていると言い換える事が出来そうです。この時、 \displaystyle k(x^{ (i) }, x)が次の式だった場合は何が想像できるでしょうか。

 \displaystyle k(x^{ (i) }, x) = \exp{ (-\gamma (x^{ (i) } - x )^2 ) }

この \displaystyle k(x^{ (i) }, x)はガウシアンカーネルと呼ばれており、 \displaystyle x^{ (i) }によって中心の位置が変わるガウス分布の形を呈しています。ある \displaystyle xを入力した時に返ってくる \displaystyle k(x^{ (i) }, x)の値は中心の位置に依存する上に、ガウス分布の形は山なりなので、そんな総和を取ってしまえば線も面もぐにゃぐにゃになる事が容易に想像できると思います。

ああもうわからん!図で説明してくれ!

例えば、次の回帰式を想定します。

 \displaystyle y = \alpha_1 k(x^{ (1) }, x) + \alpha_2 k(x^{ (2) }, x) + \alpha_3 k(x^{ (3) }, x)

この式は、ある1つの \displaystyle yを求めるのに、 \displaystyle xから新たに3つの特徴量を抜き出していると言えそうです。ここで、カーネルはガウシアンカーネルを用いるものとし、 \displaystyle \gamma = 1.0, x^{ (1) } = 0.2, x^{ (2) } = 0.5, x^{ (3) } = 0.8とした時のカーネルと計算方法を示します。

f:id:seibibibi:20191210235204p:plain

いかがでしょうか?ここで説明したのは1次元の特徴量ですが、同じことが \displaystyle n次元の特徴量にも言えます。再度、 \displaystyle n次元の特徴量の場合のカーネル回帰を見てみましょう。

 \displaystyle \alpha_1 k( {\bf x^{ (1) }, x} ) + \alpha_2 k( {\bf x^{ (2) }, x} ) + ... + \alpha_N k( {\bf x^{ (N) }, x} )

イメージが湧きましたか?

やっとロジスティック回帰のカーネル

というわけで長々と説明してしまったわけですが、いよいよ本題に入ります。とは言ってもここまでわかってしまえば話は非常に単純なので、あっさり終わらせてしまう事にします。

通常のロジスティック回帰は次の式で表されました。

 \displaystyle \ln{ \frac{p}{1-p} } =  \omega_0 + \omega_1 x_1 + ... + \omega_n x_n

この式の右辺を次のように書き換えます。

 \displaystyle \ln{ \frac{p}{1-p} } = \alpha_1 k( {\bf x^{ (1) }, x} ) + \alpha_2 k( {\bf x^{ (2) }, x} ) + ... + \alpha_N k( {\bf x^{ (N) }, x} )

これで完成です。意外と簡単!


3.実装

というわけで実装してみましょう。ロジスティック回帰のカーネル化は先ほども述べた通り、 \displaystyle  \omega_0 + \omega_1 x_1 + ... + \omega_n x_nの所をカーネル化するだけなので、そこをカーネル化してsklearnのLogisticRegressionに突っ込んでいきたいと思います。ここで、用いるカーネルはガウシアンカーネルとし、識別で用いるデータは冒頭で述べた線形分離不可能なデータセットにします。そうすると、ハイパーパラメータはガウシアンカーネルに現れる \displaystyle \gamma正則化 \displaystyle Cの2つになります。先に実験結果を述べる前に断りたいのですが、僕のクソ雑魚プログラミング能力のせいで高速化までは間に合わなかったので、パラメータのチューニングは諦めました*5

実験結果1 XORパターンで実験

ロジスティック回帰で線を引っ張ってみよう!のところで太刀打ちできなかったあのデータセットは、XORパターンと呼ばれているようです。カーネル化したロジスティック回帰で識別してみましょう。結果は次のようになりました。ここで、 \displaystyle \gamma=1.0, C=1.0としました。クソ雑魚プログラムのせいで処理がめっちゃ重いので、メッシュは荒くしてあります(識別面が歪んで見えるのはごめんなさい)。

f:id:seibibibi:20191210124026p:plain

見事に分離できている事が分かります。直線で分離不可能な物も、カーネル化によって高次元の特徴空間の下ではまあまあ線形分離できた結果という事でしょう。

実験結果2 二重丸パターンで実験

XORパターンだけだと面白くないので、データが二重丸になっている時の場合も実験してみました。

f:id:seibibibi:20191215144637p:plain

左側は普通のロジスティック回帰で、ぴえんな結果になっていますが、右側は良い感じに分類できています。ハイパーパラメータは \displaystyle \gamma=0.1, C=1.0としておきました。良い感じですね。

4.おまけ

と、色々と述べては実験をしてみたわけですが、プログラムを組んでいる時に既存のライブラリにグラム行列Kをぶっこむだけで本当に良いの?正しいの?という疑問が生まれました。そこで、sklearnのリッジ回帰とカーネルリッジ回帰は一体どういう仕組みになっているのかなあとライブラリを少し読んでみました。なぜリッジ回帰を対象にしたかというと、リッジ回帰の重みベクトルは最小二乗法によって解析的に求める事ができ、どうせ最小二乗法のメソッドかなんかに突っ込んで終わるだろうから、引数を見れば分かるだろう……とか思ったからです。で、調べてみた結果ですが…

どちらもnumpy.linalg.lstsqなるものを使っているようです。ただし、if-elseだtry-catchだ色々と分岐が入っているので条件によって色々と変えているので、全部が全部こいつに突っ込んでいる訳でも無さそうですが、とりあえずこいつは一体何者を調べてみました。

Solves the equation  \displaystyle a x = b by computing a vector x that minimizes the squared Euclidean 2-norm  \displaystyle \| b - a x \|^2_2

docs.scipy.org


要するに最小二乗法をやってくれているっぽいですね。通常のリッジ回帰の場合はXとyを引数に渡し、カーネルリッジ回帰の場合はKとyを引数に渡しているだけの違いなようです。これなら手動でグラム行列を渡してやっても同じ結果が出る気がしますね*6


5.まとめ

今回はロジスティック回帰をカーネル化して遊んでみました。他の人が割と趣味の内容だったり本学の話だったりを書いている中で、なんだか小難しい話をしてしまいました。申し訳程度に本学関連の話でもしておこうかと思いましたが、そうですね…特に話すことが思いつきませんね。強いて言えば、純粋な医療系を諦めてしまった方は僕の所属する研究室に来るととても楽しいと思います。医療系を取り扱っている研究室は割とありますが、どちらかと言えば情報技術寄りなので、ヒトの身体の仕組みとかそっちの方ではありません。しかし、うちの研究室(というより僕の所属する研究領域)は生理学や神経科学、運動学などの基本的な領域から、片麻痺を中心に理学療法作業療法、種々の評価学など、普通に医療系の専門学校や大学で学ぶことを、全てを網羅しきれてはいませんがある程度は勉強しています*7。加えて、いくつかの医療機関と共同研究を行っており、日々勉強をさせて頂いています。


医療系をやりたかったけど本学は情報系だからもうだめだ…と思っている皆さま。ぜひ本研究室にいらしてください(どんな本学関連の話をしようかと2分ほど考えた結果、研究室勧誘に落ち着きました)。一緒に楽しい事をしましょう。


というわけで、僕の記事はこれで終わりにしたいと思います。良いクリスマスをお過ごしください!それでは!


参考文献

赤穂昭太郎,カーネル多変量解析(非線形データ解析の新しい展開),岩波書店,2008.

ソースコード

ロジスティック回帰をカーネル化したクソ雑魚プログラムを貼っておきます。高速化がいかに大事かを思い知りました。ちなみに、今回は識別精度を求める事が目的ではない上に、クソ雑魚プログラムによってチューニングもしていない(遅すぎて何度も回せなかった…つらい)ので、train/val/testには分けていません。

# seibisi

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
from sklearn.preprocessing import StandardScaler
from matplotlib.colors import ListedColormap
from sklearn.linear_model import LogisticRegression


class kernelization:
    def __init__(self):
        pass
    # ガウシアンカーネルを返す
    def gauss_kernel(self, x1, x2, gamma):
        return np.exp(gamma * np.linalg.norm(x1-x2) ** 2)    
    # グラム行列を作成する
    def gram_matrix(self, X, gamma):
        return [[self.gauss_kernel(x1, x2, gamma) for x1 in X] for x2 in X]
    
class pl:
    def __init__(self):
        pass
    
    def main(self, X, Y, lr, gamma):
        markers = ("x", "s")
        cmap = ListedColormap(("red", "blue"))
        X = X.values
        # メッシュの作成
        x1_min, x1_max = X[:, 0].min()-1, X[:, 0].max()+1
        x2_min, x2_max = X[:, 1].min()-1, X[:, 1].max()+1
        x1_mesh, x2_mesh = np.meshgrid(np.arange(x1_min, x1_max, 0.1), np.arange(x2_min, x2_max, 0.1))
        # メッシュの分類
        k = kernelization()
        tst = np.array([x1_mesh.ravel(), x2_mesh.ravel()]).T
        ker = [[k.gauss_kernel(np.array(t), np.array(xx), gamma) for xx in X] for t in tst]
        z = np.array([lr.predict(pd.DataFrame([kk])) for kk in ker])
        z = z.reshape(x1_mesh.shape)
        # プロット
        plt.contourf(x1_mesh, x2_mesh, z, alpha=0.4, cmap=cmap)
        plt.xlim(x1_mesh.min(), x1_mesh.max())
        plt.ylim(x2_mesh.min(), x2_mesh.max())
        l = ["C1","C2"]
        for idx, cl in enumerate(np.unique(Y)):
            plt.scatter(x=X[Y == cl, 0],
                        y=X[Y == cl, 1],
                        alpha=0.6,
                        c=cmap(idx),
                        edgecolors='black',
                        marker=markers[idx],
                        label=l[idx])
        plt.xlabel("X1")
        plt.ylabel("X2")
        plt.legend()
        

def main():
    """
    # 線形分離不能
    p1 = pd.DataFrame(np.random.multivariate_normal([10, 10], [[5, 0], [0, 5]], 50), columns=["X1", "X2"])
    p2 = pd.DataFrame(np.random.multivariate_normal([30, 30], [[5, 0], [0, 5]], 50), columns=["X1", "X2"])
    p = pd.concat([p1, p2], axis=0).reset_index(drop=True)
    n1 = pd.DataFrame(np.random.multivariate_normal([10, 30], [[5, 0], [0, 5]], 50), columns=["X1", "X2"])
    n2 = pd.DataFrame(np.random.multivariate_normal([30, 10], [[5, 0], [0, 5]], 50), columns=["X1", "X2"])
    n = pd.concat([n1, n2], axis=0).reset_index(drop=True)
    yp = pd.Series(np.zeros(len(p)))
    yn = pd.Series(np.ones(len(n)))
    df = pd.concat( [pd.concat([p, yp], axis=1), pd.concat([n, yn], axis=1)], axis=0)
    df = df.rename(columns={0:"Y"}).reset_index(drop=True)
    df.to_csv("df.csv", index=False)
    """ 
    
    # 1. 前処理
    df = pd.read_csv("df.csv")
    X = df.drop(columns="Y")
    Y = df.Y
    sc = StandardScaler()
    X = pd.DataFrame(sc.fit_transform(X), columns=X.columns)
    # 2. カーネル化
    gamma = 1.0
    k = kernelization()
    K = pd.DataFrame(k.gram_matrix(X.values, gamma))
    # 3. カーネル化したらロジスティック回帰
    lr = LogisticRegression(C=1.0)
    lr.fit(K, Y)
    # 4. 可視化
    p = pl()
    p.main(X, Y, lr, gamma)
  

*1:最近は知りませんが過去に後輩からそんな雰囲気の事を言われました笑

*2:専門の方から見たら根本的な間違い等を犯しているかもしれません。その時は僕の不勉強です、申し訳ありません。

*3:さっさと他の識別器を使えよって話ですが

*4:重回帰モデルで重みを計算する時にL2正則化という処理を行った回帰です

*5:プログラミングは高専(しかも情報工学科)の1年生の時に諦めました

*6:実際のところどうなんですかね?研究で手が回らなくて検証する気が起きていないのですが…笑

*7:というよりもそれを知らなければ話にならないので

2標本の検定・尺度母数の推定(ノンパラメトリック検定)

以下の記事で作成したプログラムです。グラフの部分は省略しました。
seibibibi.hatenadiary.jp

# seibisi
sampling <- function(pi, mus, sig, N){
	# 一様乱数をサンプル
	u <- runif(N, 0, 1)
	data <- numeric(N)
	dens <- numeric(N)
	for (i in 1:N){
		for (j in 1:length(pi)){
			if (u[i] < pi[j]){
				break
			}
		}
		data[i] <- rnorm(1, mus[j], sig[j])
	}
	return (data)
}

# 0. いじるパラメータ
pi <- c(0.3, 0.7)
mus1 <- c(10, 20)
mus2 <- c(10, 20)
sig1 <- c(sqrt(5), sqrt(5))
sig2 <- c(sqrt(10), sqrt(10))
m <- 100
n <- 50
N <- 10000
# 1. データをサンプルして検定統計量の計算
MO1 <- data.frame()
MO2 <- data.frame()
F <- data.frame()
for (i in 1:N){
	data1 <- sampling(pi, mus1, sig1, m)
	data2 <- sampling(pi, mus2, sig2, n)
	MO1 <- rbind(MO1, mood.test(data1, data2)$statistic)
	MO2 <- rbind(MO2, mood.test(data2, data1)$statistic)
	F <- rbind(F, var.test(data1, data2)$statistic)
}
MO <- data.frame(MO1, MO2)
colnames(MO) <- c("MO1","MO2")
colnames(F) <- "F"
data <- data.frame(data1, data2)
# 2. 棄却域に落ちた個数
p_u <- qf(df1=100-1, df2=50-1, p=0.05, lower.tail=FALSE)
F_num <- length(F[F$F > p_u,])
p_l <- qnorm(p=0.025, mean=0, sd=1, lower.tail=TRUE)
p_u <- qnorm(p=0.025, mean=0, sd=1, lower.tail=FALSE)
MO_num <- length(MO$MO1[MO$MO1 < p_l]) + length(MO$MO1[MO$MO1 > p_u])
F_num
MO_num

【各種検定の理論】2標本の検定・尺度母数の推定(ノンパラメトリック検定)

導入

本記事では、ノンパラメトリック検定における2標本の検定(尺度母数の推定)について取り扱いたいと思います。尺度母数はデータの広がりを示すパラメータであり、分散のようなものです。後でも述べていきますが、今回取り扱う検定はMoodの検定と呼ばれており、尺度母数を検定するために2群の中央値が大体同じ位置にある事を前提としています。ここで1つの疑問が生まれます。中央値をずらせばMoodの検定は一体どういう振る舞いをするのか。早速遊んでいきましょう。

(※追記:シミュレーションの部分での検定統計量の山が2つ現れる謎については単に僕の勉強不足というだけでした。詳細は以下の記事の追記をご覧ください(2019年12月10日 追記))
seibibibi.hatenadiary.jp



1.Moodの検定の概要

しくみ

Moodの検定の仕組みは以下に示す図です。

f:id:seibibibi:20191109110823p:plain

中央値が大体同じな2つの母集団分布を2パターン用意しました。左は大体同じ形の分布関数が2つ存在していますが、右は少しだけ様子が違います。具体的にどう違うかというと、 \displaystyle F_2(x) \displaystyle F_1(x)を横に引き伸ばした感じになっています。この状況だとどんなことが起きるでしょうか。横に引き伸ばされると、より小さな値やより大きな値がサンプルされることが予想できます。つまり、横に引き伸ばされた母集団分布からは、中央値から遠いデータが観測されるのではないでしょうか。以上の事から、2つの標本があり、中央値からの距離がどちらも同じ位なら同じようなばらつきがあると言えそうです。しかし、横に引き伸ばされた例の様に、どちらかの標本の中に中央値から遠い距離が割と測定されているなら2つの標本は同じようなばらつきを持っているとは言えないのではないでしょうか。これを順位を用いていい感じにしたのがMoodの検定になります

仮説

仮説は以下のようになります。

 \displaystyle H_0:F_1(x)=F_2(x)

 \displaystyle H_1:F_2(x)=F_1(\Delta x)

帰無仮説が何だか謎のように見えますが、イメージ的にはしくみで述べた通りです。また、 \displaystyle \Delta xはxの微小変化のように見えますが、 \displaystyle \Delta (\Delta>1)は尺度母数でただの数字と同じ扱いである事に注意です。引き伸ばし度合いのようなものです。これを踏まえると、帰無仮説は2つの母集団分布は大体同じようなばらつき(引きのばしの程度)*1であり、対立仮説は2つの母集団は同じばらつきは持っていないという事になります。

検定統計量

検定統計量の考え方を図で示します。以下の図を見てみましょう。

f:id:seibibibi:20191109205240p:plain

2標本の検定・位置母数の推定と同じように、まずは全てのデータを1つにまとめて昇順にした後で、群1に属するデータに1、群2に属するデータに0を割り当てます。その後、平均順位を求め、平均順位と1を割り当てたデータとの差の2乗を累積します。操作だけを見ると、普通に分散を求めているような気がしますよね。検定統計量の式は以下になります。
 \displaystyle MO = \sum_{i=1}^{m} (R_i-\frac{m+n+1}{2})^2

ここで、 \displaystyle R_iは群1に属する \displaystyle i番目のデータの順位、 \displaystyle \frac{m+n+1}{2}は平均順位になります。ちなみに、Moodの検定では差の二乗を累積させていますが、ここを絶対値を累積させると、アンサリー・ブラッドレー検定という物が現れます。詳細は省略します。

2.シミュレーション

やること

それでは実験に入っていきましょう。今回検証する事は導入でも述べた通り、中央値をずらすとMoodの検定はどのような振る舞いをするのかです。行う実験の手順の図を以下に示します。

f:id:seibibibi:20191109224852p:plain

中央値をずらしたり一致させたりした2つの母集団分布から、赤の分布から100個のデータ、青の分布から50個のデータをサンプリングします。サンプリングした後は検定統計量を求め、ヒストグラムを作ります。標本抽出は全体で10000回行う事にします。加えて、今回用いる母集団分布の例を以下の図に示します。

f:id:seibibibi:20191109221151p:plain

位置母数の推定と同じように混合正規分布からデータを発生させます。この2つの混合正規分布について、分散を大きくしたり小さくしたり、左に動かしたり右に動かしたりして検定統計量がどのように推移するのかを検証します。

検定統計量はRに存在するmood.testを利用したいと思います。以下にdocumentationのリンクを貼っておきます。また、mood.testで計算される検定統計量は標準化されています。Moodの検定もまた正規近似できる検定統計量である事が知られています。
mood.test function | R Documentation
ちなみにですが…

There are more useful tests for this problem.

という何やら不穏な1文が書いてありますが、気にしない事にします。何を問題にしているのかはリンク先に飛べば分かりますが、何も分からないという事が分かります。ノンパラと言ってるのに中央値を利用するんじゃねえよって事なんですかね。真相は謎のままです*2

今回は(今回から)、僕がまだノンパラ検定の理論を深く知っていない*3事の反省を込めて、あまり深入りをしない事にしました。ノンパラ検定の理論をしっかり勉強したら、しれっと修正しておきたいと思います*4

実験結果

実験1 2つの母集団で分散と中央値を完全に一致させる

まずは中央値を完全に一致させるように混合正規分布を配置しました。2つの母集団分布は以下のような様子です。2つの混合正規分布の平均と分散は完全に一致させました。*5

f:id:seibibibi:20191109221533p:plain

この時、mood.test(data1, data2)とmood.test(data2, data1)と2パターンの検定統計量を求めた時のヒストグラムを以下に示します。

f:id:seibibibi:20191109231926p:plain

全く同じ母集団分布を2つ用意したなら、検定量が棄却域に落ちないのは納得ですね。見事に2つとも一致しています。

実験2 中央値を一致させて分散を変える

今後は中央値を一致させて分散を変えてみましょう。これが本来のMoodの検定の前提条件となっています。用いる2つの母集団分布*6を以下に示します。

f:id:seibibibi:20191109223118p:plain

検定統計量のヒストグラムを以下に示します。

f:id:seibibibi:20191109233125p:plain

Mann-Whitneyの時と同じく2つに分かれましたね。ちなみに、山が二つに分かれてはいますがp値はどちらも同じ値を示します。この2つの山の平均を求めれば平均は0になります。

実験3 分散は同じにして中央値を変える

今度は分散を同じにして中央値を変えてみましょう。以下に2つの母集団分布*7を示します。

f:id:seibibibi:20191109223840p:plain

検定統計量はどうなるのでしょうか。以下にヒストグラムを示します。

f:id:seibibibi:20191109233328p:plain

両方とも棄却域に値が集中していますね。中央値が離れてしまうと正しく検定できないんですね。当たり前のことですが…

追加実験 F検定との比較

F検定はパラメトリック検定の枠組みの中で等分散かどうかを調べる検定です。Mood検定とF検定を比較してみましょう。実験条件は、実験2と同じデータセットを用います。

f:id:seibibibi:20191109234719p:plain

今回の実験では、F検定における検定統計量は自由度が(99, 49)のカイ二乗分布に従います。有意水準 \displaystyle \alpha0.05とした時、F分布では35個、Moodの検定では1656個が棄却域に落ちていました。やっぱり約束を守らないと変な結果になっちゃうんですね。

感想

やくそくをやぶるってたのしい

まとめ

今回はMoodの検定で遊びました。やはりMoodの検定も引数の渡し方(どちらの標本を基準にするか→1を割り当てるか)によって結果は変わるようですね。2標本のノンパラ検定の全てにこういう傾向が見られるのでしょう。この辺りはノンパラ検定の理論を深く追求した時に、まとめてどこかで書きたいと思います。

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

*1:イコールと書いてしまうと全く同じ関数であると勘違いしてしまいますよね…有意差が出なかった時の解釈に困るなあと個人的に思います。伸びてるか伸びてないかで言えば伸びてはいるんだろうけど何とも言えないよねえ…的な事を言える仮説があればなあ…

*2:何を問題にしているのかが原文には明記されていない

*3:2019年11月9日現在

*4:とりあえず現在は医療統計が入り用になっているので

*5:両方とも \displaystyle N(10, 5), N(20, 5)を0.3と0.7の割合で混合させました。

*6: \displaystyle N(10, 5), N(20, 5)を0.3と0.7の割合で混合したものと、 \displaystyle N(10, 10), N(20, 10)を0.3と0.7の割合で混合したものです。

*7: \displaystyle N(10, 5), N(20, 5)を0.3と0.7で混合した物、 \displaystyle N(20, 5), N(30, 5)を0.3と0.7で混合した物です。

【各種検定の理論】2標本の検定・位置母数の推定(ノンパラメトリック検定)

導入

本記事では、ノンパラメトリック検定における2標本の検定について取り扱いたいと思います。パラメトリック検定で言う所の、対応なしt検定に相当します。今回はウィルコクソン順位和検定とマン・ホイットニーのU検定を紹介します。後で述べていきますが、ウィルコクソン順位和検定とマン・ホイットニーのU検定は実質的に同じことをしており、ウィルコクソン順位和検定の検定統計量からマン・ホイットニーの検定統計量に数理的に繋げることができます。ここで次のような疑問が生まれます。本当に何もかも完全に一致するのか。僕は初めて同じである事を知った時に本当か?と疑問に思いました。だって同じことをするなら、どっちかの存在意義が無くなりませんか?なぜ2つ残っているんでしょうかね?気になったら即手を動かして遊んでみましょう。


(※追記:本当はどちらもぴったり一致する事を確かめようと思ったのですが、マン・ホイットニーのU検定が意味分からんなので、主にマン・ホイットニーのU検定の話をしていきます…勉強不足ですみません…とりあえずやろうとしていたことはメモとして残しておきます。できなかったところは線分を引っ張っておきます。)

(※追記2:後で示しますが、期待値の謎問題については現在進行形で検証中です。深い理論の話になってしまうので医療従事者のための統計学には書きませんが、また改めて数理統計学の枠を設けて遊んでいきたいと思います。理論が完全に明らかになり次第、しれっと修正しておきます(2019年11月9日 追記))

(※追記3:期待値の謎問題は解決しました。帰無仮説 \displaystyle F_1(x)=F_2(x)の下での期待値である事をすっかり忘れていました。従って、2つ山ができるのは、単にwilcox.test(data1, data2)とwilcox.test(data2, data1)で順番を変えただけであってそりゃそうという結論になりました。ちなみに、このブログで述べているMann-WhitneyやMoodの検定などは、帰無仮説の下での検定統計量の分布は左右対称になっています。従って、例えばwilcox.test(data1, data2)とwilcox.test(data2, data1)で算出される検定統計量を足してみると期待値の2倍になる事が確認できます。不勉強のために一人で騒いでしまいました…とりあえず勉強の後を残すために、一人で騒いでいる様子はそのまま残しておきたいと思います(2019年12月10日 追記))

1.ウィルコクソン順位和検定の概要

しくみ

1標本の時と同じように、初めに仕組みの話をします。まずは簡単に前提条件と文字の定義をしておきましょう。
 \displaystyle X_1 = (X_{11}, ..., X_{1m})を分布関数 \displaystyle F_1(x)からの確率標本、 \displaystyle X_2 = (X_{21}, ..., X_{2n})を分布関数 \displaystyle F_2(x)からの確率標本とします。ここで、分布関数 \displaystyle F_1(x), F_2(x)は連続である事以外には何も分からないものとします。次の図を見てみましょう。

f:id:seibibibi:20191107191020p:plain

母集団の分布関数は点線で、データは丸で表しています。データを見て、左の状態と右の状態で、どちらがデータに差が表れているように見えるでしょうか。例えば、この2群のデータに対して順位をつけてみましょう。左側の状態では各群で順位がばらばらであり、右の状態では全体的に青の分布の方が順位が高い事が分かるのではないでしょうか。これが本記事で述べる検定の仕組みです。もし、2群間でデータに差が表れているなら、どちらかの群のデータはもう一方の群のデータよりも順位が高くなるはずです。では、この仕組みをどうやって実現しているのでしょうか。

仮説と検定統計量

仮説

まず初めに仮説を示します。

 \displaystyle H_0:F_1(x)=F_2(x)

 \displaystyle H_1:F_2(x)=F1(x-\Delta)

ここで、謎の対立仮説が現れました。この仮説を図で表すなら、以下のようになります。
f:id:seibibibi:20191107191231p:plain

 \displaystyle \Deltaは移動量のようなものであり、位置母数と呼ばれています。つまり、一方の分布を \displaystyle \Deltaだけ移動させると完全に重なります。この事から、ウィルコクソン順位和検定は、データに差が表れているなら分布関数が大きくシフトしているはず(差が大きい)という事を検証します。

検定統計量

次に、検定統計量を述べていきます。しくみの所でも述べたように、基本的には全てのデータを1つにまとめた後で順位で並べ替えて、各群の順位を見るという処理を行います。式を書く前にイメージを図で示します。

f:id:seibibibi:20191107200523p:plain

ごくシンプルに、順位と割り当てた数の積和になっています。この事を踏まえて検定統計量を眺めてみましょう。検定統計量は次の式で表されます。
 \displaystyle W = \sum_{i=1}^{m}R_i

ここで、 \displaystyle R_iは群1に所属するデータの順位です。0を割り当てた物は総和しても仕方が無いので、1を割り当てた物だけを総和しています。とてもシンプルですよね。

2.マン・ホイットニーのU検定

ウィルコクソン順位和検定との違い

基本的に、ウィルコクソン順位和検定と同じ仕組み、同じ仮説を設定しています。従って、特にこれと言って説明する事は無いのですが、強いて言うなら検定統計量の導出が少し異なっています。

検定統計量

初めに図を示します。

f:id:seibibibi:20191107200412p:plain

マン・ホイットニーのU検定の場合は、網羅的にデータの大きさを比較します。群1のデータの方が大きければ1を割り当て、そうでなければ0を割り当てます。最終的に割り当てた数字を全て総和してあげれば、マン・ホイットニーのU検定の検定統計量が導出されます。式は以下の通りです。
 \displaystyle MW = \sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)

ここで、 \displaystyle U(i,j)は、 \displaystyle X_{1i} \geq X_{2j}の時に1、それ以外が0となる量です。要するに上で示した処理を行っています。

何が同じ

マン・ホイットニーのU検定の検定統計量はウィルコクソン順位和検定の検定統計量から導出できます。具体的には、

 \displaystyle W=\sum_{i=1}^{m}R_i = \sum_{i=1}^{m} (\sum_{j=1}^nU(i,j)+i) = \sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}

イメージ的には、ぜんぶ1つにまとめて並べた時に、母集団1の任意のサンプルから次のサンプルまでに入っている母集団2のサンプルの個数を足し合わせている感じでしょうか。図で言えば以下のようになります。


f:id:seibibibi:20191108113648p:plain


やってることが同じな気がしてきましたね、なんだかとてもつまらない予感がしてきました…

3.シミュレーション

やること

さて、ここからシミュレーションに入ります。導入でも述べた通り、確かめたい事はウィルコクソン順位和検定とマン・ホイットニーのU検定は完全に一致するのかのはずでした。以下にやろうとしていた実験の手順を示します。

f:id:seibibibi:20191107202142p:plain

母集団分布を2つ用意し、各群でデータをサンプルします。サンプルした後で各検定の検定統計量を求めます。この手順を複数回繰り返し、検定統計量に関するヒストグラムを作成します。ここで、今回完全に一致するのかどうかは検定統計量を正規近似した下で行う事とします*1。どちらの検定統計量も、サンプル数が大きい時には正規近似できることが知られています。従って、この2つの検定統計量を標準化し、標準正規分布上で比較したいと思います。もし本当に何もかも一致するなら、ヒストグラムは完全に一致するはずです。

ここで、用いる母集団分布を示しておきます。今までは分布関数について話をしてきましたが、以下で示す分布は確率密度関数である事に注意です。

f:id:seibibibi:20191107203845p:plain

2つの分布*2から、サンプルするデータの個数は赤の分布から100個、青の分布から50個、標本抽出は10000回行います。

正規近似のシミュレーション

理論的背景

1標本のノンパラ検定でも述べましたが、この検定においても正規近似を適用する事が出来ます。ある検定の検定統計量を \displaystyle Uとした時、標準化後の検定統計量 \displaystyle Zは次の式で正規近似できます。

 \displaystyle Z = \frac{U-E[U]}{\sqrt{V[U]}}~N(0, 1)

ここで、それぞれの検定量の期待値と分散は次の式で示されます。
 \displaystyle E[W]=\frac{m(m+n+1)}{2}, V[W]=\frac{mn(m+n+1)}{12}

 \displaystyle E[MW]=\frac{mn}{2}, V[MW]=\frac{mn(m+n+1)}{12}

という事で、この期待値と分散により正規近似していきます。

実験結果1 自作した検定統計量とパッケージとの対応

「R Mann Whitney」で調べてもwilcox.testを使っている物ばかりだし、「R Mann Whitney documentation」で調べても意味不明な物が出てきて萎えたので、それなら一から作ってしまおうという事で作りました。以下に作成したコードの1部分を載せます。

# 検定統計量の計算
calc <- function(data1, data2){
	# 両方向の比較を書いてみる(data1→data2とdata2→data1)
	# 1. data1→data2
	# Mann-Whitney
	MW1 <- 0
	d1 <- sort(data1)
	d2 <- sort(data2)
	for (i in 1:length(d1)){
		MW1 = MW1 + length(d2[d2 < d1[i]])
	}
	# たしかめ
	wilcox.test(data1, data2)
	# 2. data2→data1
	MW2 <- 0
	for (i in 1:length(d2)){
		MW2 = MW2 + length(d1[d1 < d2[i]])
	}
	MW2
	# たしかめ
	wilcox.test(data2, data1)
	return (c(MW1, MW2))
}

自分で作成したは良いものの、合ってるかどうかが肝心なので、たしかめの部分を実行するとこんな感じになりました。

> # Mann-Whitney
> MW1 <- 0
> d1 <- sort(data1)
> d2 <- sort(data2)
> for (i in 1:length(d1)){
+ MW1 = MW1 + length(d2[d2 < d1[i]])
+ }
> # たしかめ
> wilcox.test(data1, data2)

        Wilcoxon rank sum test with continuity correction

data:  data1 and data2
W = 1199, p-value = 2.163e-07
alternative hypothesis: true location shift is not equal to 0

> MW1
[1] 1199
> 

良い感じに一致していました。で、作った後で気付いたのですが、wilcox.testのdocumentationに次のような記述がありました。

R's value can also be computed as the number of all pairs (x[i], y[j]) for which y[j] is not greater than x[i], the most common definition of the Mann-Whitney test.

wilcox.test function | R Documentation
おう、まじか…ちゃんと調べれば良かった…と思いました(萎えないでちゃんと調べよう)。でも、逆に言えばウィルコクソンの検定統計量が計算されている訳では無いという事です。マン・ホイットニーの検定統計量とウィルコクソンの検定統計量の関係式を用いて計算したかったのですが、次に示す謎により計算できませんでした。

実験結果2 謎

とりあえずマン・ホイットニーのU検定の検定統計量のヒストグラムを書いてみます。以下の様になりました。

f:id:seibibibi:20191108160552p:plain

赤が群1→群2で比較した検定統計量、青が群2→群1で比較した検定統計量です。上で示したマン・ホイットニーの検定統計量の期待値の式に \displaystyle m=100, n=50を入れると、2500になります。しかし、どちらの山も2500にはならず、なんなら2つの山の真ん中が2500なのかなという感じがします。期待値が山のピークになる事は必ずしもそうではない(指数分布とか連続一様分布と考えると)ので、2つの山の真ん中が期待値と言われても別に違和感がない(むしろ期待値の式に突っ込めば真ん中になるのでしょう)と言えばないですが、うーん…わからない。両方考えなきゃなのかなあ、定義は片方だけなのに。こんな結果なので、ウィルコクソンの検定統計量を計算しようにも、 \displaystyle nを入れるのか \displaystyle mを入れるのか分からないので、残念ながら比較はできませんでした。

考察

つらい

まとめ

今回は、ウィルコクソン順位和検定とマン・ホイットニーのU検定で遊びました。本当にどっちも同じなのか?という問いに対しては、これも後から気付いたのですが、以下の式で証明できそうです。
 \displaystyle E[W] = E[\sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}] = E[MW+\frac{ m(m+1) }{2}]
 \displaystyle = \frac{mn}{2}+\frac{ m(m+1) }{2} = \frac{m(m+n+1)}{2} = E[MW]
 \displaystyle V[W] = V[\sum_{i=1}^{m}\sum_{j=1}^{n}U(i,j)+\frac{ m(m+1) }{2}] = V[MW]
同じじゃんか…

後から気付いて証明っぽいのもできてしまった上に、検定統計量の振る舞いも良く分からないし……なんもわからんです…がんばりたい…

*1:サンプル数が少ない時の帰無分布を求めるのがとても面倒くさいので、サンプル数が少ない時の検証は魔が差した時にします。

*2:赤、青ともに混合比が0.3と0.7の混合正規分布です。赤は \displaystyle N(10, 5), N(20, 5)の2つ、青は \displaystyle N(15, 5), N(50, 5)の2つを混ぜています。

1標本の検定(ノンパラメトリック検定)のプログラム

以下の記事で作成したプログラムです。
グラフの描画のソースは省略します。
seibibibi.hatenadiary.jp

# 混合正規分布からデータのサンプリング、検定統計量の取得
sampling <- function(pi, mus, sig, N){
	# 一様乱数をサンプル
	u <- runif(N, 0, 1)
	data <- numeric(N)
	dens <- numeric(N)
	for (i in 1:N){
		for (j in 1:length(pi)){
			if (u[i] < pi[j]){
				break
			}
		}
		data[i] <- rnorm(1, mus[j], sig[j])
	}
	return(data)
}


# 0. サンプリングに関するパラメータ
pi <- c(0.2, 0.8)			# 混合比(中央値は4.3)
#pi <- c(0.5, 0.5)			# 混合比(中央値は0)
mus <- c(-5, 5)
sig <- c(sqrt(5), sqrt(5))
n = 100	# 混合正規分布からの乱数の数
N = 10000	# シミュレーション回数

# 1. データをサンプルして平均、分散、検定統計量を取得
U <- data.frame()
for (i in 1:N){
	data <- sampling(pi, mus, sig, n)
	# 検定統計量の計算
	# wilcox.testは帰無仮説はΔ=0,符号検定の帰無仮説はΔ(中央値)=0
	w_n <- wilcox.test(data)$statistic
	s_n <- SIGN.test(data, md=0)$statistic	# BSDAパッケージのダウンロードが必要
	U <- rbind(U, c(w_n, s_n))
}
colnames(U) <- c("w_n","s_n")
# 標準化
U2 <- U
exp_w_n <- n * (n+1) / 4
var_w_n <- n * (n+1) * (2*n+1) / 24
exp_s_n <- n / 2
var_s_n <- n / 4
U2$w_n <- (U$w_n-exp_w_n) / sqrt(var_w_n)
U2$s_n <- (U$s_n-exp_s_n) / sqrt(var_s_n)

# 描画
graph(U)
graph2(U2)

【各種検定の理論】1標本の検定(ノンパラメトリック検定)

導入

本記事では、ノンパラメトリック検定における1標本の検定について扱っていきます。パラメトリック検定で言えば、1標本のt検定だったり対応のある2標本のt検定に相当します。ここではウィルコクソン符号付順位検定符号検定の2つを紹介しますが、それぞれの理論で前提条件が存在しています。後でも細かく述べますが、データを発生させる母集団となる分布関数に対称性が仮定されるならウィルコクソン符号付順位検定対称性が保証されなければ符号検定という物です。ここで2つの疑問が生じます。1つ目は、対称性が仮定されている所で対称性を仮定できない検定を用いるとどうなるのか。もう1つは、対称性を仮定できない所で対称性を仮定した検定を利用するとどうなるのか。後から述べていきますが、この話は中央値が鍵となります。対称な分布からデータが発生すれば、対称の中心と中央値は一致しますが、そうでない場合はそもそも対称という考え方が存在していないので、中央値で考えるしかありません。で、この中央値は一体どんな重要な役割を果たしていくのでしょうか。早速遊んでいきます。

1.ウィルコクソン符号付順位検定の概要

仕組み

初めに仕組みをしっかり押さえると、後々が入ってきやすい気がするので、先に結論となる仕組みの部分を述べます。以下の図を見てみましょう。

f:id:seibibibi:20191106141133p:plain

未知の分布関数が点線で表され、この分布に従ってデータがx軸上の丸としてサンプリングされたとします。この状態から母数である \Deltaはどの辺りかな?という推定を行う事がこの検定の本質になります。分布関数に対称性を仮定している場合、理想的にはデータは \displaystyle \Deltaを中心に左右対称にサンプル*1されます。この性質を踏まえて、多分 \displaystyle \Deltaはここだ!と \displaystyle \Delta = \Delta_0と予想した時に、次の2つのパターンが考えられます。

  • 予想した値 \displaystyle \Delta_0と母数 \displaystyle \Deltaかなり近い値を取る

データの出方は \displaystyle \Delta_0を中心に左右対称に見えるはずです。

  • 予想した値 \displaystyle \Delta_0と母数 \displaystyle \Delta大きくかけ離れている

データの出方は \displaystyle \Delta_0を中心に左右対称にはならないはずです。上の図が典型的な例になります。

 \displaystyle \Deltaはここだ!と予想した後の分布帰無仮説の下での分布)に対してそのデータは「「「図形的には」」」左右対称に表れているのかどうかを見るのがこの検定の仕組みです(厳密には違う気がするので注釈*2に回します)。図で言えばという言葉を強調したのは、厳密な説明ではない気がしたからです。

前提条件

仕組みの所で述べてしまいましたが、ウィルコクソン符号付順位検定では、母集団の分布関数が①連続であり②対称性を有する事が前提条件となります。

仮説

帰無仮説と対立仮説は以下のようになります。

 \displaystyle H_0:\Delta = \Delta_0

 \displaystyle H_0:\Delta \neq \Delta_0

上で示した対立仮説は両側仮説ですが、当然片側仮説も存在しています。

検定統計量

検定統計量の考え方を図に示します。

f:id:seibibibi:20191106141525p:plain

 \displaystyle Z_j = |X_j-\Delta_0| (j=1..n)とする時、検定統計量は次の式で定義されます。
 \displaystyle W_n = \sum_{i=1}^{n}iV_i

ここで、 \displaystyle V_i \displaystyle | Z_{(i)} |における符号が正の時に1、負の時に0になる量です。もし、帰無分布が母集団分布から大きく外れてしまった場合は、 \displaystyle W_nが極端に大きくなる( \displaystyle \Delta_0を小さく見積もってしまった=対立仮説( \displaystyle H_1:\Delta>\Delta_0)が成立)極端に小さくなるか( \displaystyle \Delta_0を大きく見積もってしまった=対立仮説( \displaystyle H_1:\Delta<\Delta_0)が成立)のどちらかになります。後は棄却域*3を設けて、基準となる値よりも大きければ(小さければ)帰無仮説を棄却します。これがウィルコクソンの符号付順位検定です。

2.符号検定の概要

ウィルコクソン符号付順位検定との違い

ウィルコクソン符号付順位検定では、分布関数に対称性を仮定していましたが、対称性が保証されない場合はどのような問題が生じるでしょうか。以下の図を見てみましょう。

f:id:seibibibi:20191106142745p:plain

対称でない母集団の分布関数に対して、対称である分布を仮定すると、帰無仮説の下での分布においてサンプルに偏りが生じます。上で示した図は母集団の分布関数を図で示しているため、対称じゃねえ!とかそこは中央じゃねえ!とか突っ込めますが、データだけを見ると意外に分からない物です。以下の図を見てみましょう。
f:id:seibibibi:20191106144520p:plain

いかがでしょうか?あんまり良く分からないと思います。*4

符号検定の考え方

先に仮説を示します。

 \displaystyle H_0:\Delta = \Delta_0または \displaystyle Pr(X>\Delta_0)=Pr(X<\Delta_0)=\frac{1}{2}

 \displaystyle H_1:\Delta \neq \Delta_0または \displaystyle H_1:Pr(X>\Delta_0) \neq \frac{1}{2}

ここで、 \displaystyle \Delta中央値であり、上で示した対立仮説は両側仮説ですが、片側仮説の場合は不等式になります。

なんだか難しく見えますが、符号検定では中央値を基準にしてデータの現れる確率が同じかどうかについて検討します。以下の図を見てみましょう。

f:id:seibibibi:20191106150732p:plain

実は、符号検定の初めの所で紹介した図における母集団の分布関数は、こちらの確率密度関数から生成した物でした*5。具体的には、この確率密度関数から10000点のデータをサンプルし、累積させることで分布関数を作成しました。分布関数を作成するのと同時に中央値を求めてみた所、14.27という値が出てきました。図で言えば赤いライン付近に当たる部分です。
中央値を基準にしてデータの現れる確率が同じかどうかという事は、上の確率密度関数で説明するなら*6、赤いラインより左側の面積が0.5、赤いラインより右側の面積も0.5という事です。例えば、予想した中央値が母集団の中央値に対して小さく見積もってしまった場合、予想した分布からデータを見ると、やたらと中央値より大きいデータが出てくるはずです。これを確かめるのが符号検定の考え方になります。

検定統計量

検定統計量は以下の通りです。

 \displaystyle S_n = \sum_{i=1}^{n} I(X_i)

ここで、 \displaystyle I(x) \displaystyle x>\Delta_0の時に1であり、それ以外は0となる量です。つまり、中央値よりも大きくなるデータの個数を数え上げたのが符号検定の検定統計量*7です。先ほどの予想した中央値が母集団の中央値に対して小さく見積もってしまった例の場合、この検定統計量は非常に大きな値を取る事が考えられます。

3.シミュレーション

やること

さて、説明が長くなってしまいましたが、本題であるシミュレーションの話に入りましょう。導入でも述べた通り、今回行う実験は、対称性が仮定されている所で対称性を仮定できない検定を用いるとどうなるのか(実験1)を確かめる事と、対称性を仮定できない所で対称性を仮定した検定を利用するとどうなるのか(実験2)を確かめる事でした。今回実験で用いる2つの分布を以下に示します。
(これまでの話は分布関数の話をしてきましたが、以下の図は確率密度関数です。違いに注意してください)

f:id:seibibibi:20191107143048p:plain

左側の確率密度関数は対称であり、サンプルはどちらの山からも均等にサンプルされます。しかし、右側の確率分布は非対称であり、サンプルは右側の山から主にサンプルされます。つまり、右側の確率分布は対称の中心と中央値が一致していないケースです。この状況で、まあ対称の中心は0だろうどうせと適当に検定*8してみたらどうなるのかという事を実験してみます。

具体的な実験の手順としては、この山から100個のデータをサンプルし検定統計量を求めます。サンプルして検定統計量を求めて、またサンプルして検定統計量を求めて…を各実験で10000回行い、検定統計量に関するヒストグラムを作成します。
また、今まで述べてきた検定統計量はサンプル数が大きい時、検定統計量から平均を引いて標準偏差で割った量は標準正規分布に近づく事に知られています。これを利用してパーセント点を求め、正しい仮定を行った検定統計量のヒストグラムと誤った過程を行った検定統計量のヒストグラムでどのような違いが見られるかを検討していきます。

実験1では、帰無仮説として0を置くべきを中央値を置いてしまった事を想定します。実験1で用いる分布は対称の中心と中央値がだいたい一致しているので、さほど結果は変わらないことが予想できます。
実験2では、帰無仮説として中央値を置くべきを0としてしまった事を想定します。実験2はそもそも対称という考え方が無いので、おかしい結果が出てくることは目に見えています。

2つの実験をするにあたり、予め母数である中央値は既知であるとします。具体的には、左側の確率密度関数の中央値は0、右側の確率密度関数の中央値は4.3としておきます。

実験0 検定統計量の正規近似

理論的背景

本題に入る前に、予備実験として検定統計量の正規近似の実験を行います。ウィルコクソン符号付順位検定と符号検定では、検定統計量が従う分布の導出が少しだけ異なります(詳細は注釈*9に回します)。棄却域が2つの検定で変わられては困るので、検定統計量の正規近似を利用します。サンプルサイズが大きい時には、帰無仮説の下での検定統計量 \displaystyle U中心極限定理を利用して次のように正規近似する事が出来ます。

 \displaystyle Z = \frac{U-E[U]}{\sqrt{V[U]}}~N(0, 1)

ウィルコクソン符号付順位検定の検定統計量 \displaystyle W_nと符号検定の検定統計量 \displaystyle S_nもサンプルサイズが大きい時は正規近似されます。ここで、それぞれの検定統計量の期待値と分散は次の式で表されることが知られています。
 \displaystyle E[W_n] = \frac{n(n+1)}{4}, V[W_n] = \frac{n(n+1)(2n+1)}{24}

 \displaystyle E[S_n] = \frac{n}{2}, V[S_n] = \frac{n}{4}

以上より、取ってきたデータに対して、①平均と分散を求めて理論値になるかどうか②標準化した後の検定統計量は標準正規分布に本当に従うのかの2点を確かめていきたいと思います。

ちなみに、「サンプルサイズが大きい時」というのは一体どこから?という話ですが、例えばRのwilcox.testのdocumentationには以下のような記述があります。

… an exact p-value is computed if the samples contain less than 50 finite values and there are no ties. Otherwise, a normal approximation is used.

wilcox.test function | R Documentation
どうやら50より大きければ正規近似が適用されるようです。符号検定については特に記述が無かったので、とりあえずはサンプル数を50で……とは言わずに、冒頭で述べたようにでっかく100くらいにしておきました。どうせシミュレーションですし。

実験結果

初めに、それぞれの検定統計量の平均と分散のヒストグラムを示します。実験条件は、サンプルサイズを100、標本抽出を10000回行いました。

f:id:seibibibi:20191107141332p:plain

赤のグラフが \displaystyle W_n、緑のグラフが \displaystyle S_nです。手に入れた検定統計量から計算した平均と分散はそれぞれ、 \displaystyle W_nが平均2521.41、分散が83157.5で、 \displaystyle S_nが平均49.94、分散が24.49でした。理論的には \displaystyle W_nの期待値は2525、分散が84587.5で、 \displaystyle S_nが平均50、分散が25なので、かなり理論値に近い値が出ている事が分かります。

もう上のヒストグラムから察しは付いているような気はしますが、それぞれの検定統計量を標準化してみましょう。以下の図がその結果になります。

f:id:seibibibi:20191107142407p:plain

標準正規分布になりました。中心極限定理おそるべし…

実験1 対称性を仮定すべきを非対称を仮定した場合

標準化後の検定統計量のヒストグラムを以下に示します。

f:id:seibibibi:20191107155231p:plain

対称の中心と中央値が一致しているので、2つの山はほとんど重なっている事が分かります。

実験2 非対称性を仮定すべきを対称性を仮定した場合

こちらも標準化後の検定量ヒストグラムを以下に示します。

f:id:seibibibi:20191107161359p:plain

中央値が4.3である所を、帰無仮説として0を仮定してしまったウィルコクソンはほとんどの検定量が棄却域に落ちている事が分かります。

考察

そりゃそう

結論

今回は、ウィルコクソン符号付順位検定と符号検定で遊んでみました。ここまで理論を勉強してプログラムを考えてみて、そもそも対称だ非対称だを考える事は無いんじゃないかなあと思いました。その上、どうやったら検定の適用違いを起こすことができるかなあと思って考えていたのですが、この2つの検定の適用違いの例は思いつきませんでした。なので、当たり前の実験しかできませんでした…残念…

と、ここまで偉そうに記事を書いてしまいましたが、とても自信がありません。というのも、ノンパラ検定の理屈の話を勉強して大して経っていない*10上に、手を動かしながら勉強したい、勉強ノートのようなものがあればモチベが上がるという動機で始めたからです…間違いがたくさんある気がします…今回は以上でした!


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

*1:対称な分布関数を微分すると対称な確率密度関数が出てくると思います。ノンパラメトリックの話をしている所で持ってくる分布ではありませんが、例えば正規分布とか一様分布とか

*2:厳密には違う気がしています。これは僕がノンパラ検定の理論の勉強をし始めたばかりという事で確証が無いので注釈に回しました。ウィルコクソン符号付順位検定は符号検定と比較すると、検定統計量に順位の情報を用いている所が特徴的です。つまり、全体に対するデータの現れ方のバランスを定量化しているのではなくて、順位で重みづけた情報を見ています。この事から、ウィルコクソン符号付順位検定は、順位の和がある一定の値を超えるかどうかを見ているに過ぎないような気がしています。一方で、符号検定は後で述べる検定統計量と帰無仮説から考察するに、中央値を基準として確率が左右対称かどうかを見ているのだと思います。なので、左右対称という言葉の定義がまず難しいのですが、確率的に左右対称と条件づけるなら、左右対称かどうかを見ているのは符号検定な気がしています。間違えていたらごめんなさい…

*3:考え得る全ての順位の組み合わせを考えた後でその順位が実現される個数を数え上げます。全ての個数に対して該当する順位が実現される個数の割合が、その順位の確率となります。全ての順位で個の確立を設ける事で帰無分布を作成する事が出来ます。

*4:値が小さいデータの所に若干寄せました。この状態で例えば対称性を有する分布を仮定すると、値が大きい所のデータが帰無仮説の下での分布では値が大きすぎる事になると思います。

*5:混合正規分布

*6:一応「上で作成した確率密度関数が母集団分布」と定義しておきます。結局は上で示した図も標本でしかないので、そこを議論してしまえば先に進まなくなってしまうので…

*7:よって、二項分布に従います。

*8:現実ではありえない話ですが、シミュレーションならあり得ない話も実現できます。

*9:ウィルコクソン符号付順位検定では、全ての順位の組み合わせから分布を導出します。一方で、符号検定の場合は検定統計量が二項分布に従うので、帰無仮説の下での分布は仮定せず、単に二項分布表から棄却域を求める事が出来ます。

*10:なんなら1週間経ったか経たないかくらい

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

導入

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

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

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