統計的機械学習の数理100問（12/20）

　いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として

統計的機械学習の数理100問 with R (with R) (機械学習の数理100問シリーズ)

統計的機械学習の数理100問 with Python (機械学習の数理100問シリーズ)

を用いることにする。

前回

power-of-awareness.com

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)

次回

power-of-awareness.com