最尤法による一般化線形モデルのパラメータ推論

【概要】

一般化線形モデル(ロジスティック回帰、ポアソン回帰)のパラメータ推論を最尤法を使って実装してみます
確率モデルとして考えることで統一した考え方ができます（わざわざ「〇〇回帰」みたいな名称を覚える必要ない）

【目次】

はじめに
統計モデリング
- 一般化線形モデル
- 最尤法
ロジスティック回帰
ポアソン回帰
実装(jupyter notebook)
まとめ
参考文献

はじめに

「一般化線形モデル」の実用例として頻繁に使われる（と思う）ロジスティック回帰とポアソン回帰のパラメータを最尤法によって推論してみます。「ロジスティック回帰」、「ポアソン回帰」のようにそれぞれ名前が付いているように思いますが、これらは確率モデルの一つです。確率モデルとして考えれば統一したアプローチでモデルの推論ができます。

本記事は以下のセミナーで話した内容を整理したものです。

liberal-arts-for-tech.connpass.com

次回はこちらです。

liberal-arts-for-tech.connpass.com

なおここでの話題は、こちらのテキストで取り扱っています（宣伝）。

【トップに戻る】

統計モデリング

「統計モデリング」とは、データの構造を確率分布を組み合わせて表現する活動だと理解しています。そして、確率分布を組み合わせて構築するモデルのことは「確率モデル」と呼ぶことが多いと思います。（解釈が間違っているかも。間違っていたらご指摘お願いします。）

例えば、私が趣味でやっているもう一つのブログ「チョコボール統計」では、エンゼルの出現有無というデータをベルヌーイ分布に従うと仮定し、データからそのパラメータを推論しています。

他にも、線形関数にガウスノイズが乗っていると考える「線形回帰*1」などがよく使われると思います。

一般化線形モデル

あくまで説明上の都合ですが、”一般化”ではない「線形モデル」は次のような確率モデルです。

$\begin{aligned} &y_n \sim p(y_n | x_n, W) = \mathcal{N}(y_n | f(x, w), \sigma^2) \\ &f(x, w) = w_0 + w_1 x_n \end{aligned}$

変数 $y$ がガウス分布に従い、ガウス分布のパラメータが重みベクトル $\mathbf{w}=\{w_0, w_1\}$ と変数 $\mathbf{x}=\{1, x\}$ の線型結合になっているモデルです。

次に、「一般化線形モデル」では、変数 $y$ が従う分布がガウス分布以外にも様々な分布になり得ます。「一般化」とはこれだけです。例えば、あるイベントが発生するかしないかといった2値のデータを扱う場合には、ベルヌーイ分布に従うと仮定します。ベルヌーイ分布のパラメータはイベントの発生確率なので、0~1の実数です。上記の「線形モデル」のように $\mathbf{w}$ と $\mathbf{x}$ の線型結合とすると0~1の範囲を超えてしまう可能性があるので、ロジスティック・シグモイド関数をかけるような調整を行います。これが「ロジスティック回帰」と呼ばれているものです。

ということで、「一般化線形モデル」のように名称がついていますが、あくまでも確率モデルです。確率変数がどのような分布に従うと仮定するか？、その分布のパラメータの定義域は？を考えればわざわざ名前を付けて別に扱う必要ないです*2。

最尤法

確率モデルの未知のパラメータを推論するための考え方として、尤度を基準にしてパラメータを推論する手法がよく使われます。

尤度 $L(\theta)$ というのは、データと確率モデルがどの程度フィットしているかを測るための指標です。データに対しての同時確率で表現されます。

$\begin{aligned} &L(\theta) = \prod^N_{n=1} p(y_n | \theta_n) \\ &\theta_n = f(x_n, W) \end{aligned}$

ここで、 $\theta$ は確率モデルのパラメータであり、変数xの関数です。先にあげたチョコボールの例では、xは無く常に一定のパラメータと考えています。線形回帰モデルの場合は $f(x, w_0, w_1)=w_0 + w_1x$ のような関数が使われます。ここで、 $W=(w_0, w_1)$ がデータから求めたい未知のパラメータです。

なお、上式のように同時確率を総積で表現できるのは、独立同分布(i.i.d.)を仮定しているためです。つまり、データの順番は関係ないということです。時系列データのように順番に依存する場合はこのようには書けません。

尤度は推論対象のパラメータ $\mathbf{W}$ の関数です。尤度が高いほどデータとモデルがフィットしているので、尤度が最も大きくなるようにパラメータ $\mathbf{W}$ を調整してやるのが「最尤推定」という考え方になります。

最大値を求めるということで $\mathbf{W}$ について微分して0となる $\mathbf{W}$ を導出すれば良いと考えると思いますが、後述のように活性化関数として非線形な関数を使っており解析解が導出できないケースがあります。そこで、勾配法などの探索的なアルゴリズムによってパラメータを推論します。

勾配法とは、誤差関数（対数尤度の符号反転）の微分に基づいて、少しずつパラメータを最適値に近づけていこうとする方法です。手順としては次のようになります。

尤度関数を定義する
- $L(\theta) = \prod^N_{n=1} p(y_n | \theta_n)$
誤差関数を定義する
- $E(\theta)$ は負の対数尤度
- $E(\theta) = - \log L(\theta) = - \sum \log(p(y_n | \theta_n))$ を最大化する
$E(\theta)$ の微分 $\nabla E(\theta)$ を算出
に小さい係数（学習係数:）を掛けてパラメータを更新する
- $\mathbf{W}_{t+1} = \mathbf{W}_{t} + \alpha \nabla E(\theta)$

