ロジスティック回帰

サポートベクターマシンは、線形しきい素子を用いた識別器であるが、 Rosenblatt が提案したパーセプトロンも、同様に、線形しきい素子を用いて、 訓練サンプルから学習する識別機械である。サポートベクターマシンと同様に、 単純パーセプトロンでは、入力 $\mbox{\boldmath$x$}=(x_1,\ldots,x_M)^T$ に対する出力 $y$ を

$$y = f(\eta)$$

$$\eta = \mbox{\boldmath$w$}^T \mbox{\boldmath$x$} - h = \tilde{\mbox{\boldmath$w$}}^T \tilde{\mbox{\boldmath$x$}}$$ (37)

のように計算する。ここで、$w_i$ は、$i$番目の入力から出力への結合荷重で あり、$h$ はバイアスである。以下では、簡単のために、これらをまとめて、 $\tilde{\mbox{\boldmath$w$}}=(h,w_1,\ldots,w_M)^T$ のように表すものとする。また、入力 特徴ベクトルに定数項を加えたベクトルを $\tilde{\mbox{\boldmath$x$}}=(-1,x_1,\ldots,x_M)^T$ と表す。出力ユニットの入出力関数 $f$ は、Rosenblatt のオリジナルなモデルではしきい関数

$$f(\eta) = \mbox{sign}(\eta) = \left\{ \begin{array}{ll} 1 ... ... $\eta \geq 0$} \\ 0 & \mbox{otherwise} \end{array} \right.$$ (38)

が用いられた。この他の入出力関数としてロジスティック関数

$$f(\eta) = \frac{\exp(\eta)}{1+\exp(\eta)}$$ (39)

がよく使われる。多変量データ解析の用語を用いれば、出力ユニットの入出力関 数がロジスティック関数の単純パーセプトロンは、ロジスティック回帰モデルで ある。ロジスティック回帰では、普通、最尤法によりパラメータが推定される。

今、訓練サンプル集合を $\{(\mbox{\boldmath$x$}_i,u_i)\vert i=1,\ldots,N\}$ とする。ここで は、教師信号 $u_i$ は、$0$ か $1$ の2値で与えられるものとする。

入力 $\mbox{\boldmath$x$}$ を与えたときの出力 $y$ を、入力 $\mbox{\boldmath$x$}$ のもとで教師信号 $u$ が $1$ である確率の推定値と考えると、訓練サンプル集合に対するネット ワークの尤度は、

$$L = \prod_{i=1}^N y_i ^{u_i} (1 - y_i)^{(1 - u_i)}$$ (40)

で与えられる。従って、その対数(対数尤度)は、

$$l = \sum_{i=1}^N \{u_i \log y_i + (1 - u_i) \log (1 - y_i) \}$$

$$= \sum_{i=1}^N \{u_i \log \{ \frac{\exp(\eta_i)}{1+\exp(\eta_i)} \}$$

$$+ (1-u_i) \log \{ \frac{1}{1+\exp(\eta_i)} \} \}$$

$$= \sum_{i=1}^N \{u_i \eta_i - \log \{1+\exp(\eta_i) \} \}$$ (41)

となる。これを最大とするパラメータがネットワークの最尤推定値である。

まずは、最急降下法によりパラメータ $\tilde{\mbox{\boldmath$w$}}$を遂次更新することで最 適なパラメータを求める方法について考えてみよう。対数尤度のパラメータ $w_j$に関する偏微分は、

$$\frac{\partial l}{\partial w_j} = \sum_{i=1}^N (u_i - y_i) x_{ij} = \sum_{i=1}^N \delta_i x_{ij}$$ (42)

のようになる。ここで $\delta_i = (u_i- y_i)$ である。一方、対数尤度のパ ラメータ$h$に関する偏微分は、

$$\frac{\partial l}{\partial h} = \sum_{i=1}^N (u_i - y_i) (-1) = \sum_{i=1}^N \delta_i (-1)$$ (43)

となる。したがって、パラメータの更新式は、

$$w_j \Leftarrow w_j + \alpha (\sum_{i=1}^N \delta_i x_{ij})$$ (44)

$$h \Leftarrow h + \alpha (\sum_{i=1}^N \delta_i (-1))$$ (45)

のようになる。

最尤推定においては、Fisher 情報行列が重要な役割を演じる。一般に、 データ $y$ がパラメータ $\theta_1,\ldots,\theta_M$ をもつ 密度関数 $f(y,\theta_1,\ldots,\theta_M)$ をもつ分布に従うとき、

$$F_{ij} = - E \left( \frac{\partial^2}{\partial \theta_i \partial \theta_j} \log f(y,\theta_1,\ldots,\theta_M) \right)$$ (46)

を Fisher 情報量と呼び、行列 $F = [ F_{ij} ]$ を Fisher 情報行列という。 Fisher 情報量は不変推定量の分散と密接に関係している。

ここでは、ロジスティック回帰の Fisher 情報量を具体的に計算する。そのため には、式 (41) の対数尤度の2次微分を計算する必要がある。対数尤 度の2次微分は、

$$\frac{\partial^2 l}{\partial \tilde{w}_k \partial \tilde{w}_j} = - \sum_{i=1}^N \omega_i \tilde{x}_{ik} \tilde{x}_{ij}$$ (47)

となる。ただし、 $\omega_i = y_i (1 - y_i)$ である。1次微分と2次微分をま るとめて行列表現すると

