More than 5 years have passed since last update.

色んな統計量のバラツキを求めたい！再標本化 (Bootstrap)のススメ

Last updated at 2017-02-19Posted at 2017-02-06

　データ解析はまだまだブームのようです。機械学習であれ何であれ、まずデータが与えられてすることは、データを視覚化して、解析のとっかかりとなる仮説をつかむことですね。仮説を検証するため、様々な解析をして様々な統計量を計算します。　

　そんなことをしていると、ある統計量の母集団におけるバラツキが知りたい、という事態に出くわすことがあります。

　それがデータの平均値であれば、標準誤差 (standard error of the mean)というものを計算すればOKです。標準偏差を、データ数の平方根で割ったものですね。

　ところが、解析では平均だけでなく、中央値 (median)、百分率 (centile)、相関係数 (correlation coefficient)など様々な統計量を見ていくことになります。そして、これらには標準誤差や信頼区間 (confidence interval)といったバラツキを定量化できる簡単な公式が存在しません。　

　そんなときに役に立つのが、再標本化 (resampling)、特に、ブートストラップ (Bootstrap) という手法です。

ブートストラップとは　

　ブートストラップではまず、再標本化 (resampling)によって元データから何千もの似たデータセットを作り出します。そして、それぞれのデータセットで、欲しい統計量を計算します。すると、その統計量の分布ができるので、その分布の標準偏差を計算することにより、その統計量の標準誤差を求めることができます。

　例として、20人のIQのデータがあったとします。値はそれぞれ、61, 88, 89, 89, 90, 92, 93, 94, 98, 98, 101, 102, 105, 108, 109, 113, 114, 115, 120, 138だとしましょう。分布は以下のような感じです。

　このとき、平均値 (mean) は100.85、平均値の標準誤差 (standard error of the mean) は3.45になります。

　ただ、こうした値の分布って、所得とかもそうですけど、平均値 (mean) よりも中央値 (median) の方が人の実感を反映してるんですよね。なので、今回はブートストラップを使って、中央値のバラツキ、すなわち中央値の標準誤差 (standard error of the median) を求めてみたいと思います。

ブートストラップのアルゴリズム

ブートストラップは、以下の流れで行います。

1: 何回再標本化するか決めます。
　500, 1000, 10000, 100000など、キリが良く、大きい数にしましょう。今回は1000でやってみます。

2: 再標本化して、人工的にデータを生成します。このとき、重複を許して元データの標本数だけ、データ標本を取り出します。

　例えば、元データが $X_{1}, X_{2}, . . ., X_{5}$ と $5$ 個あったとしたら、再標本化するたびに生成されるデータセットの例は以下のようになります。

再標本1 = $X 1, X 2, X 3, X 3, X 4$
再標本2 = $X 1, X 1, X 3, X 4, X 4$
再標本3 = $X 2, X 2, X 5, X 5, X 5$
......
再標本998 = $X 1, X 2, X 3, X 4, X 5$
再標本999 = $X 1, X 3, X 4, X 5, X 5$
再標本1000 = $X 3, X 3, X 3, X 3, X 4$

　こうして、ステップ１で決めた数だけ、重複を許しデータ標本を元データ数だけ取り出すことを繰り返し、人工的にデータを作ってしまいます。

　各再標本にどのデータ標本が入るかはランダムです。5個すべて同じ値が入るときもあるかもしれませんし、5個全て違う値が入るときもあるでしょう。

3: 各再標本で、求めたい統計量を計算する。
　中央値なりなんなり、それぞれの再標本データを使って統計量を求めます。1000回再標本化したら、統計量も1000個とれますね。

再標本1 = $X 1, X 2, X 3, X 3, X 4$ 　 $\to$ 統計量 = $S 1$
再標本2 = $X 1, X 1, X 3, X 4, X 4$ $\to$ 統計量 = $S 2$
再標本3 = $X 2, X 2, X 5, X 5, X 5$ $\to$ 統計量 = $S 3$
......
再標本998 = $X 1, X 2, X 3, X 4, X 5$ $\to$ 統計量 = $S 998$
再標本999 = $X 1, X 3, X 4, X 5, X 5$ $\to$ 統計量 = $S 999$
再標本1000 = $X 3, X 3, X 3, X 3, X 4$ $\to$ 統計量 = $S 1000$

