いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として
を用いることにする。
3. リサンプリング
- 真の統計モデルよりも複雑なモデルを選択してしまう過学習による影響を受けずに学習の性能を評価する方法であるクロスバリデーションを学ぶ。
- 学習結果のバラつきを評価する汎用的な方法であるブートストラップを学ぶ。
3.3 線形回帰におけるCV
-
クロスバリデーションで特に分割数
と標本数と一致させると処理に時間が掛かる。以下、線形回帰の場合に高速なクロスバリデーションの実現方法を議論する。
相互に重ならないの
の部分集合を考える。そのうちの1つをテストデータとし、それ以外を訓練データとするのがクロスバリデーションである。テストデータとなる部分集合を変えて
回行って得られる評価値の総和は以下で計算できることが知られている:
線形回帰のクロスバリデーション評価値
で与えられる。ここでは
の
に所属する行と列からなる部分行列、
は残差
の
に含まれている行からなるベクトルである。また
は各成分の二乗和を表すものとする。
および
が成立する。
次に以下が成り立つ:
が成り立つ。
が成り立つ。 )
そしてを適用すれば、
が成立する。ここで以下の定理を用いた。
部分行列の正則性
である。このとき
であることを踏まえれば、
を得る。なお最後の変形ではを用いた。以上から、
および
が正則ならば、
もまた正則である。
)
したがって以上からに属するデータを用いずに推定した回帰係数の推定値
は
を満たす。ここでは全データを用いた
の推定値であり、
は、その
を用いて
に属するデータで評価した場合の誤差に相当する。
に属さないデータで推定した回帰係数
に
に属するデータを適用して残差および残差平方和を評価すると
と計算できる。
個の標本すべてを用いた推定値
の残差がそれぞれ
であるのに対して
に属さない標本を用いて推定された
の残差はそれぞれ
と
倍されることになる。
予測誤差を求めるには、事前におよび
を求めておく。そして各
に対して対応する行および列を選んで
が計算できる。そして
の
について和を取ればよい。
3.3.1 Pythonでの実行例
import numpy as np import time import matplotlib.pyplot as plt %matplotlib inline import japanize_matplotlib def cv_fast(X,y,k): n = len(y) m = int(n/k) H = X@np.linalg.inv(X.T@X)@X.T I = np.diag(np.repeat(1,n)) e = (I - H)@y I = np.diag(np.repeat(1,m)) S = 0 for j in range(k): test = np.arange(j * m, (j + 1) * m, 1, dtype = int) S = S + (np.linalg.inv(I-H[test,test])@e[test]).T@np.linalg.inv(I-H[test,test])@e[test] return S/n cv_fast(X,y,10) # テストデータの生成 n = 1000 p = 5 beta = randn(p + 1) x= randn(n,p) X = np.insert(x,0,1,axis = 1) y = X@beta + randn(n) # U_1 = [] V_1 = [] U_f = [] V_f = [] for k in range(2, n+1, 1): if n%k ==0: # 普通のもの t1 = time.time() # 処理前の時点 cv_linear(X,y,k) t2 = time.time() # 処理後の時点 U_1.append(k); V_1.append(t2-t1) # 高速のもの t1 = time.time() # 処理前の時点 cv_fast(X,y,k) t2 = time.time() # 処理後の時点 U_f.append(k); V_f.append(t2-t1) plt.plot(U_1, V_1, c = "red", label = "cv_linear") plt.plot(U_f, V_f, c = "blue", label = "cv_fast") plt.legend() plt.xlabel("k") plt.ylabel("実行時間(秒)") plt.title("cv_linearとcv_fastの比較")
cv_linearとcv_fastの比較

3.3 ブートストラッピング
変量の標本
を格納したデータセット
があるときにその各行
から統計量
の推定量
を用いて
が得られるものとする。このとき、標本から重複を許して無作為に選んだ
個からなる新しいデータセット
から
が得られる。これを回繰り返すことで
の不偏分散として
が得られる。
3.3.1 より詳細な理論
ブートストラップ法の実行プロセスを、推定量のバイアス、分散および推定量の確率分布とパーセント点の推定を通じて述べる。
未知の確率分布にしたがって生成された
個のデータを
とする。これらは互いに独立に分布
に従う確率変数
の実現値とする。いま確率分布
に従う確率変数
に関する母数
をとある推定量
を用いて推定する。推定の誤差を捉える評価尺度として以下のバイアスおよび分散
とが用いられる。
またの分布およびパーセント点
も評価に用いることができる。
(1) | ブートストラップ標本の反復抽出 | |
観測データ |
||
(2) | ブートストラップ誤差推定 | |
観測データに基づく推定値 |
||
(3) | ブートストラップ確率分布推定 | |
推定量の分布と |
3.3.2 Pythonによる実装例
import numpy as np def bootstrap(df, f, r): m = df.shape[0] org = f(df,np.arange(0,m,1)) u = [] for j in range(r): idx = np.random.choice(m, m, replace = True) # 復元抽出 u.append(f(df,idx)) return {'original':org,'bias':np.mean(u) - org, 'stderr':np.std(u)} # alpha =(V[Y]-V[X])/(V[X]+V[Y]-2Cov[X,Y]) をbootstrapを用いて推定 portfolio = np.loadtxt("C:/Users/Julie/Downloads/TOPIXdiv_SP500.csv", delimiter = ",") def calc_alpha(data, idx): x = data[idx,0] y = data[idx,1] return (np.var(y) - np.var(x))/(np.var(x) + np.var(y)-2 *np.cov(x,y)[0,1]) bootstrap(portfolio, calc_alpha, 10000)