反復重み付き最小二乗法で遊ぶ

1．はじめに

積読が増えてきてしまった*1ので、最近はDobson本を読み始めたのですが、反復重み付き最小二乗法（IRLS）なるとても興味深いものを見つけました。面白そうだったからとりあえず実装してみたというのが今回の記事です。本記事は、IRLSを実装するために最低限必要な数式しか載せていないので、文章や式の展開がだいぶ飛躍しています。これに関連して、記事中で出てくる記号について先に述べておくと、 $\displaystyle {\boldsymbol{\beta}}$ は推定したいパラメータであり、 $\displaystyle {\bf{b}}$ は $\displaystyle {\boldsymbol{\beta}}$ の推定値となります。

2．指数型分布族と一般化線形モデル

指数型分布族から述べます。指数型分布族は次の式で表されます。

$\displaystyle f(y;\theta) = exp(a(y)b(\theta)+c(\theta)+d(y))$
独立な確率変数 $\displaystyle Y_1, ..., Y_N$ が従う分布を、指数型分布族に属する分布であると仮定します。さらに、以下に示す関係を仮定します。

$\displaystyle E[Y_i] = \mu_i$
$\displaystyle g(\mu_i) = {\bf{x}}_i^T {\boldsymbol{\beta}} = \eta_i$

ここで、 $\displaystyle {\bf{x}}_i^T$ は $\displaystyle i$ 番目のデータ $\displaystyle Y_i$ に対する説明変数（p×1ベクトル）であり、 $\displaystyle {\boldsymbol{\beta}}$ は推定したいパラメータです。一般化線形モデルでは、連結関数 $\displaystyle g(\mu_i)$ が非常に重要なパーツとなります。

3．パラメータの推定（IRLS）

3.1　とりあえず結論

IRLSの本題に移りますが、冒頭にも述べた通り実装に不要な数式は一気にすっ飛ばしていくので、とりあえずIRLSの反復式を示します。反復式は次のようになります。

$\displaystyle {\bf{X}}^T {\bf{WXb}}^{(m)} = {\bf{X}}^T {\bf{Wz}}$

ここで、 $\displaystyle {\bf{X}}$ は説明変数の行列でありN×p行列です。 $\displaystyle {\bf{b}}$ はパラメータの推定値ベクトルでp×1ベクトルです。 $\displaystyle {\bf{W}}$ はN×Nの対角行列であり、 $\displaystyle (i, i)$ 成分は次の式で表されます（式中の偏微分の箇所は後で説明します）。

$\displaystyle \omega_{ii} = \frac{1}{Var(Y_i)} (\frac{\partial \mu_i}{ \partial \eta_i})^2$

最後に、 $\displaystyle {\bf{z}}$ はN×1ベクトルであり、 $\displaystyle i$ 番目の要素は次の式で表されます。

$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + (y_i - \mu_i)(\frac{\partial \eta_i}{\partial \mu_i})$

IRLSを実装する時は、 $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ を逐一計算していきます。次に、この反復式を導出するに当たって必要な流れを述べていきます。

3.2　導出

最尤推定の流れから導出していきます。指数型分布族と一般化線形モデルの部分で述べた仮定に基づいて、まずは対数尤度関数を求めていきます。 $\displaystyle i$ 番目のデータに対する対数尤度関数 $\displaystyle l_i$ は次の式で表されます。

$\displaystyle l_i = y_ib(\theta_i) + c(\theta_i) + d(y_i)$

$\displaystyle a(y_i)$ が $\displaystyle y_i$ となっているのは、一般化線形モデルの場合は正準形を持つ分布が前提となっているからです。ここからパラメータの最尤推定値を求めていくわけですが、さっくり流れだけ追っていきます。パラメータの最尤推定値を得るために、対数尤度関数を微分します。とりあえず簡単のために、 $\displaystyle i$ 番目のデータに対する対数尤度関数を微分します。

$\displaystyle \frac{\partial l_i}{\partial \beta_j} = \frac{\partial l_i}{\partial \theta_i} \frac{\partial \theta_i}{\partial \mu_i} \frac{\partial \mu_i}{\partial \beta_j}$

対数尤度関数は合成関数になっている*2ため、微分の連鎖法則を用いて微分していきます。さらに、右辺の3つ目の偏微分は次のように微分の連鎖法則で微分できます。

$\displaystyle \frac{\partial \mu_i}{\partial \beta_j} = \frac{\partial \mu_i}{\partial \eta_i} \frac{\partial \eta_i}{\partial \beta_j} = \frac{\partial \mu_i}{\partial \eta_i}x_{ij}$

