ロジスティック回帰分析!
Updated:
説明変数K個・目的変数1個のロジスティック回帰分析のアルゴリズムについて、プログラムとして実装できるようにするために自分なりに理解してまとめたものです。
ロジスティック回帰分析の実装について調べると、ほとんどが基本的な内容とライブラリを使った実装方法の説明であり、ライブラリに頼らずに自分で実装する方法について知りたいのに、それを詳しく説明しているサイトがなかったので。
0. 前提条件Permalink
- 回帰分析(線型回帰、曲線回帰、重回帰等)についての基本的な知識があること。
1. アルゴリズムPermalink
ロジスティック回帰分析では、求める回帰式として次の式を想定する。
(目的変数1個、説明変数K個、データ数N個を想定)
(ここの p は確率(目的変数; y でもよいが確率なので p とした)、 e はネイピア数)
また、「シグモイド関数」を次の式で定義すると、
σ(x)=11+e(−x)式(1)は次のように表される。
p=σ(β0+β1x1+β2x2+⋯+βKxK)さて、式(1)を変形すると、
p1−p=eβ0+β1x1+β2x2+⋯+βKxK(ここの p1−p は「オッズ(odds)」と呼ばれるもの)
対数をとると、
logp1−p=β0+β1x1+β2x2+⋯+βKxK(ここの logp1−p は「ロジット(logit)」と呼ばれるもの)
また、式(5)(もしくは元の式(1))を変形し、 i 番目のデータについて表現すると、
(また、 e で表現すると指数部分が小さくなるので、以降 exp で表現する)
以下、「最尤法」を用いて係数を推定する。
まず、i番目のデータについて、予測値との当てはまりを考えてみる。
Pi=ytii(1−yi)1−ti(yi:i番目のデータにおける、ロジスティック回帰式から算出される確率(式(1)のp)、
ti:i番目のデータにおける実際の値(0 or 1))
Piが大きい方が予測式の当てはまりが良いとすると、
・ti=1の時:Pi=yiなので、yiは大きい方がいい
・ti=0の時:Pi=1−yiなので、yiは小さい方がいい
ということになる。
この式(7)を全テータに対して適用すると、
L(β)=N∏i=1ytii(1−yi)1−tiこのL(β)を「尤度関数」と呼び、この式が最大になるようにβの値を定めれば、当てはまりが最も良くなると言える。
式(8)の最大化が目的だが、乗算は計算が複雑な上に、確率の総乗を計算しようとすると、値が限りなく0に近付いてしまう。(アンダーフロー) よって、この式(8)の対数をとることによって対処する。
−logL(β)=−N∑i=1{tilogyi+(1−ti)log(1−yi)}対数をとったこの負の式(9)の最小化で元の式(8)の最大化となる。
式(9)の最小値は「勾配降下法」により求める。
まず、次の式を「交差エントロピー誤差関数(⊂損失関数)」とよぶ。(式(9)をE(β)と置いただけ)
E(β)=−logL(β)=−N∑i=1{tilogyi+(1−ti)log(1−yi)}この関数をパラメータβiについて偏微分する。(計算過程省略)
(β0はバイアス(定数)なので、k=0とk∈{1,2,⋯,K}(今回説明変数はK個を想定)で分ける)
∂E(β)∂β0=N∑i=1(yi−ti) (k=0の場合)この式(11),(12)が0になるようにパラメータβを算出すればよいが、解析的には算出できないので、勾配降下法を用いる。
以下のようにして、パラメータβkを更新していく。αは(機会学習で言うところの)学習率で、0.01等の数値を使用する。要は、βkで偏微分した値に学習率を乗じて元のβkから減算していき、値がほとんど変化しなくなったら、値の更新を終了させる、ということ。
βk←βk−α∂E(β)∂βkちなみに、データ数が多いと式(11),(12)の値が大きくなってしまうので、式(11),(12)を次のようにすることもあるようだ。
∂E(β)∂β0=1NN∑i=1(yi−ti) (k=0の場合)近いうちに、プログラムに実装してみたいです。
以上。
Comments