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

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

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

　それがデータの平均値であれば、標準誤差 (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](https://qiita-image-store.s3.amazonaws.com/0/122354/2a74fd77-b748-55dc-1d42-3ece76e9e0fd.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$


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

　各再標本にどのデータ標本が入るかはランダムです。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も下にあります）。

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

```python

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

```

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

```python

# 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])

```

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

```python

# 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の連続した整数値を重複ありでランダムに並び替えるには、以下のようにします。

```python

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

```

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

```python

# 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](https://qiita-image-store.s3.amazonaws.com/0/122354/493968ab-5e21-6d6b-94d7-94142e5750d7.png)


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

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

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

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

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

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

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

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

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

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

```python: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は以下です。

```matlab: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](http://www.dummies.com/education/science/biology/the-bootstrap-method-for-standard-errors-and-confidence-intervals/)