GLMをもう少し理解したい②
前回の記事からだいぶ間が空いてしまいましたが続きを書いてみます。なおこの記事は主にDobsonの「一般化線形モデル入門」の第3・4章を参考にしていますので、そちらも合わせてご確認ください。良書です。
- 作者: Annette J.Dobson,田中豊,森川敏彦,山中竹春,冨田誠
- 出版社/メーカー: 共立出版
- 発売日: 2008/09/08
- メディア: 単行本
- 購入: 15人 クリック: 152回
- この商品を含むブログ (13件) を見る
さて、前回の記事では、以下のように結びました:
つまり
glm
は最尤法によって解を推定していると思われている(そしてそれは正しい)のだけれど、実際には最小二乗解をHouseholder法によって得ているのだということがわかりました。
これは果たしてどういうことなのでしょうか?これについて答えるために以下の流れで見ていきます。
※1/22 記事を少し修正しました
ニュートン・ラプソン法
初めにニュートン・ラプソン法についてです。やや天下り的で申し訳ないのですが、一般にGLMでは解を推定する際に数値的な解法を必要とします(理由については後述します追記:考えていた説明が誤っていたので削除します)。その際に用いられる手法としてニュートン・ラプソン法というものがあり、これは非線形な方程式の解を反復によって求める手法です。
あまり詳しい説明はできませんが、原理的には以下の式を反復的に更新することで求めたい解を得ることができます:
ここでおよびはそれぞれ更新前後の解となります。 は対象となる関数、はその導関数を表します。以下の例で簡単に見てみましょう:
x <- seq(-3, 3, 0.1) y <- sin(x) plot(x, y, type = "l")
上記のようなsin関数においてとなる点を見つけます。先ほど示したニュートン・ラプソンのアルゴリズムに従い、
- 初期値を入力
- 関数およびその導関数による返り値を計算する
- 解を更新する
- 更新後の関数値が十分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回の反復でほぼに収束している様子がわかります。この例の場合はにおいてつまりが0となることが示されました1。
スコア法
上記では一般的なニュートン・ラプソン法の原理を紹介しましたが、本記事のテーマであるGLMに置き換えると、およびはそれぞれおよびが対象となります。ここではパラメータ、即ち推定したい対象であり、は対数尤度関数の導関数です。
なぜ対数尤度関数そのものではなく対数尤度関数の導関数なのでしょうか?これについては以下のような図を考えるとわかりやすいかもしれません。
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")
ここで黒い線は対象となる関数(二次関数)、赤い線はその導関数、青い点線はy = 0を示しています。
先ほどニュートン・ラプソン法の節において
となる点を見つけます
と書きましたが、このアルゴリズムでは返り値が0となる点を見つける一方、我々が得たいのは対数尤度関数を最大にする点です。このグラフで言えば黒い線を最大にするようなです。
もしこの黒い線に対してニュートン・ラプソン法を当てはめるとどうなるでしょうか?おそらく黒い線と青い線が交わるところに収束するでしょう。先ほどのスクリプトを少し修正してみます:
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
となるため、確かにが0となる点を見つけているようです。しかしこれは最大となる値ではありません。
それでは赤い線で試してみるとどうでしょうか?答えはグラフからも明らかですが、元々の目的である黒い線の最大値で収束しそうです。
上記のような理由により、GLMにニュートン・ラプソンを当てはめる場合には、対数尤度関数そのものではなく対数尤度関数の導関数を用いると言えそうです。
それでは改めて、GLMにおけるニュートン・ラプソン法の推定方程式を見てみましょう:
ここでは対数尤度関数をで微分したもので、スコア関数、スコア統計量、または単にスコアなどと呼ばれます。
はスコアを更に微分したものとなりますが、最尤推定においてはそのものではなくにより近似したものを用いることがあります。このとき、
という関係があるため2、先程の方程式は以下のように書き直されます:
さらに推定対象をパラメータのベクトルに置き換えると:
と一般化され、この式を用いて解を得る方法をスコア法と言います。
最尤推定
それでは上の式がそれぞれどのように得られるのかを確認していきますが、ここで一度切りたいと思います。