1. katsu1110

    No comment

    katsu1110
Changes in body
Source | HTML | Preview
@@ -1,268 +1,267 @@
 データ解析はまだまだブームのようです。機械学習であれ何であれ、まずデータが与えられてすることは、データを視覚化して、解析のとっかかりとなる仮説をつかむことですね。仮説を検証するため、様々な解析をして様々な統計量を計算します。 
 そんなことをしていると、**ある統計量の母集団におけるバラツキ**が知りたい、という事態に出くわすことがあります。
 それがデータの平均値であれば、標準誤差 (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$
 こうして、ステップ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も下にあります)。
 まず、必要なライブラリをゲットします。
```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 allowing overlapping
- re_x = np.zeros(len(x))
- for re in np.arange(len(x)):
- # get one data point after shuffling
- np.random.shuffle(x)
- re_x[re] = x[0]
+ 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 allowing overlapping
- re_x = np.zeros(len(x))
- for re in np.arange(len(x)):
- # get one data point after shuffling
- np.random.shuffle(x)
- re_x[re] = x[0]
+ # 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 allowing overlapping
- re_x = zeros(1,length(x));
- for re = 1:length(x)
- % get one data point after shuffling
- shu = randperm(length(x));
- re_x(re) = x(shu(1));
- end
-
+ % 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/)