手順4を何度も繰り返し、パラメータを更新します。

【トップに戻る】

ロジスティック回帰

ロジスティック回帰は、二値変数 $y \in \{0, 1\}$ を予測するモデルです。 $y=1$ の出現確率を $\theta$ とし、 $\theta$ は変数 $x$ に依存する構造になります。確率を回帰しているモデルと見做せます。

例えば、二値変数 $y$ としてWebのバナー広告がクリックされるか否かを考えてみます。「クリック」というイベントが発生する確率を $\theta$ としますが、広告の商材を扱うWebサイトを何日前に訪れたかといった変数 $x$ に依存して $\theta$ は変わるのではないか？と考えられます。このような関係をモデリングするのに、ロジスティック回帰はよく使われると思います。実際は $x$ として多変量の値が用いられることが多いですね。

f:id:hippy-hikky:20200529152944p:plain — ロジスティック回帰の例．確率θのベルヌーイ分布に従う変数y(青点)，θは変数xに依存して変化する(赤点)．

モデル

$y \in \{0, 1\}$ なので、 $y$ が従う確率分布としてベルヌーイ分布を考えます。ベルヌーイ分布のパラメータ $\theta$ を変数 $x$ とパラメータ $\mathbf{w}$ の線形結合で予測します。ここで、 $\theta$ は確率なので0~1の範囲に収まっている必要があります。なので、ロジスティック・シグモイド関数で0~1の範囲に収めてしまいます。

f:id:hippy-hikky:20200529153706p:plain:w200

ここで、 $\mathbf{w}$ は重みベクトルで、ここでは $\mathbf{w}=[w_0, w_1]^T$ の二次元ベクトルとします。 $\mathbf{x}$ は $\mathbf{x}=[1, x]^T$ とします。つまり、 $\mathbf{w}\mathbf{x}=w_0+w_1x$ という線形回帰の式になります*3。

ちなみに、一般化線形モデルの説明でよく出てくるリンク関数というのは、上記の $\theta$ の式を変形することで現れます。ロジスティック回帰の場合は下式の通り、ロジット関数になります。

$\begin{aligned} \mathrm{logit}(\theta) = \log\frac{\theta}{1-\theta} = \mathbf{w}\mathbf{x} \end{aligned}$

このようになるのですが、私としては「リンク関数」を考えるよりも、パラメータを確率モデルに合うように押し込めるための変換と考えた方がわかりやすい気がします。このような見方は「活性化関数」という表現をされることが多いです。

パラメータ推論

先に説明した勾配法を使ってパラメータを推論します。

まずは、勾配を求めるために、誤差関数を定義します。誤差関数は、負の対数尤度ですので以下の通りです。

$\begin{aligned} E(\mathbf{w}) &= - \log L(\theta) \\ &= - \sum^{N}_{n=1} \left\{ y_n \log \theta_n + (1-y_n) \log (1-\theta_n) \right\} \end{aligned}$

ということでクロスエントロピーが出てきました。

あとはこれをパラメータ $\mathbf{w}$ について微分すれば良いです。微分の結果は以下のようになります。

f:id:hippy-hikky:20200529171812p:plain:w300

ということで、ここで出てきた勾配に基づいて、パラメータを更新していけば良いわけです。

モデルの定義と勾配の導出までの手書きメモを貼っておきます。汚い字なので読めないかも。

f:id:hippy-hikky:20200530223555p:plain:w200 f:id:hippy-hikky:20200530223611p:plain:w200 f:id:hippy-hikky:20200530223625p:plain:w200

実装

勾配が計算できればあとは単純にパラメータを更新していくだけです。具体的な実装については、下記のnotebookを参照ください。

先に挙げたサンプルデータを推論した結果は以下のようになりました。真の関数をかなり忠実に推論できています。

f:id:hippy-hikky:20200529172323p:plain — ロジスティック回帰の最尤推定結果．青点がデータ(x, yのペア)，赤点が真の関数．緑点が推定結果．

学習曲線はこんな感じです。

f:id:hippy-hikky:20200529172526p:plain — ロジスティック回帰の学習曲線．

【トップに戻る】

ポアソン回帰

ポアソン回帰は、0を含んだ正の整数( $y$ )を予測するモデルです。後述の通りyはポアソン分布に従うとします。ポアソン分布のパラメータはyの平均 $\lambda$ なので、平均 $\lambda$ が $x$ に依存する構造になります。 yはイベントの発生回数などカウント系のデータです。

例えば、かき氷の販売数を $y$ として、気温を $x$ として気温とかき氷の販売数の関係をモデリングするような応用が考えられます*4。

f:id:hippy-hikky:20200529173640p:plain — ポアソン回帰のデータ例．平均回数λのポアソン分布に従う変数y(青点)．λは変数xに依存して変化する(赤点)．

モデル

$y$ は0を含んだ正の整数、つまり、カウント系のデータですので、 $y$ が従う確率分布としてポアソン分布を考えます。ポアソン分布のパラメータは上記の通り平均 $\lambda$ です。 $\lambda$ が変数 $x$ に依存するとして、重み $\mathbf{w}$ との線型結合で予測をします。ここで、 $\lambda$ は平均回数ということなので、0以上の実数であることが必要です。そこで、指数関数に押し込めてしまいましょうというのがポアソン回帰です。

f:id:hippy-hikky:20200529180614p:plain:w200