いい加減時代の潮流に乗ろうということで機械学習を学びたいと思う。またRはともかくとしてPythonは未経験であるため、丁度良い書籍として
を用いることにする。
4. 情報量基準
- モデル選択において適合性と簡潔性のバランスが取れた統計モデルの指標である情報量基準を扱う。
4.1 情報量基準
情報量基準は観測データから統計モデルの妥当性を評価する場合の指標として定義される。情報量基準は統計モデルによるデータの説明力(適合力)と統計モデルの簡潔度合い(簡潔性)*1の両方を評価する。
線形回帰を想定する*2と、組の観測値
を用いて
個の変数候補から説明変数として実際に採用する変数を選ぶ問題として定式化される。すなわち
について
を選ぶ問題として定式化される。この
を用いた場合の残差変動
が適合性を、変数の個数
が簡潔性を表す。このとき、
はそれぞれ
により定義される。
4.2 有効推定量とFisher情報量行列
を導出すべく最小二乗法による線形回帰の推定方法が不偏推定量の中でも分散を最小にする有効推定量であることを導く。
観測データが定数
および確率変数
を用いて
で表されるものとする。すなわち
である。
ここで最尤法、すなわち尤度
を最大にするようなを考える。計算を簡単にすべく、また対数
が
に関して単調増加であることを踏まえ、対数尤度
を最大化することを考えることにする。すると
であり、を固定すれば、対数尤度の最大化と
の最小化は同値となる。また対数尤度
を
で偏微分すれば
であるから、を代入することで
の最尤推定量
は上記の偏微分を
とおいて
を得る。
次にが有効推定量であることを示す。一般に対数尤度
を
の各成分で微分した値から成り立つベクトル
の分散共分散行列を
で割った
を
情報量行列という。
の両辺をで偏微分することで
を得る。
また
および
が成り立つ。これは
を意味する。
(の両辺をで偏微分すると
である。これを分散共分散行列の形式で書けば
となる。
またより、これは
を意味する。そこでと
をつなぎ合わせたサイズ
の確率ベクトルの分散共分散行列は
である。ここではともに分散共分散行列であるから非負である。したがって
である。右辺が非負定値であるから、左辺も非負定値である。任意のについて
ならば任意の
と
について
が成立する。これはが非負定値であることを意味する。実際、
について
ならば、においても不等号が成り立つ。
)
4.3 Kullback-Leibler情報量
上の確率密度関数
に対して
を-
情報量といい、任意の
上の
集合について
すなわち絶対連続ならばは定義できる。
線形回帰の
は任意のについて
が成り立つ。
である。したがって
が成立し、に関して総和を取ることで
が成り立つ。 )
がそれぞれ
から発生したと仮定するとき、
に関する対数尤度の平均は、
の平均が
、
の平均が
であることに注意すれば
を得る。この値は
と書け、これはに依らない定数であるから、
-
情報量の和
を最小するようにを選べばよい。
4.4 AICの導出
ここでは赤池情報量基準を導出する。
真の母数(ベクトル)をとおき、その不偏推定量の中で
-
情報量の和
を最小にするような母数の推定量を選ぶ。
確率変数に対して
の不等式から
が成り立つ。-
の不等式から任意の
の推定量
について対数尤度を
とおけば、実際、
に関する2次方程式
において解は存在しても高々1つであるから、判別式を考えれば、とおけば
が得られる。
において
であり、またこの右辺のトレースを取ればであるから
が成り立つ。さらに
である。以上から
が得られ、等号が成立する。
は、第2項をその最小値で置き換えた
を最小にすることを目的とする。特に変数選択の問題では、個すべての変数を用いるのではなく、
個の変数を選択する。すなわち
を最小にするようなを選択する。
このときは未知である。そこで変数の部分集合
に対して
で代用する方法が考えられる。しかし
は
と比較して平均的に小さくなる。
標本誤差の平均値の評価
が成り立つ。
が成り立つ。またであるから、各
において
である。さらにの定理から
が成立する。ここでより
が成り立つ。
この式の右辺の第項について二項定理から
が得られる。
まずを
に関する多項式と見れば、最高次の項の係数は
で
次の係数は
であるから、
の
次および
次の項の係数は
を得る。したがって
が成り立つ。
またおよび
の定理の適用結果から、
を適用して
であるから、
が得られる。 )
以上から、の項を除いて
が成り立つ。そこでを
に置き換えることで
を得、これを最小にするようなを選ぶのが
である。
4.5 実証
4.5.1 Rスクリプト
さまざまなパッケージ内にを求める関数が存在する。今回は参考書P.81を基にしたものと、
に内蔵された
とで計算してみる*3。
# 残差変動を最小にするような変数の組み合わせを取得する RSS_min <- function(mt_x,vc_y,mt_t){ m <- ncol(mt_t) RSS_min <- Inf # 初期値 for(j in 1:m){ q <- mt_t[,j] # 回帰 ls_reg <- lm(vc_y~mt_x[,q]) RSS <- sum((ls_reg$fitted.values -vc_y)^2)/n # もし残差変動がより小さいならばその値を控えておく if(RSS < RSS_min){ RSS_min <- RSS set_q <- q set_AIC <- AIC(ls_reg) } } return(list(value = RSS_min, set = set_q, AIC = set_AIC)) } ################################ ### Bostonデータを用いた実証 ### ################################ library(MASS) df <- Boston # テストデータ mt_x <- as.matrix(df[,c(1,3,5:8,10:13)]) vc_y <- df[[14]] p <- ncol(mt_x) n <- length(vc_y) AIC_min <- Inf df_AIC_results <- data.frame() # 各結果(AICの値)を保存 ls_variables <- list() # 各結果(変数) for(k in 1:p){ mt_t <- combn(1:p, k) res <- RSS_min(mt_x, vc_y, mt_t) nm_AIC <- n * log(res$value / n) + 2 * k # AIC if(nm_AIC < AIC_min){ AIC_min <- nm_AIC set_min <- res$set } # 各結果を保存しておく df_AIC_results <- rbind(df_AIC_results,c(k,nm_AIC,res$AIC)) ls_variables[[k]] <- res$set } print(AIC_min) # 最小のときのAIC print(set_min) # 最小のときの説明変数
4.5.2 Pythonスクリプト
################################################### ### データセットからAICを用いたモデル選択を行う ### ################################################### import numpy as np from sklearn.linear_model import LinearRegression import itertools res = LinearRegression() ### 残差変動が最も小さくなる説明変数の組み合わせを計算する def RSS_min(X, y, T): S_min = np.inf # 初期値 m = len(T) for j in range(m): q = T[j] res.fit(X[:,q], y) # 回帰 y_hat = res.predict(X[:,q]) # 予測値 S = np.linalg.norm(y_hat - y) ** 2 # 残差変動 if S < S_min: S_min = S set_q = q return(S_min,set_q) ### 実証 # テストデータを取得 from sklearn.datasets import load_boston boston = load_boston() X = boston.data[:,[0,2,4,5,6,7,9,10,11,12]] y = boston.target n, p = X.shape # Xの行数と列数 AIC_min = np.inf # 初期値 for k in range(1,p+1,1): T = list(itertools.combinations(range(p), k)) # p個からk個を選ぶ組み合わせ S_min, set_q = RSS_min(X, y, T) AIC = n * np.log(S_min/n)+2 * k # AIC if AIC < AIC_min: AIC_min = AIC set_min = set_q print(AIC_min, set_min)
*1:具体的には説明変数の数。
*2:これは説明上の都合でそうしたのであり、情報量基準は一般のモデルに対して適用できる。
*3:前述したようには真のモデルからの距離に近いもの(Kullback-Leibler情報量)として測っている。が、真のモデルの位置合いが不明なので相対値でしかない。そのため、計算パッケージによっては算出ロジックが相違する可能性がある。そのため絶対的な値を異なるパッケージの出力を比較しても意味が無い。同じパッケージ(計算ロジック)で計算したものを比較すること。実際、以下のスクリプトでは、自作した
の計算結果とdf_AIC_resultsの3列目に保存した値は値が異なる。しかしこれらで選択したモデルは同一である。