　たくさん統計量がとれましたね。そして、こう考えます。

「このたくさんの統計量の分布って、その統計量の母集団分布とほぼ同じと考えていいよね！だってどうせもっと実験してデータをたくさん取ったところで、似たような値のデータが増えるだけだからね！！」

　つまりブートストラップでは、実験的にではなく、統計のトリックを使ってデータ数を増やしているのですね。そしてデータ数がある程度大きければ、その分布は限りなく母集団の分布に近いはずだ、という想定をしています。
　
4: 求めた統計量の分布を使って、バラツキ（標準誤差、信頼区間など）を計算する。
　例えば標準誤差は、母集団分布の標準偏差のことなので、単純に、求められた統計量分布の標準偏差を求めればいいですね。

　こうして、めでたく気になる統計量のバラツキを求めることができます。

ブートストラップの実践

　pythonでやってみましょう (Matlab/Octaveも下にあります）。

　まず、必要なライブラリをゲットします。


# import librarys
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

　Numpyを使って、データをアレイに入れます。


# get data
iq = np.array([61, 88, 89, 89, 90, 92, 93, 
               94, 98, 98, 101, 102, 105, 108,
               109, 113, 114, 115, 120, 138])

　基本統計量として、平均値、平均値の標準誤差、中央値を計算します。


# compute mean, SEM (standard error of the mean) and median
mean_iq = np.average(iq)
sem_iq = np.std(iq)/np.sqrt(len(iq))
median_iq = np.median(iq)

　それぞれ、平均値: 100.85、平均値の標準誤差: 3.45、中央値: 99.5 と計算できました。

　中央値の標準誤差を求めるため、ブートストラップを実行します。ブートストラップが正しくできてるかの sanity check として、一応平均値の標準誤差も一緒に計算しています。

　pythonで、1~nの連続した整数値を重複ありでランダムに並び替えるには、以下のようにします。


np.random.choice(n,n,replace=True)

　これを使って、bootstrapの再標本化を行います。


# bootstrap to compute sem of the median
def bootstrap(x,repeats):
    # placeholder (column1: mean, column2: median)
    vec = np.zeros((2,repeats))
    for i in np.arange(repeats):
        # resample data with replacement
        re = np.random.choice(len(x),len(x),replace=True)
        re_x = x[re]
            
        # compute mean and median of the "new" dataset
        vec[0,i] = np.mean(re_x)
        vec[1,i] = np.median(re_x)
    
    # histogram of median from resampled datasets
    sns.distplot(vec[1,:], kde=False)
    
    # compute bootstrapped standard error of the mean,
    # and standard error of the median
    b_mean_sem = np.std(vec[0,:])
    b_median_sem = np.std(vec[1,:])
    
    return b_mean_sem, b_median_sem

　この関数を x = iq, repeats = 1000 として実行すると、ブートストラップによって求められた1000個の中央値の分布をみることができます。

　……もっと繰り返しを増やした方がよかったかもですね。

　とりあえず、この分布から中央値の標準誤差を求めることができます。単に1000個の中央値の標準偏差を計算するだけです。

ブートストラップで求められた中央値の標準誤差: 4.22

　同じように、平均値の標準誤差も計算します。

ブートストラップで求められた平均値の標準誤差: 3.42

　最初に、元データから求めた平均値の標準誤差は3.45なので、結構近いですね。ブートストラップは正しく統計量のバラツキを推定したようです。

注意点
　ブートストラップではランダムに再標本化するため、実行するたびに結果が少しずつ異なります。結果のバラツキを減らし、より正確な推定をするには、繰り返し数を増やす (~10,000) 必要があります。

まとめのまとめ

　ブートストラップでは再標本化により、興味のある統計量のバラツキを推定することができる。

終わりに

　ブートストラップをはじめ、統計学の中でも再標本化 (resampling) を考えた人は天才だと思います。データを人工的に増やして母集団分布を近似してしまうなんて、とりあえずデータをひたすら集め続ける苦行から解放されるわけですね！もちろん、ある程度データを実際にとって母集団の代表例を集めなければ、いくら再標本化しても永遠に母集団には近づけませんが……。

　ソースコードは以下に載せておきます。

bootstrap_demo.py


# import librarys
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns

# get data
iq = np.array([61, 88, 89, 89, 90, 92, 93, 
               94, 98, 98, 101, 102, 105, 108,
               109, 113, 114, 115, 120, 138])

# compute mean, SEM (standard error of the mean) and median
mean_iq = np.average(iq)
sem_iq = np.std(iq)/np.sqrt(len(iq))
median_iq = np.median(iq)

# bootstrap to compute sem of the median
def bootstrap(x,repeats):
    # placeholder (column1: mean, column2: median)
    vec = np.zeros((2,repeats))
    for i in np.arange(repeats):
        # resample data with replacement
        re = np.random.choice(len(x),len(x),replace=True)
        re_x = x[re]
            
        # compute mean and median of the "new" dataset
        vec[0,i] = np.mean(re_x)
        vec[1,i] = np.median(re_x)
    
    # histogram of median from resampled datasets
    sns.distplot(vec[1,:], kde=False)
    
    # compute bootstrapped standard error of the mean,
    # and standard error of the median
    b_mean_sem = np.std(vec[0,:])
    b_median_sem = np.std(vec[1,:])
    
    return b_mean_sem, b_median_sem   

# execute bootstrap
bootstrapped_sem = bootstrap(iq,1000)

Matlab/Ovctaveは以下です。

bootstrap_demo.m


function bootstrap_demo

% data
iq = [61, 88, 89, 89, 90, 92, 93,94, 98, 98, 101, 102, 105, 108,109, 113, 114, 115, 120, 138];

% compute mean, SEM (standard error of the mean) and median
mean_iq = mean(iq);
sem_iq = std(iq)/sqrt(length(iq));
median_iq = median(iq);
disp(['the mean: ' num2str(mean_iq)])
disp(['the SE of the mean: ' num2str(sem_iq)])
disp(['the median: ' num2str(median_iq)])
disp('---------------------------------')

[b_mean_sem, b_median_sem] = bootstrap(iq, 1000);
disp(['bootstrapped SE of the mean: ' num2str(b_mean_sem)])
disp(['bootstrapped SE of the median: ' num2str(b_median_sem)])

% bootstrap to compute sem of the median
function [b_mean_sem, b_median_sem] = bootstrap(x, repeats)

% placeholder (column1: mean, column2: median)
vec = zeros(2,repeats);
for i = 1:repeats
    % resample data with replacement
    re_x = x(datasample(1:length(x),length(x),'Replace',True));
    
    % compute mean and median of the "new" dataset
    vec(1,i) = mean(re_x);
    vec(2,i) = median(re_x);
    
end

% histogram of median from resampled dataset
histogram(vec(2,:))

% compute bootstrapped standard error of the mean, and standard error of
% the median
b_mean_sem = std(vec(1,:));
b_median_sem = std(vec(2,:));

参考

The bootstrap method for standard errors and confidence intervals

You get articles that match your needs
You can efficiently read back useful information
You can use dark theme

What you can do with signing up

Why not login to Qiita and try out its useful features?

色んな統計量のバラツキを求めたい！再標本化 (Bootstrap)のススメ

ブートストラップとは

ブートストラップのアルゴリズム

ブートストラップの実践

まとめのまとめ

終わりに

参考

@katsu1110

Today's trending articles

メモリの動的確保はなぜ嫌われるのか？

セキュリティエンジニアの知識地図：読んでみた！

プログラムに1文字加えて本番環境のバッチ処理を4倍遅くした

「7日間でハッキングをはじめる本」をやってみての感想

とりあえずGUIが欲しいならStreamlitがおすすめ

Comments

Qiita Advent Calendar is held!

Login to continue?

Login or Sign up with social account

Login or Sign up with your email address

Why not login to Qiita and try out its useful features?

色んな統計量のバラツキを求めたい！再標本化 (Bootstrap)のススメ

ブートストラップとは

ブートストラップのアルゴリズム

ブートストラップの実践

まとめのまとめ

終わりに

参考

@katsu1110

Today's trending articles

メモリの動的確保はなぜ嫌われるのか？

セキュリティエンジニアの知識地図：読んでみた！

プログラムに1文字加えて本番環境のバッチ処理を4倍遅くした

「7日間でハッキングをはじめる本」をやってみての感想

とりあえずGUIが欲しいならStreamlitがおすすめ

Comments

parkQiita Advent Calendar is held!

ブートストラップとは　

Qiita Advent Calendar is held!