ノンパラメトリックな方法を用いた探索的データ分析

はじめに

千葉大学/Nospareの米倉です．今回はノンパラメトリックな方法（分析で用いる分布を仮定しない）を用いた，探索的データ分析について解説します．ノンパラメトリックな方法を用いることにより，モデルの特定化の失敗の影響を抑えたり，またパラメトリックな方法を用いる際の予備的分析になると思います．また簡単なA/Bテストにもいいのかなと思います．

Empirical Cumulative Distribution Function(ECDF)

Empirical Cumulative Distribution Function(ECDF)とは，累積分布関数を経験分布を用いて推定（データから推定）したものです．確率変数$X$が$x$以下になる確率$P(X\leq x)$を累積分布関数といいます．これをデータから推定したものがECDFです．推定方法は簡単で，例えば次のJuliaのコードではsort(data)でdataを小さい順に並べ替え，y=collect(1:1:n)/nで1/nからnまで等間隔で並ぶ数列を作り推定しています．

function ECDF(data)

    n = length(data)

    x = sort(data)

    y = collect(1:1:n)/n

    return(x,y)
end

ECDFはデータが従う分布を確認する際に有用な手法です．例えば多くの場合，データが何かしらの正規分布に従うことを想定することが多いと思います．この場合，データから標本平均と標本分散を求め，それらをパラメータに持つ正規分布から乱数を発生させECDFを計算し，それとデータのECDFを比較することによって，データが正規分布に従ってそうか否かを判断することができます．

例えば次の図は，平均１，分散１の正規分布から１０００個データを発生させそれを仮想的なデータとみなしてECDFを計算し，その仮想的データの標本平均・分散をもとめそれらをパラメータにもつ正規分布から１０００個乱数を発生させ，それぞれのECDFをプロットしたものです．dataはデータのECDF,Normalは推定した標本平均・分散をパラメータに持つ正規分布からサンプリングしたECDFを表します．明らかに両者は同じようなECDFをもつので，この場合データは正規分布していると判断していると良さそうです．

一方で次の図は同じことを自由度５のカイ二乗分布からデータを発生させた時に行ったものです．前回と違い，今回は両者で違うECDFをもつので，この場合はデータは正規分布していないと判断できます．

なぜヒストグラムではなくECDFなのか

データが従う分布の形を大まかに知りたいとき，ヒストグラムを使うことを考える方もいると思います．確かにヒストグラムも有用なノンパラメトリックな方法なのですが，弱点として用いる棒（bin）の数によって，異なった情報を得る点がります．棒の数次第では正規分布に見えたり，見えなかったりするので，このような問題がないECDFを用いて確認する方がベターだと思います．例えば次の図は自由度５のカイ二乗分布の例で用いたデータと，データから推定した標本平均・分散を持つ正規分布からサンプリングしたデータのヒストグラムを，棒の数を１０としてプロットしたものです．大げさな例ですが，この様にどちらともそんなに形状に大差がないので，データが正規分布していると判断してしまってもおかしくありません．

ブートストラップ法

例えばサンプルサイズが$n$のデータがあり，その平均を求めたいとします．$n$が十分大きい時は単純に標本平均を求めれば一致性があるのでよいですが，小さい場合，次のような方法が考えられます．
1重複を許して同じ確率でデータを再抽出（リサンプリング）して新しいサイズが$n$の標本を作る．
２求めた標本の標本平均を求める．
３これらを十分大きい$B$回繰り返し，$B$個の標本平均の標本平均を計算し，これを平均の推定値とする．

この手法を（平均を計算する）ブートストラップ法といい，典型的なノンパラメトリックな手法として広く使われています．ここでは標本平均を求める例を用いましたが，回帰係数や分散など基本的にはブートストラップは幅広く適用できます．

以下はJuliaによるブートストラップ法のサンプルコードです．関数の引数はdata=適用するデータ，func=適用する関数，size=ブートストラップ標本の個数です．

function Bootstrap(data,func,size)

    result = zeros(size)

    for i in 1:size
        bsample = sample(data,length(data))
        result[i] = func(bsample)
    end

    return(result)
end

ブートストラップ法を用いた仮説検定

ある統計量$T$が帰無仮説$H_0$が正しいと仮定した時に，統計量の観測値$t$よりも大きくなる確率$P(t\leq T\mid H_0)$を（片側検定の時）みなさん大好きp値と呼ぶのでした．p値の計算・仮説検定は色々とややこしかったり難しかったりするのですが，ブートストラップ法を用いると，帰無仮説・統計量に対応するp値の計算が比較的簡単に行えます．

今二つの分布，$F$と$G$からサンプリングされた，２つのデータセット$x\sim F$と$y\sim G$があるとして，それぞれのサンプルサイズを$n$と$m$と仮定します．また帰無仮説$H_0$は$F=G$とします．それぞれのサンプルサイズは同じでも違くても構いません．$z=(x,y)$と置き，考える統計量を$t()$とします．

まず帰無仮説が真の下で$z$に対して重複を許してリサンプリングを$n+m$回行い，最初の$n$個を$x_b$，残りの$m$個を$y_b$とし，$z_b=(x_b,y_b)$と置き，$t(z_b)$を計算します．これらを$B$回繰り返し，$B$個の$t(z_b)$を求めます．$v_b$を$t(z)\leq t(z_b)$の時$1$を取る変数とすると，対応するp値は

$$p_{bootstrap} = \frac{1}{B}\sum_{b=1}^B v_b$$
で計算できます．
以下がJuliaを用いた，２つのデータセットがあるときに，両者の平均が等しいか否かをブートストラップ法を用いて検定するコードです．

#理論値の計算
diff = mean(data1) - mean(data2)

#結合されたデータの平均の計算
data = append!(data1,data2)
totalave = mean(data)

#二つのデータの平均が同じようにする（帰無仮説が真の時のデータセット）
same_ave_data1 = data1 .- mean(data1) .+ totalave
same_ave_data2 = data2 .- mean(data2) .+ totalave


#ブートストラップ法を用いた帰無仮説が真の下での平均の推定値とその差
B_ave_data1 = Bootstrap(same_ave_data1,mean,100000)
B_ave_data2 = Bootstrap(same_ave_data2,mean,100000)
est_diff = (B_ave_data1.-B_ave_data2)

#p値の計算
p = sum(est_diff.>=diff)/length(est_diff)

例として，ECDFで用いた平均１，分散１の正規分布に従うデータと，自由度５のカイ二乗分布に従うデータに対して，解説したブートストラップを用いた平均の検定を行うと，値は0.00で有意に帰無仮説（2つのデータの平均は等しい）を棄却できます．

おわりに

株式会社Nospareでは統計学の様々な分野を専門とする研究者が所属しております．統計アドバイザリーやビジネスデータの分析につきましては株式会社Nospare までお問い合わせください．