「大人の教養・知識・気付き」を伸ばすブログ

一流の大人(ビジネスマン、政治家、リーダー…)として知っておきたい、教養・社会動向を意外なところから取り上げ学ぶことで“気付く力”を伸ばすブログです。目下、データ分析・語学に力点を置いています。今月(2022年10月)からは多忙につき、日々の投稿数を減らします。

MENU

統計的機械学習の数理100問(12/20)

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

を用いることにする。

3. リサンプリング

  • 真の統計モデルよりも複雑なモデルを選択してしまう過学習による影響を受けずに学習の性能を評価する方法であるクロスバリデーションを学ぶ。
  • 学習結果のバラつきを評価する汎用的な方法であるブートストラップを学ぶ。

3.3 線形回帰におけるCV

 k-\mathrm{fold}クロスバリデーションで特に分割数k=nと標本数と一致させると処理に時間が掛かる。以下、線形回帰の場合に高速なクロスバリデーションの実現方法を議論する。
 相互に重ならない\{1,\cdots,n\}kの部分集合を考える。そのうちの1つをテストデータとし、それ以外を訓練データとするのがクロスバリデーションである。テストデータとなる部分集合を変えてk回行って得られる評価値の総和は以下で計算できることが知られている:


線形回帰のクロスバリデーション評価値
 \{1,\cdots,n\}の分割S_1,\cdots,S_k,S_i\cap S_j=\emptyset,i\neq j,\displaystyle{\bigcup_i S_i}=Sに対して、そのクロスバリデーションの評価値の和は

\begin{aligned}
\displaystyle{\sum_{S_k}\left\|(I-H_{S_k})^{-1}e_{S_k}\right\|^2}
\end{aligned}

で与えられる。ここでH_{S_k}=X_{S_k}({}^{t}XX)^{-1}{}^{t}X_{S_k}H=X({}^{t}XX)^{-1}{}^{t}XS_kに所属する行と列からなる部分行列、e_{S_k}は残差\boldsymbol{e}=Y-\boldsymbol{X}\hat{\boldsymbol{\beta}}S_kに含まれている行からなるベクトルである。また\|\cdot\|^2は各成分の二乗和を表すものとする。

 この命題の成否について検討する。まずX\in\mathbb{R}^{n\times(p+1)}Sに含まれる行から成る行列をX_S\in\mathbb{R}^{r\times p}、それ以外のXの行から成る行列をX_{-S}\in\mathbb{R}^{(n-r)\times p}とする。ここでrは集合Sの要素数とする。同様に\boldsymbol{Y}\in\mathbb{R}^nについても\boldsymbol{Y}_S,\boldsymbol{Y}_{-S}を定義する。このとき

\begin{aligned}
{}^{t}XX&=\displaystyle{\sum_{j=1}^{n}\boldsymbol{X}_j{}^{t}\boldsymbol{X}_j}\\
&=\displaystyle{\sum_{j\in S}\boldsymbol{X}_j{}^{t}\boldsymbol{X}_j}+\displaystyle{\sum_{j\notin S}\boldsymbol{X}_j{}^{t}\boldsymbol{X}_j}\\
&={}^{t}X_SX_S+{}^{t}X_{-S}X_{-S}
\end{aligned}

および


\begin{aligned}
{}^{t}X\boldsymbol{Y}&=\displaystyle{\sum_{j=1}^{n}{}^{t}x_j y_j}=\displaystyle{\sum_{j\in S}{}^{t}x_j y_j}+\displaystyle{\sum_{j\notin S}{}^{t}x_j y_j}\\
&={}^{t}X_S\boldsymbol{Y}_S+{}^{t}X_{-S}\boldsymbol{Y}_{-S}
\end{aligned}

が成立する。
 次に以下が成り立つ:


\mathrm{Sherman}-\mathrm{Morroson}-\mathrm{Woodbury} m,n\geq1、行列A\in\mathbb{R}^{n\times n},U\in\mathbb{R}^{n\times m},C\in\mathbb{R}^{m\times m},V\in\mathbb{R}^{m\times n}について

\begin{aligned}
(A+UCV)^{-1}=A^{-1}-A^{-1}U\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}
\end{aligned}

が成り立つ。

(\because 

\begin{aligned}
&(A+UCV)\left(A^{-1}-A^{-1}U\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}\right)\\
=&I+UCVA^{-1}-U\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}+\\
&UCVA^{-1}U\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}\\
=&I+UCVA^{-1}-UC(C^{-1})\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}+\\
&UC(VA^{-1}U)\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}\\
=&I+UCVA^{-1}-UC\left(C^{-1}-VA^{-1}U\right)\left(C^{-1}+VA^{-1}U\right)^{-1}VA^{-1}\\
=&I+UCVA^{-1}-UCVA^{-1}\\
=&I
\end{aligned}

が成り立つ。 \blacksquare)

 そしてn=p+1,m=r,A={}^{t}XX,C=I,U={}^{t}X_S,V=-X_Sを適用すれば、