$$\nabla l = \sum_{i=1}^N \delta_i \tilde{\mbox{\boldmath$x$}}_i = X^T \mbox{\boldmath$\delta$},$$ (48)

$$\nabla^2 l = - \sum_{i=1}^N \omega_i \tilde{\mbox{\boldmath$x$}}_i \tilde{\mbox{\boldmath$x$}}_i^T = - X^T W X$$

となる。ただし、 $X^T = [\tilde{\mbox{\boldmath$x$}}_1,\ldots,\tilde{\mbox{\boldmath$x$}}_N]$ , $W = diag(\omega_1,\ldots,\omega_N)$ および $\mbox{\boldmath$\delta$}=(\delta_1,\ldots,\delta_N)^T$ である。

これらを用いて、パラメータ $\tilde{\mbox{\boldmath$w$}}$ に対する Fisher 情報行列、す なわち、Hessian 行列の期待値のマイナスは、

$$F = -E(\nabla^2 l) = X^T W X$$ (49)

となる。これは、入力ベクトル $\{\tilde{\mbox{\boldmath$x$}}_i\}$ の $\omega_i$ で重み 付けた相関行列である。そのときの重み $\omega_i$ は、図6に示 すような２次関数で、ニューロンの出力が確定している（ $0$ あるいは $1$ に近い）場合には小さくなり、出力が不確定な（ $0.5$ に近い）場合には大 きくなる。従って、Fisher 情報行列は、主に、出力が不確定な入力ベクトル の相関行列であると考えることができる。

図 6: 重み $\omega _p$

対数尤度(41)を最大とするようなパラメータを求めるためには、非線形最適 化法を用いる必要がある。ロジスティック回帰では、このために Fisher のスコアリ ングアルゴリズムが使われる [68]。これは、一種のニュートン法で、 Hessian 行列のかわりに Fisher 情報行列を用いる。ニューロン１個のみからなるネッ トワークの場合、Fisher 情報行列とHessian 行列は単に符合が異なるだけなので、 Fisher のスコアリングアルゴリズムはニュートン法そのものとなる。

今、現時点でのパラメータの推定値を $\mbox{\boldmath$w$}$ とし、それを修正ベクトル $\delta\tilde{\mbox{\boldmath$w$}}$ により、

$$\tilde{\mbox{\boldmath$w$}}^* = \tilde{\mbox{\boldmath$w$}} + \delta \tilde{\mbox{\boldmath$w$}}$$ (50)

のように更新する ものとする。修正ベクトル $\delta\tilde{\mbox{\boldmath$w$}}$ は、線形方程式

$$F \delta \tilde{\mbox{\boldmath$w$}} = \nabla l$$ (51)

を解くことにより求められる。パラメータの更新式(50) に左から $F$ を掛けると、

$$F \tilde{\mbox{\boldmath$w$}}^* = F \tilde{\mbox{\boldmath$... ...mbox{\boldmath$w$}} = F \tilde{\mbox{\boldmath$w$}} + \nabla l$$ (52)

となる。今、 $F \tilde{\mbox{\boldmath$w$}}$ は、

$$F \tilde{\mbox{\boldmath$w$}} = X^T W \mbox{\boldmath$\eta$}$$ (53)

となる。ただし、 $\mbox{\boldmath$\eta$} = (\eta_1,\ldots,\eta_N)^T$ である。 従って、新しい推定値 $\tilde{\mbox{\boldmath$w$}}^*$ は、

$$\tilde{\mbox{\boldmath$w$}}^* = F^{-1} (F \tilde{\mbox{\boldmath$w$}} + \nabla l)$$

$$= (X^T W X)^{-1} (X^T W \mbox{\boldmath$\eta$} + X^T \mbox{\boldmath$\delta$})$$

$$= (X^T W X)^{-1} X^T W (\mbox{\boldmath$\eta$} + W^{-1} \mbox{\boldmath$\delta$})$$ (54)

により求めることができる。ただし、 $\mbox{\boldmath$\delta$}=(\delta_1,\ldots,\delta_N)^T$ である。この式は、入力データ $X$ から目的 変数 $\mbox{\boldmath$\eta$} + W^{-1} \mbox{\boldmath$\delta$}$ への重み付き最小２乗法の正規方程式 とみなすことができる。従って、最尤推定値を求めるには、ある初期値からはじ めて、この重み付き最小２乗法を繰り返せばよいことになる。

上記のアルゴリズムは、繰り返しアルゴリズムであるためパラメータの初期値が 必要である。これは、例えば、以下のような簡単な方法で推定することが可能で ある。今、結合重みがすべて $0$ 、つまり、 $\mbox{\boldmath$w$} = \mbox{\boldmath$0$}$ とする。この とき、 $W = \frac{1}{4} I$, $\mbox{\boldmath$\eta$} = \mbox{\boldmath$0$}$ および $\mbox{\boldmath$\delta$} = \mbox{\boldmath$u$} - \frac{1}{2} \mbox{\boldmath$1$}$ である。従って、これらを (54) の計算 式に代入すると、初期パラメータの推定値 $\tilde{\mbox{\boldmath$w$}}_0$ は、

$$\tilde{\mbox{\boldmath$w$}}_0 = 4 (X^TX)^{-1} X^T (\mbox{\boldmath$u$} - \frac{1}{2} \mbox{\boldmath$1$})$$ (55)

となる。これは、初期パラメータを入力から教師信号 $\mbox{\boldmath$t$}-\frac{1}{2}\mbox{\boldmath$1$}$ への線形回帰により求めることに対応している。