一番最後の式の偏微分はとりあえず結論で述べた所にも現れています。この偏微分は（記事には載せませんが）スコアや情報行列にも現れてきており、解析的に反復の式を導出する際はこれを計算する必要があります。

あとは偏微分から良い感じになんやかんやしてなんやかんやする*3と、IRLSの反復計算式を求める事が出来ます。途中経過は省略です。

3.3　例．ロジスティック回帰

実験ではロジスティック回帰で遊んでいきます。Dobson本では第7章で触れられています*4が、スコアと情報行列を出して終わっているので、IRLSで計算すべき $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ までは述べてくれていません。従って、自力で導出していきます。計算が間違えている可能性もあるので、あらかじめごめんなさいしておきます。

ひとまず重要なのは、次の2つです。

$\displaystyle \omega_{ii} = \frac{1}{Var(Y_i)} (\frac{\partial \mu_i}{ \partial \eta_i})^2$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + (y_i - \mu_i)(\frac{\partial \eta_i}{\partial \mu_i})$

ここで、 $\displaystyle Y_i～binomial(n_i, \pi_i)$ であり、ロジスティック回帰では期待値と連結関数は次のようになります。

$\displaystyle E[Y_i] = \mu_i = n_i \pi_i$
$\displaystyle g(\pi_i) = \eta_i = ln \frac{\pi_i}{1-\pi_i} = {\bf{x}}_i^T {\boldsymbol{\beta}}$

$\displaystyle w_{ii}$ 中の $\displaystyle Var(Y_i)$ は $\displaystyle n_i \pi_i (1-\pi_i)$ なのでとりあえずOKです。次に、どちらにも共通して現れる偏微分を求めていきます。偏微分の結果は次のようになります。

$\displaystyle \frac{\partial \mu_i}{\partial \eta_i} = \frac{\partial (n_i \pi_i) }{\partial \eta_i} = n_i \frac{\partial \pi_i}{\partial \eta_i} = n_i (\frac{\partial \eta_i}{\partial \pi_i})^{-1} = n_i \pi_i (1-\pi_i)$

第3式から第4式は逆関数の微分を利用しました。これで、 $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ は次のようにすっきり求める事が出来ました。

$\displaystyle \omega_{ii} = n_i \pi_i (1-\pi_i)$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + \frac{y_i - n_i \pi_i }{ n_i \pi_i (1-\pi_i) }$

これで完成で、あとは反復式に入れるだけです。

4．実験

4.1　IRLSの実装

さて、長々と述べてきたわけですが、IRLSの実装に重要なポイントをまとめていきます。

反復式

$\displaystyle {\bf{X}}^T {\bf{WXb}}^{(m)} = {\bf{X}}^T {\bf{Wz}}$
反復式では、 $\displaystyle {\bf{W}}$ と $\displaystyle {\bf{z}}$ を計算する必要があります。

ロジスティック回帰のIRLS

ロジスティック回帰では、 $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ は次の式で表す事が出来ました。

$\displaystyle \omega_{ii} = n_i \pi_i (1-\pi_i)$
$\displaystyle z_i = \sum_{k=1}^{p} x_{ik}b_k^{(m-1)} + \frac{y_i - n_i \pi_i }{ n_i \pi_i (1-\pi_i) }$
あとは上手く反復の手順を考えるだけです。

反復の手順

パラメータが多いと可視化に困るので、 $\displaystyle {\bf{x}}_i$ と $\displaystyle {\bf{b}}$ を次のように決めておきます。

$\displaystyle {\bf{x}}_i^T = [ 1　x_i], {\bf{b}}^T = [ b_1　b_2]$

$\displaystyle m$ ステップ目の推定したパラメータを $\displaystyle \hat{\bf{b}}^{(m)}$ とした時、反復の手順を次に示します。

手順1：初期値を決める→パラメータの初期値を $\displaystyle \hat{ \bf{b} }^{(0)}$ 、確率の最尤推定値の初期値を $\displaystyle {\hat{\pi_i}}^{(0)}$ として与えます。
手順2： $\displaystyle \omega_{ii}$ と $\displaystyle z_i$ を計算
手順3：反復式を次のように変形して、パラメータを推定

