統計コンサルの議事メモ

統計や機械学習の話題を中心に、思うがままに

GLMをもう少し理解したい②

前回の記事からだいぶ間が空いてしまいましたが続きを書いてみます。なおこの記事は主にDobsonの「一般化線形モデル入門」の第3・4章を参考にしていますので、そちらも合わせてご確認ください。良書です。

ushi-goroshi.hatenablog.com

一般化線形モデル入門 原著第2版

一般化線形モデル入門 原著第2版

さて、前回の記事では、以下のように結びました:

つまりglmは最尤法によって解を推定していると思われている(そしてそれは正しい)のだけれど、実際には最小二乗解をHouseholder法によって得ているのだということがわかりました。

これは果たしてどういうことなのでしょうか?これについて答えるために以下の流れで見ていきます。

※1/22 記事を少し修正しました

ニュートン・ラプソン法

初めにニュートン・ラプソン法についてです。やや天下り的で申し訳ないのですが、一般にGLMでは解を推定する際に数値的な解法を必要とします(理由については後述します追記:考えていた説明が誤っていたので削除します)。その際に用いられる手法としてニュートン・ラプソン法というものがあり、これは非線形な方程式の解を反復によって求める手法です。

あまり詳しい説明はできませんが、原理的には以下の式を反復的に更新することで求めたい解を得ることができます:


x^ {(m)} = x^ {(m-1)} - \frac{f(x^ {(m-1)})}{f'(x^ {(m-1)})}

ここで x^ {(m-1)}および x^ {(m)}はそれぞれ更新前後の解となります。  f(x)は対象となる関数、 f'(x)はその導関数を表します。以下の例で簡単に見てみましょう:

x <- seq(-3, 3, 0.1)
y <- sin(x)
plot(x, y, type = "l")

f:id:ushi-goroshi:20190116233606p:plain

上記のようなsin関数において sin(x) = 0となる点を見つけます。先ほど示したニュートン・ラプソンのアルゴリズムに従い、

  1. 初期値を入力
  2. 関数およびその導関数による返り値を計算する
  3. 解を更新する
  4. 更新後の関数値が十分0に近くなったら反復を停止する

ことで最終的な解を得ます。Rで書くと以下のようになるでしょう:

iteration <- 0 # イテレータ
threshold <- 1e-02 # 反復の停止を判断する閾値
x_old <- -1 # 更新前の解
updates <- 1 # 更新後の解による関数値(の絶対値)

while (updates > threshold) {
   x_new <- x_old - sin(x_old)/cos(x_old) # 更新後の解を得る
   updates <- abs(sin(x_new)) # 更新後の解による関数値を得る
   x_old <- x_new # 解を更新する
   cat(sprintf("Iterations: %i, X_New: %f \n", iteration, x_new)) # 表示
   iteration <- iteration + 1 # イテレータを増やす
}

上記のスクリプトを実行すると下のような結果が得られました。

Iterations: 1, X_New: 0.557408 
Iterations: 2, X_New: -0.065936 
Iterations: 3, X_New: 0.000096 

3回の反復でほぼ f(x) = 0に収束している様子がわかります。この例の場合は x = 0において f(x)つまり sin(x)が0となることが示されました1

スコア法

上記では一般的なニュートン・ラプソン法の原理を紹介しましたが、本記事のテーマであるGLMに置き換えると、 xおよび f(x)はそれぞれ \thetaおよび l'(\theta)が対象となります。ここで \thetaはパラメータ、即ち推定したい対象であり、 l'は対数尤度関数導関数です。

なぜ対数尤度関数そのものではなく対数尤度関数の導関数なのでしょうか?これについては以下のような図を考えるとわかりやすいかもしれません。

x <- seq(-2, 3, 0.1)
quad <- function(x) -(x-1)^2 + 3
derive <- function(x) -2 * x + 2
plot(x, quad(x), type = "l", ylim = c(-10, 10), ylab = "")
par(new = T)
plot(x, derive(x), type = "l", col = "red", ylim = c(-10, 10), ylab = "f(x) or f'(x)")
lines(x, rep(0, length(x)), lty = 2, col = "blue")

f:id:ushi-goroshi:20190117003859p:plain

ここで黒い線は対象となる関数(二次関数)、赤い線はその導関数、青い点線はy = 0を示しています。

先ほどニュートン・ラプソン法の節において

 sin(x) = 0となる点を見つけます

と書きましたが、このアルゴリズムでは返り値が0となる点を見つける一方、我々が得たいのは対数尤度関数を最大にする点です。このグラフで言えば黒い線を最大にするような xです。

もしこの黒い線に対してニュートン・ラプソン法を当てはめるとどうなるでしょうか?おそらく黒い線と青い線が交わるところに収束するでしょう。先ほどのスクリプトを少し修正してみます:

iteration <- 0
threshold <- 1e-02
x_old <- 0 # 初期値を変更
updates <- 1

while (updates > threshold) {
   x_new <- x_old - quad(x_old)/derive(x_old) # 関数を変更
   updates <- abs(quad(x_new))
   x_old <- x_new
   iteration <- iteration + 1
   cat(sprintf("Iterations: %i, X_New: %f \n", iteration, x_new))
}
Iterations: 1, X_New: -1.000000 
Iterations: 2, X_New: -0.750000 
Iterations: 3, X_New: -0.732143

-0.732143という点で収束しましたが、quad(-0.732143)-0.0003193724となるため、確かに f(x)が0となる点を見つけているようです。しかしこれは最大となる値ではありません。

それでは赤い線で試してみるとどうでしょうか?答えはグラフからも明らかですが、元々の目的である黒い線の最大値で収束しそうです。

上記のような理由により、GLMにニュートン・ラプソンを当てはめる場合には、対数尤度関数そのものではなく対数尤度関数の導関数を用いると言えそうです。

それでは改めて、GLMにおけるニュートン・ラプソン法の推定方程式を見てみましょう:


\theta^ {(m)} = \theta^ {(m-1)} - \frac{U^ {(m-1)}}{U'^ {(m-1)}}

ここで Uは対数尤度関数を \theta微分したもので、スコア関数スコア統計量、または単にスコアなどと呼ばれます。

 U'はスコアを更に微分したものとなりますが、最尤推定においては U'そのものではなく E(U')により近似したものを用いることがあります。このとき、


E(U') = -Var(U) = -\mathfrak{J}

という関係があるため2、先程の方程式は以下のように書き直されます:


\theta^ {(m)} = \theta^ {(m-1)} + \frac{U^ {(m-1)}}{\mathfrak{J}^ {(m-1)}}

さらに推定対象をパラメータのベクトル \mathbf{\beta}に置き換えると:


\mathbf{b}^ {(m)} = \mathbf{b}^ {(m-1)}  + [\mathfrak{J}^ {(m-1)}]^ {-1} \mathbf{U}^ {(m-1)}

と一般化され、この式を用いて解を得る方法をスコア法と言います。

最尤推定

それでは上の式がそれぞれどのように得られるのかを確認していきますが、ここで一度切りたいと思います。


  1. 仮に初期値を2としていた場合、グラフの右端に向かって値が更新されることになるためx = 3.1415…で収束します。

  2. Dobsonの式(3.16)