一般化線形モデルのリンク関数になぜ正準パラメータ関数を用いるのか

　一般化線形モデルの正準リンク関数がよくわからなかったので調べてみました。

# 一般化線形モデル
　**一般化線形モデル（Generalized Linear Model; GLM）** とは、線形予測子 $Z = \sum_i \beta_i X_i $ と全単射な単調関数 $g$ を用いて

$$g\left(\mathbb{E}\left[Y\right]\right) = Z$$

と表されるようなモデルのことです。$g$ を**リンク関数**と呼びます。
　確率モデル $f\left(x\mid\theta\right)$ を用いるときは、基本的にリンク関数を $f\left(x\mid\theta\right)$ の正準パラメータ関数 $\eta$ としてモデルを作ります。**正準パラメータ**とは、確率モデルが
$$f\left(x\mid\theta\right) = h\left(x\right) \exp{\left( \eta\left( 
\theta \right) x - A\left(\theta\right) \right)}$$
と書けるときの $\eta\left(\theta\right)$ のことです。例えば、ベルヌーイ分布は

$$
\begin{align}
f\left(x\mid p\right) & = p^x{\left(1-p\right)}^{1-x} \\\\
& = \exp{\left(\log{\left( p^x{\left(1-p\right)}^{1-x} \right)}\right)} \\\\
& = \exp{\left(x\log{\left(p\right)}+\left({1-x}\right)\log{\left(1-p\right)}\right)} \\\\
& = \exp{\left(x\log{\left(p\right)}+\log{\left(1-p\right)} -x\log{\left(1-p\right)}\right)}\\\\
& = \exp{\left(x\log{\left(\frac{p}{1-p}\right)}+\log{\left(1-p\right)}\right)}\\\\
& = \exp{\left(\log{\left(\frac{p}{1-p}\right)} x - \left(-\log{\left(1-p\right)}\right)\right)}
\end{align}
$$

と変形でき、正準パラメータは $p$ をロジット関数に通した
$$\eta\left(p\right) = \log{\left(\frac{p}{1-p}\right)}$$
となります。よって、$0$ から $1$ の間の確率値を求めるような GLM（ロジスティック回帰モデル）はロジット関数を用いて
$$\log{\left( \frac{\mathbb{E}\left[Y\right]}{1-\mathbb{E}\left[Y\right]} \right)} = Z$$
と書けます。ロジット関数の逆関数はロジスティック関数（シグモイド関数）なので、
$$ \mathbb{E}\left[Y\right] = \frac{1}{1 - \exp{\left(-Z\right)}} $$
とも書けます。

# ブレグマン・ダイバージェンス
　しかし、単に確率値を求めるだけならこのようなモデルの作り方をする必然性はなさそうです。リンク関数を正準パラメータから定めたモデルを推定すると、どのようなメリットがあるのでしょうか。確率モデル

$$f\left(x\mid\theta\right) = h\left(x\right) \exp{\left( \eta\left( 
\theta \right) x - A\left(\theta\right) \right)}$$

の対数尤度関数は

$$
\begin{align}
\log{f\left(x\mid \theta\right)} & = \log{\left(h\left(x\right) \exp{\left( \eta\left( 
\theta \right) x - A\left(\theta\right) \right)}\right)} \\\\
& = \log{\left(h\left(x\right)\right)} + \eta\left( 
\theta \right) x - A\left(\theta\right) \\\\
& = \log{\left(h\left(x\right)\right)} + \eta\left( 
\theta \right) x - A'\left(\eta\left(\theta\right)\right)
\end{align}
$$

となります。ただし $A' = A\circ\eta^{-1}$ です。ここで、強凸関数 $\varphi$ についての**ブレグマン・ダイバージェンス（Bregman Divergence）** を以下で定めます。
$$BD_\varphi\left(a||b\right) := \varphi\left(a\right) + \sup_{a}{\left(ab - \varphi\left(b\right) \right)} -ab$$

この式を変形すると、
$$ab - \varphi \left(a \right) = \sup_{a}{\left(ab - \varphi\left(b\right) \right)} - BD_\varphi\left(a||b\right)$$
が得られます。$\varphi = A', a = \eta\left(\theta\right), b = x $ として上式を用いると、

$$
\begin{align}
\log{f\left(x\mid \theta\right)} 
& = \log{\left(h\left(x\right)\right)} + \eta\left( 
\theta \right) x - A'\left(\eta\left(\theta\right)\right) \\\\
& = \log{\left(h\left(x\right)\right)} + \sup_{\eta\left(\theta\right)}{\left(\eta\left(\theta\right) x - A'\left(x\right) \right)} - BD_{A'}\left(\eta\left(\theta\right)|| x\right)
\end{align}
$$

と書き換えることができます。第1項はパラメータ $\theta$ と関係なく、第2項も $\eta\left(\theta\right)$ の定義域さえ定まれば $x$ にしか依存しないため、対数尤度関数を最大化するには $\eta\left(\theta\right)$ と $x$ のブレグマン・ダイバージェンス $BD_{A'}\left(\eta\left(\theta\right)|| x\right)$ を最小化すればよいということがわかります。

　ブレグマン・ダイバージェンスは、我々がよく知っているような損失関数を導きます。例えば正規分布の $x$ と $\eta\left(\mu\right) = \sigma^{-2} \mu$ のブレグマン・ダイバージェンスは二乗誤差
$$\frac{1}{2\sigma^2}{\left(x - \mu\right)}^2$$
となり、ベルヌーイ分布の $x$ と $\eta\left(p\right) = \log{\left({p} \ /\left({1-p}\right)\right)}$ のブレグマン・ダイバージェンスはロジスティック損失
$$ x \log{\left(\frac{x}{p}\right)} + \left(1-x\right) \log{\left(\frac{1-x}{1-p}\right)} = \sum_{i=1}^{2}x_i \log{\left(\frac{x_i}{p_i} \right)}$$
（ただし $x_1 = x, x_2 = 1 - x, p_1 = p, p_2 = 1-p$）となります。また、多項分布からは KL ダイバージェンス
$$ \sum_{i=1}^{n} x_i \log{\left(\frac{x_i}{p_i} \right)}$$
が得られます。

　以上より、「正準パラメータ関数 $\eta$ から作った一般化線形モデルを適切な損失関数（ブレグマン・ダイバージェンス）の最小化により推定する」イコール「最尤推定」ということがわかります。


# 参考文献
- I. Rish, G. Grabarnik. Sparse Modeling: Theory, Algorithms, and Applications. CRC Press, Inc., 2014.