\begin{aligned}
\left({}^{t}X_{-S}X_{-S}\right)^{-1}=({}^{t}XX)^{-1}+({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}X_S({}^{t}XX)^{-1}
\end{aligned}

が成立する。ここで以下の定理を用いた。


部分行列の正則性 {}^{t}XXが正則であるとき、各S\subset\{1,\cdots,n\}に対して[tex:{}^{t}X_{-S}X_{-S}が正則であれば、I-H_Sも正則である。
(\because m,n\geq1,U\in\mathbb{R}^{m\times n},V\in\mathbb{R}^{n\times m}について

\begin{aligned}
\begin{bmatrix}
I&O\\
V&I
\end{bmatrix}
\begin{bmatrix}
I+UV&U\\
O&I
\end{bmatrix}
\begin{bmatrix}
I&O\\
\left.-V\right.&I
\end{bmatrix}&=\begin{bmatrix}
I+UV&U\\
V+VUV&VU+I
\end{bmatrix}
\begin{bmatrix}
I&O\\
\left.-V\right.&I
\end{bmatrix}\\&=
\begin{bmatrix}
I&U\\
O&I+VU
\end{bmatrix}
\end{aligned}

である。このとき


\begin{aligned}
\mathrm{det}(I+UV)=\mathrm{det}(I+VU)
\end{aligned}

であることを踏まえれば、


\begin{aligned}
\mathrm{det}\left({}^{t}X_{-S}X_{-S}\right)&=\mathrm{det}\left({}^{t}XX-{}^{t}X_{-S}X_{-S}\right)\\
&=\mathrm{det}({}^{t}XX)\mathrm{det}(I-({}^{t}XX)^{-1}{}^{t}X_S X_S)\\
&=\mathrm{det}({}^{t}XX)\mathrm{det}(I-X_S({}^{t}XX)^{-1}{}^{t}X_S)
\end{aligned}

を得る。なお最後の変形ではU=({}^{t}XX)^{-1}{}^{t}X_S,V=-X_Sを用いた。以上から、{}^{t}X_{-S}X_{-S}および{}^{t}XXが正則ならば、I-H_Sもまた正則である。 \blacksquare)

 したがって以上からSに属するデータを用いずに推定した回帰係数の推定値\hat{\boldsymbol{\beta}}_{-S}=\left({}^{t}X_{-S}X_{-S}\right)^{-1}{}^{t}X_{-S}\boldsymbol{Y}_{-S}


\begin{aligned}
\hat{\boldsymbol{\beta}}_{-S}&=\left({}^{t}X_{-S}X_{-S}\right)^{-1}\left({}^{t}X\boldsymbol{Y}-{}^{t}X_S\boldsymbol{Y}_S\right)\\
&=\left\{({}^{t}XX)^{-1}+({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}X_S({}^{t}XX)^{-1}\right\}
\left({}^{t}X\boldsymbol{Y}-{}^{t}X_S\boldsymbol{Y}_S\right)\\
&=\hat{\boldsymbol{\beta}}-({}^{t}XX)^{-1}{}^{t}X_S\boldsymbol{Y}_S-({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}(X_S\hat{\boldsymbol{\beta}}-H_S\boldsymbol{Y}_S)\\
&=\hat{\boldsymbol{\beta}}-({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}\{(I-H_S)\boldsymbol{Y}_S-X_S\hat{\boldsymbol{\beta}}+H_S\boldsymbol{Y}_S\}\\
&=\hat{\boldsymbol{\beta}}-({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}\boldsymbol{e}_S
\end{aligned}

を満たす。ここで\hat{\boldsymbol{\beta}}=({}^{t}XX)^{-1}{}^{t}X\boldsymbol{Y}は全データを用いた\boldsymbol{\beta}の推定値であり、\boldsymbol{e}_S=\boldsymbol{Y}_S-X_S\hat{\boldsymbol{\beta}}は、その\hat{\boldsymbol{\beta}}を用いてSに属するデータで評価した場合の誤差に相当する。
 Sに属さないデータで推定した回帰係数\hat{\boldsymbol{\beta}}_{-S}Sに属するデータを適用して残差および残差平方和を評価すると


\begin{aligned}
\boldsymbol{Y}_S-X_S\hat{\boldsymbol{\beta}}_{-S}&=\boldsymbol{Y}_S-X_S\{\hat{\boldsymbol{\beta}}-({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}\boldsymbol{e}_S\}\\
&=\boldsymbol{Y}_S-X_S\hat{\boldsymbol{\beta}}+X_S({}^{t}XX)^{-1}{}^{t}X_S(I-H_S)^{-1}\boldsymbol{e}_S\\
&=\boldsymbol{e}_S+H_S(I-H_S)^{-1}\boldsymbol{e}_S\\
&=(I-H_S)^{-1}\boldsymbol{e}_S
\end{aligned}

と計算できる。
 n個の標本すべてを用いた推定値\hat{\boldsymbol{\beta}}の残差がそれぞれ\boldsymbol{e}_Sであるのに対してSに属さない標本を用いて推定された\hat{\boldsymbol{\beta}}_{-S}の残差はそれぞれ(I-H_S)^{-1}\boldsymbol{e}_S(I-H_S)^{-1}倍されることになる。
 予測誤差を求めるには、事前にH=X{}^{t}XX)^{-1}{}^{t}Xおよび\boldsymbol{e}=(I-H)\boldsymbol{Y}を求めておく。そして各Sに対して対応する行および列を選んでH_S,\boldsymbol{e}_Sが計算できる。そして\|(I-H_S)^{-1}\boldsymbol{e}_S\|^2Sについて和を取ればよい。

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 ブートストラッピング

 p変量の標本i=1,2,\cdots,nを格納したデータセット\boldsymbol{df}={}^{t}(df_{1},\cdots,df_{n})\in\mathbb{R}^{n\times p}があるときにその各行df_{1},\cdots,df_{n}から統計量\alphaの推定量fを用いて


\begin{aligned}
\hat{\alpha}=f(df_{1},\cdots,df_{n})
\end{aligned}

が得られるものとする。このとき、標本df_{1},\cdots,df_{n}から重複を許して無作為に選んだn個からなる新しいデータセットdf_{i_1},\cdots,df_{i_n},i_n\in\{1,2,\cdots,n\}から


\begin{aligned}
\hat{\alpha}_1=f(df_{i_1},\cdots,df_{i_n})
\end{aligned}

が得られる。これをr回繰り返すことで\hat{\alpha}_h,h=1,2,\cdots,rの不偏分散として


\begin{aligned}
\displaystyle{\frac{1}{r-1}\sum_{h=1}^{r}(\hat{\alpha}_h-\hat{\alpha})^2}
\end{aligned}

が得られる。

3.3.1 より詳細な理論

 ブートストラップ法の実行プロセスを、推定量のバイアス、分散および推定量の確率分布とパーセント点の推定を通じて述べる。
 未知の確率分布Fにしたがって生成されたn個のデータを\boldsymbol{x}=\{x_1,\cdots,x_n\}とする。これらは互いに独立に分布Fに従う確率変数\boldsymbol{X}=\{X_1,\cdots,X_n\}の実現値とする。いま確率分布Fに従う確率変数Fに関する母数\thetaをとある推定量\hat{\theta}=\hat{\theta}(\boldsymbol{X})を用いて推定する。推定の誤差を捉える評価尺度として以下のバイアスおよび分散


\begin{aligned}
b(F)&=E_{F}\left[\hat{\theta}\right]-\theta,\\
\sigma^2(F)&=E_{F}\left[\left\{\hat{\theta}-E_F\left[\hat{\theta}\right]\right\}^2\right]
\end{aligned}

とが用いられる。
 また\hat{\theta}-\thetaの分布およびパーセント点


\begin{aligned}
G(x)&=P_{F}\{\hat{\theta}-\theta\leq x\},\\
c_{\alpha}&\in\mathbb{R}\ s.t.\ P_{F}\{\hat{\theta}-\theta\leq c_{\alpha}\}=\alpha
\end{aligned}

も評価に用いることができる。

   (1) ブートストラップ標本の反復抽出
      観測データ\boldsymbol{x}=\{x_1,\cdots,x_n\}からn個のデータを復元抽出によって取り出したブートストラップ標本を\boldsymbol{x}^{*}=\{x_1^{*},\cdots,x_n^{*}\}とする。このようなブートストラップ標本をB\{\boldsymbol{x}_i^{*};i=1,2,\cdots,B\}生成する。これらを用いてB個の推定値を計算し、\{\hat{\theta}^{*}(i);i=1,2,\cdots,B\}とする。
   (2) ブートストラップ誤差推定
      観測データに基づく推定値\hat{\theta}B個のブートストラップ推定値\{\hat{\theta}^{*};i=1,2,\cdots,B\}を用いて
\begin{aligned}b(\hat{F})&\approx\displaystyle{\frac{1}{B}\sum_{i=1}^{B}\hat{\theta}^{*}(i)-\hat{\theta}},\\\sigma^2(\hat{F})&\approx \displaystyle{\frac{1}{B-1}\sum_{i=1}^{B}\left\{\hat{\theta}^{*}(i)-\hat{\theta}^{*}(\cdot) \right\}^2}\end{aligned}
ここで\hat{\theta}^{*}(\cdot)=\displaystyle{\frac{1}{B}\sum_{i=1}^{B}\hat{\theta}^{*}(i)}とする。
   (3) ブートストラップ確率分布推定
      推定量の分布と100\alpha%点を以下で数値的に近似計算する:
\begin{aligned}\hat{G}(x)\approx&\displaystyle{\frac{1}{B}\left\{B個の\hat{\theta}^{*}(i)-\hat{\theta}のうちx以下の値の数\right\}},\\\hat{G}^{-1}(\alpha)\approx&\left\{昇順で並べたB個の\hat{\theta}^{*}(i)-\hat{\theta}のうち\right.\\&\left.B\alpha番目の値\right\}\end{aligned}
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)
プライバシーポリシー お問い合わせ