Edited at

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

More than 1 year has passed since last update.

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

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

 それがデータの平均値であれば、標準誤差 (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だとしましょう。分布は以下のような感じです。

figure1.png

 このとき、平均値 (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 = $X1,X2,X3,X3,X4$

再標本2 = $X1,X1,X3,X4,X4$

再標本3 = $X2,X2,X5,X5,X5$

......

再標本998 = $X1,X2,X3,X4,X5$

再標本999 = $X1,X3,X4,X5,X5$

再標本1000 = $X3,X3,X3,X3,X4$

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

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

3: 各再標本で、求めたい統計量を計算する。

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

再標本1 = $X1,X2,X3,X3,X4$ $\to$ 統計量 = $S1$

再標本2 = $X1,X1,X3,X4,X4$ $\to$ 統計量 = $S2$

再標本3 = $X2,X2,X5,X5,X5$ $\to$ 統計量 = $S3$

......

再標本998 = $X1,X2,X3,X4,X5$ $\to$ 統計量 = $S998$

再標本999 = $X1,X3,X4,X5,X5$ $\to$ 統計量 = $S999$

再標本1000 = $X3,X3,X3,X3,X4$ $\to$ 統計量 = $S1000$

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

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

 つまりブートストラップでは、実験的にではなく、統計のトリックを使ってデータ数を増やしているのですね。そしてデータ数がある程度大きければ、その分布は限りなく母集団の分布に近いはずだ、という想定をしています。

 

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個の中央値の分布をみることができます。

figure2.png

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

 とりあえず、この分布から中央値の標準誤差を求めることができます。単に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,:));



参考