$\displaystyle \hat{ \bf{b} }^{(1)} = ({\bf{X}}^T {\bf{WX}} )^{-1} {\bf{X}}^T {\bf{Wz}}$
手順4：パラメータの推定後、 $\displaystyle {\hat{\pi_i}}^{(1)}$ を次の式で計算
$\displaystyle {\hat{\pi_i}}^{(1)} = \frac{exp({\hat{b_1}}^{(1)} + \hat{b_2}^{(1)} x_i )}{1 + exp({\hat{b_1}}^{(1)} + \hat{b_2}^{(1)} x_i ) }$
手順5：手順1の $\displaystyle \hat{ \bf{b} }^{(0)}, {\hat{\pi_i}}^{(0)}$ を $\displaystyle \hat{ \bf{b} }^{(1)}, {\hat{\pi_i}}^{(1)}$ に変更して、手順2から繰り返す

これをパラメータが収束するまで繰り返していきます。

4.2　実験

まず初めにデータセットを作成します。 $\displaystyle \beta_1 = 2.0, \beta_2 = 3.0$ とし、 $\displaystyle (x_i, \pi_i)$ を20個作成しました。この時、生成したデータセットに対してロジスティック回帰モデルを当てはめ、IRLSでパラメータを求める様子をアニメーション化してみました。

左はフィットの様子で、紫のデータは作成したデータセットです。右はパラメータの推移の様子で、青があらかじめ設定したパラメータです。左右共通して赤が推定した結果です。IRLSはあっという間に収束してしまうようで、じわじわパラメータが推定される様子が見えたら良かったのですが、なんだかつまらない感じになってしまいました*5。最終的なパラメータの推定値は、 $\displaystyle b_1 = 1.99, b_2 = 3.01$ でした。

5．まとめ

ちゃんとうごかすことができてとてもたのしかったです

プログラム

こっそりあとがきですが、実はこの記事を書いた段階ではDobson本は第4章までしか読めておらず、収束の判定まではプログラムに実装できませんでした*6。安直にfor文で回しただけのプログラムになってしまいましたが、IRLSのコアとなる部分はfor文の中の4行です。もう少し上手い事プログラムを組めばより短くなるかもしれません。

# seibisi
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from numpy.random import *
from matplotlib import animation

def make_p(beta1, beta2, x):
    return np.exp(beta1 + beta2 * x) / (1 + np.exp(beta1 + beta2 * x))
    
def main():
    # パラメータ
    beta = [2.0, 3.0]
    x = np.linspace(-5, 5, 20)
    pi = make_p(beta[0], beta[1], x)
    # データセットの作成
    n = [round(normal(50, 5)) for i in range(len(x))]
    y = [round(nn * pp) for nn, pp in zip(n, pi)]
    # 初期値
    beta_hat = [0, 0]
    pi_hat = [0.5 for i in range(len(y))]
    X = np.array([np.ones(len(y)), x]).T
    # IRLS
    for i in range(10):
        # 1. zとWの計算
        z = [beta_hat[0] + beta_hat[1] * xx + (yy - nn * pp)/(nn * pp * (1-pp)) for xx, yy, nn, pp in zip(x, y, n, pi_hat)]
        W = np.diag([nn * pp * (1-pp) for nn, pp in zip(n, pi_hat)])
        # 2. パラメータの計算
        beta_hat = np.dot(np.linalg.inv(np.dot(np.dot(X.T, W), X)), np.dot(np.dot(X.T, W), z))
        # 3. 予測
        pi_hat = make_p(beta_hat[0], beta_hat[1], x)

*1:博論を書いているうちに積んでしまっていた

*2:えぐい

*3:ここからスコアと情報行列を導き、情報行列を行列の積で表すと出てきます

*4:実は第4章で読むのが止まっている

*5:7, 8ステップ目で完全に収束しました。pythonかRかなんかを使って収束するまでのステップ数を比較すればもっと面白かったかもしれない

*6:最後まで読んでから実装しよう

遊具

たくさんのおもちゃ

反復重み付き最小二乗法で遊ぶ

1．はじめに

2．指数型分布族と一般化線形モデル

3．パラメータの推定（IRLS）

3.1　とりあえず結論

3.2　導出

3.3　例．ロジスティック回帰

4．実験

4.1　IRLSの実装

反復式

ロジスティック回帰のIRLS

反復の手順

4.2　実験

5．まとめ

プログラム

1．はじめに

2．指数型分布族と一般化線形モデル

3．パラメータの推定（IRLS）

3.1 とりあえず結論

3.2 導出

3.3 例．ロジスティック回帰

4．実験

4.1 IRLSの実装

反復式

ロジスティック回帰のIRLS

反復の手順

4.2 実験

5．まとめ

プログラム

3.1　とりあえず結論

3.2　導出

3.3　例．ロジスティック回帰

4.1　IRLSの実装

4.2　実験