最尤法とGLM(Generalized Linear Model)｜問題演習で理解する統計学【7】

当記事では最尤法とGLM(Generalized Linear Model)の導出に関する演習について取り扱いました。簡単な回帰からロジスティック回帰やポアソン回帰までの計算の流れを詳しく確認できるように作成を行いました。

・標準演習$100$選
https://www.hello-statisticians.com/practice_100

基本問題

基本的な最尤法の流れ

・問題
正規分布$N(\mu, \sigma^2)$に従ってそれぞれ独立にサンプルが得られたと考える。このとき、正規分布$N(\mu, \sigma^2)$の確率密度関数$P(x|\mu,\sigma^2)$は下記のように表すことができる。
$$
\begin{align}
P(x|\mu,\sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2 \sigma^2} \right)
\end{align}
$$
以下の問題に答えよ。
i) $x=x_1, x=x_2$に対応する確率密度関数の値をそれぞれ求めよ。
ⅱ) $x=x_1$と$x=x_2$がそれぞれ独立に観測されることに着目し、同時確率の$P(x=x_1, x=x_2|\mu,\sigma^2)$を表せ。
ⅲ) $P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$と$\log P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$が同じであることを示せ。
iv) $P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$は$\displaystyle \frac{x_1+x_2}{2}$となることを示せ。
v) 同時確率を最大にするパラメータ$\mu$をどのように解釈すると良いか。最尤法の考え方に基づいて述べよ。

・解答
i)
$x=x_1, x=x_2$に対応する確率密度関数の値は下記のように表すことができる。
$$
\large
\begin{align}
P(x=x_1|\mu,\sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_1-\mu)^2}{2 \sigma^2} \right) \\
P(x=x_2|\mu,\sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_2-\mu)^2}{2 \sigma^2} \right)
\end{align}
$$

ⅱ)
AとBが独立した事象の場合$P(A,B) = P(A \cap B) = P(A)P(B)$のように表すことができる。このことを利用することで、$P(x=x_1, x=x_2|\mu,\sigma^2)$は下記のように求めることができる。
$$
\begin{align}
P(x=x_1, x=x_2|\mu,\sigma^2) &= P(x=x_1|\mu,\sigma^2)P(x=x_2|\mu,\sigma^2) \\
&= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_1-\mu)^2}{2 \sigma^2} \right) \times \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_2-\mu)^2}{2 \sigma^2} \right) \\
&= \prod_{i=1}^{2} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_i-\mu)^2}{2 \sigma^2} \right)
\end{align}
$$

ⅲ)
$\displaystyle (\log x)’ = \frac{1}{x}$により、関数の定義域の$x>0$において、$(\log x)’>0$である。これは$\log x$が単調増加関数であることを表す。単調増加関数では、$\log x$が最大になるときと$x$が最大になるときが一致するため、これより「$P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$と$\log P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$が同じであること」を示すことができる。

iv)
$\log P(x=x_1, x=x_2|\mu,\sigma^2)$は下記のように表すことができる。
$$
\begin{align}
\log P(x=x_1, & x=x_2|\mu,\sigma^2) = \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_1-\mu)^2}{2 \sigma^2} \right) \times \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_2-\mu)^2}{2 \sigma^2} \right) \right) \\
&= \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_1-\mu)^2}{2 \sigma^2} \right) \right) + \log \left( \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x_2-\mu)^2}{2 \sigma^2} \right) \right) \\
&= – \frac{(x_1-\mu)^2}{2 \sigma^2} – \frac{(x_2-\mu)^2}{2 \sigma^2} + \rm{Const.}
\end{align}
$$
上記に対し、$\mu$に関して微分を行う。
$$
\begin{align}
\frac{\partial \log P(x=x_1, x=x_2|\mu,\sigma^2)}{\partial \mu} = \frac{x_1-\mu}{\sigma^2} + \frac{x_2-\mu}{\sigma^2}
\end{align}
$$
$\log P(x=x_1, x=x_2|\mu,\sigma^2)$が最大になる場合は、$\displaystyle \frac{\partial \log P(x=x_1, x=x_2|\mu,\sigma^2)}{\partial \mu} = 0$となる。
$$
\begin{align}
\frac{\partial \log P(x=x_1, x=x_2|\mu,\sigma^2)}{\partial \mu} &= 0 \\
\frac{x_1-\mu}{\sigma^2} + \frac{x_2-\mu}{\sigma^2} &= 0 \\
x_1 – \mu + x_2 – \mu &= 0 \\
2 \mu &= x_1 + x_2 \\
\mu &= \frac{x_1+x_2}{2}
\end{align}
$$
よって、「$P(x=x_1, x=x_2|\mu,\sigma^2)$を最大にする$\mu$は$\displaystyle \frac{x_1+x_2}{2}$であること」を示すことができた。

v)
最尤法は「現実の標本は確率最大のものが実現した」と仮定する「最尤原理(principle of maximum likelihood)」という考え方に基づいている。この仮定を用いて、実際に得られた標本の同時確率を尤度とみなし、尤度の最大値問題の形式でパラメータを求める手法を最尤法という。
よって、「$x_1, x_2$は$\displaystyle \mu=\frac{x_1+x_2}{2}$の確率最大のものが実現した」と演繹的に考えることができるし、「$x_1, x_2$が得られる確率を最大にするように$\mu$を求めた」と帰納的に考えることもできる。
最尤法は帰納的にサンプルからパラメータを求める手法であると考えるとわかりやすい。

・解説
この問題では最尤法の基本的な流れについて取り扱いました。尤度や最尤法などは考え方が非常にややこしいので、流れを抑えるのを最優先すると良いと思います。
最尤法の意味合いに関しては「基礎統計学Ⅰ 統計学入門(東京大学出版会)」の$11.2.2$節の記載がわかりやすいので、こちらも合わせて参照すると良いと思います。
また、独立して生成されるサンプルの同時確率は確率密度関数の積の形で表されることが多いため、対数を取り和の形に変形し、微分しやすくするというテクニックを用います。式が複雑になるため、iv)の途中式のように、$\mu$に関係ない項は定数を表す$\rm{Const.}$などで省略すると見やすくなります。

指数型分布族の定義式と具体的な確率分布との対応

・問題
GLM(Generalized Linear Model)を考えるにあたっては、指数型分布族(exponential family)という確率分布を用いる。指数型分布族は正規分布やポアソン分布、ベルヌーイ分布などを総括するような確率分布の集合で、下記のように定義される確率分布が指数型分布族と定義される。
$$
\large
\begin{align}
P(x|\theta) = \exp(a(x)b(\theta)+c(\theta)+d(x))
\end{align}
$$
上記の$\theta$は正規分布の$\mu, \sigma^2$のような確率分布のパラメータを総称したものであると考える。ここまでで確認した指数型分布族の定義に基づいて、以下の問題に答えよ。
i) 正規分布の確率密度関数の式$P(x|\mu, \sigma^2)$を述べよ。
ⅱ) i)で確認した正規分布の式を指数型分布族の定義式の形式に直し、$a(x), b(\mu), c(\mu), d(x)$がそれぞれどのような式となるか答えよ。
ⅲ) ポアソン分布の確率分布の式$P(x|\lambda)$を述べよ。
iv) ⅲ)で確認したポアソン分布の式を指数型分布族の定義式の形式に直し、$a(x), b(\lambda), c(\lambda), d(x)$がそれぞれどのような式となるか答えよ。
v) 正規分布やポアソン分布の最尤法を考えるとき、通常の数式を用いるよりも指数型分布族の式の$P(x|\theta) = \exp(a(x)b(\theta)+c(\theta)+d(x))$に直したものを用いた方が計算がシンプルになるがなぜか。対数と指数が同時に現れる場合の計算に着目しつつ答えよ。

・解答
i)
正規分布の確率密度関数$P(x|\mu, \sigma^2)$は下記のように表される。
$$
\large
\begin{align}
P(x|\mu, \sigma^2) = \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2 \sigma^2} \right)
\end{align}
$$

ⅱ)
$x=\exp(\log(x))$のように変形できることに基づいて、正規分布の確率密度関数$P(x|\mu, \sigma^2)$は下記のように書き直すことができる。
$$
\large
\begin{align}
P(x|\mu, \sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2 \sigma^2} \right) \\
&= \exp \left( \log \left[ \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(x-\mu)^2}{2 \sigma^2} \right) \right]\right) \\
&= \exp \left( -\frac{(x-\mu)^2}{2 \sigma^2} – \frac{1}{2} \log(2 \pi \sigma^2) \right) \\
&= \exp \left( -\frac{x^2 – 2 \mu x + \mu^2}{2 \sigma^2} – \frac{1}{2} \log(2 \pi \sigma^2) \right) \\
&= \exp \left( \frac{\mu x}{\sigma^2} – \frac{x^2}{2 \sigma^2} – \frac{\mu}{2 \sigma^2} – \frac{1}{2} \log(2 \pi \sigma^2) \right)
\end{align}
$$

上記より、$a(x), b(\mu), c(\mu), d(x)$はそれぞれ下記のように表せる。
$$
\large
\begin{align}
a(x) &= x \\
b(\mu) &= \frac{\mu}{\sigma^2} \\
c(\mu) &= – \frac{\mu}{2 \sigma^2} – \frac{1}{2} \log(2 \pi \sigma^2) \\
d(x) &= – \frac{x^2}{2 \sigma^2}
\end{align}
$$

ⅲ)
ポアソン分布の確率関数$P(x|\lambda)$は下記のように表される。
$$
\large
\begin{align}
P(x|\lambda) = \frac{\lambda^{x} e^{-\lambda}}{x!}
\end{align}
$$

iv)
$x=\exp(\log(x))$のように変形できることに基づいて、ポアソン分布の確率関数$P(x|\lambda)$は下記のように書き直すことができる。
$$
\large
\begin{align}
P(x|\lambda) &= \frac{\lambda^{x} e^{-\lambda}}{x!} \\
&= \exp \left( \log \left[ \frac{\lambda^{x} e^{-\lambda}}{x!} \right] \right) \\
&= \exp \left( x \log{\lambda} – \lambda – \log{x!} \right)
\end{align}
$$

v)
下記のように対数尤度の計算の際に生じる式をシンプルに変形を行うにあたって、指数型分布族の表記方法は有用である。
$$
\large
\begin{align}
P(x_1,…,x_n|\theta) &= \prod_{i=1}^{n} P(x_i|\theta) \\
&= \prod_{i=1}^{n} \exp(a(x_i)b(\theta)+c(\theta)+d(x_i)) \\
&= \exp \left( \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i)) \right) \\
\log L(\theta) &= \log{P(x_1,…,x_n|\theta)} \\
&= \log \left[ \exp \left( \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i)) \right) \right] \\
&= \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i))
\end{align}
$$

・解説
指数型分布族を$P(x|\theta) = \exp(a(x)b(\theta)+c(\theta)+d(x))$のように表記を行うと、v)で表したように式をシンプルにまとめられるので、抑えておくと良いと思います。

発展問題

最小二乗法と最尤法

・問題
最尤法では下記のように観測されたサンプルの同時確率$P(x=x_1, …, x=x_n|\theta)$を尤度$L(\theta)$とみなし、尤度の最大化によって$\theta$の値を求める。
$$
\begin{align}
L(\theta) &= P(x=x_1, …, x=x_n|\theta) \\
&= P(x=x_1|\theta)P(x=x_2|\theta)…P(x=x_n|\theta)
\end{align}
$$
このとき、単に標本$x_i$についてのみ考える場合は共通の$\theta$を用いるが、最尤法を回帰に適用する場合は少々論理展開が異なる。パラメータ$\theta$のように抽象的に考えると難しいため、以下では正規分布$N(\mu, \sigma^2)$の$\mu$に主に着目して議論を行う。
$$
\begin{align}
\mu_i = \hat{y}_i = a x_i + b
\end{align}
$$
ここで上記のように定義した$\mu_i$に基づいてサンプル$y_i$が観測されると考えるのであれば、確率密度関数$P(y_i|\mu_i, \sigma^2)$は下記のように表すことができる。
$$
\begin{align}
P(y_i|\mu_i, \sigma^2) &= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(y_i-\mu_i)^2}{2 \sigma^2} \right) \\
&= \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} \right) \\
&= P(y_i|a x_i + b, \sigma^2)
\end{align}
$$
ここまでの内容に基づいて以下の問題に答えよ。
i) 確率密度関数$P(y_i|\mu_i, \sigma^2) = P(y_i|ax_i + b, \sigma^2)$はサンプル$(x_i, y_i)$にどのような前提を付与しているかについて答えよ。
ⅱ) サンプル$(x_1,y_1), …, (x_n,y_n)$が観測されたとき、$a,b$を変数とする尤度関数$L(a,b)$を求めよ。
ⅲ) $a,b$を変数とする対数尤度関数$\log L(a,b)$を求めよ。ただし、偏微分を行うにあたって$a,b$のみを変数と考えるので、$a,b$に関係しない項はConstでまとめて良いものとする。
iv) ⅲ)で求めた$\log L(a,b)$の$a,b$に関する最大化を考える際に、$a,b$のみに着目することでより式をシンプルに表すことができる。このとき、$a,b$に関してどの関数の最大化を行うのが良いか。
v) iv)の結果は最小二乗法の際の誤差関数に一致することに基づき、最小二乗法を行う際に設定する前提について説明せよ。

・解答
i)
$\mu_i=\hat{y}_i=ax_i+b$のように求めた$\mu_i$の周辺に$y_i$が分散$\sigma^2$で正規分布に従って分布すると考えることができる。正規分布の中心の$\mu_i$はサンプルごとに異なる一方で、中心からの分布については各サンプル全てが同様であることに注意が必要である。

ⅱ)
$a,b$を変数とする尤度関数$L(a,b)$はサンプルが観測される同時確率に等しいため、下記のようになる。
$$
\large
\begin{align}
L(a,b) &= \prod_{i=1}^{n} P(y_i|\mu_i, \sigma^2) \\
&= \prod_{i=1}^{n} P(y_i|ax_i + b, \sigma^2) \\
&= \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} \right)
\end{align}
$$

ⅲ)
$L(a,b)$は下記のように整理することができる。
$$
\large
\begin{align}
L(a,b) &= \prod_{i=1}^{n} \frac{1}{\sqrt{2 \pi \sigma^2}} \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} \right) \\
&= \prod_{i=1}^{n} \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} + \rm{Const.} \right)
\end{align}
$$
上記に対して対数尤度関数$\log L(a,b)$を計算する。
$$
\large
\begin{align}
\log L(a,b) &= \log \prod_{i=1}^{n} \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} + \rm{Const.} \right) \\
&= \sum_{i=1}^{n} \log \left( \exp \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} + \rm{Const.} \right) \right) \\
&= \sum_{i=1}^{n} \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} + \rm{Const.} \right)
\end{align}
$$

iv)
$$
\large
\begin{align}
\log L(a,b) &= \sum_{i=1}^{n} \left( -\frac{(y_i-(ax_i+b))^2}{2 \sigma^2} + \rm{Const.} \right)
\end{align}
$$
上記において、$a, b$が関連するのは$(y_i-(ax_i+b))^2$のみである。よって、下記の最大化を考えれば良いことがわかる。
$$
\large
\begin{align}
– \sum_{i=1}^{n} (y_i-(ax_i+b))^2
\end{align}
$$
この結果は最小二乗法で定義する誤差関数に一致する。

v)
最小二乗法は「観測値$y_i$が$ax_i+b$の周りに正規分布で分布する」という前提に基づいていると考える必要がある。よって、予測誤差$\epsilon_i=y_i-(ax_i+b)$が正規分布に従うことが立式の前提にあることを理解しておくことが望ましい。

・解説
ⅲ)で正規分布の式を指数型分布族の式を考慮しながら全てをexpでくくるのは計算ミスを減らすにあたっては有用です。また、iv)、v)のように、最尤法から最小二乗法を導出できるということは抑えておくと理論の理解にあたっては良いので、なるべく抑えておくと良いと思います。

最尤法とGLM

・問題
GLMでは観測値の$(x_i,y_i)$に対して、$x_i$の線形予測子$b_0+b_1x_i$などを考え、リンク関数$g$に基づいて、確率分布の代表値を表すパラメータ$\theta_i$を$\theta_i=g^{-1}(b_0+b_1x_i)$のように求め、$\theta_i$の周りに$y_i$が観測されると考える。

たとえばポアソン回帰に関しては平均を表す$\lambda$に関して$\lambda_i=\exp(b_0+b_1x_i)$を考え、その周囲に下記の式に従って$y_i$が観測されると考える。
$$
\begin{align}
P(y_i|\lambda_i=\exp(b_0+b_1x_i)) = \frac{\lambda_i^{y_i} \exp{(-\lambda_i)}}{y_i!}
\end{align}
$$

ここで上記における推定対象のパラメータは$\lambda_i$ではなく、$b_0, b_1$であることに注意が必要である。$\lambda_i$は各サンプルの$x_i$に対応したポアソン分布の代表値の予測値である。

GLMでは$y_i$の分布に$P(x|\theta) = \exp(a(x)b(\theta)+c(\theta)+d(x))$で表される指数型分布族を用いることで、同時確率に基づく対数尤度を用いた計算をシンプルに記載することが可能となる。

ここまでの内容を元に下記の問いに答えよ。
i) $P(x_i|\theta) = \exp(a(x_i)b(\theta)+c(\theta)+d(x_i))$であることを用いて、同時確率$P(x_1,x_2,…,x_n|\theta)$を計算せよ。
ⅱ) $L(\theta) = P(x_1,x_2,…,x_n|\theta)$のように尤度$L(\theta)$を考えるとき、対数尤度$\log{L(\theta)}$を計算せよ。
ⅲ) ⅱ)の$x_i$を$y_i$、$\theta$を$\theta_i=g^{-1}(b_0+b_1x_i)$で置き換えることで、$\log{L(b_0,b_1)}$を求めよ。
iv) $\displaystyle \frac{\partial \log{L(b_0,b_1)}}{\partial b_0}, \frac{\log{L(b_0,b_1)}}{\partial b_1}$が得られるとき、勾配法を用いてパラメータ推定を行う方法を示せ。
v) 下記を参考に勾配法とニュートン法を比較した際のニュートン法の利点をまとめよ。
https://www.hello-statisticians.com/optimization/newton_grad-desc.html

・解答
i)
同時確率$P(x_1,x_2,…,x_n|\theta)$は下記のように計算できる。
$$
\large
\begin{align}
P(x_1,x_2,…,x_n|\theta) &= \prod_{i=1}^{n} P(x_i|\theta) \\
&= \prod_{i=1}^{n} \exp(a(x_i)b(\theta)+c(\theta)+d(x_i)) \\
&= \exp \left( \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i)) \right)
\end{align}
$$

ⅱ)
対数尤度$\log{L(\theta)}$は下記のように計算できる。
$$
\large
\begin{align}
\log{L(\theta)} &= \log \left[ \exp \left( \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i)) \right) \right] \\
&= \sum_{i=1}^{n} (a(x_i)b(\theta)+c(\theta)+d(x_i))
\end{align}
$$

ⅲ)
$\log{L(b_0,b_1)}$は下記のように計算できる。
$$
\large
\begin{align}
\log{L(b_0,b_1)} &= \sum_{i=1}^{n} (a(y_i)b(\theta_i)+c(\theta_i)+d(y_i)) \\
&= \sum_{i=1}^{n} (a(y_i)b(g^{-1}(b_0+b_1x_i))+c(g^{-1}(b_0+b_1x_i))+d(y_i))
\end{align}
$$

iv)
下記のような漸化式に基づいて繰り返し演算によってパラメータの推定を行うことができる。
$$
\large
\begin{align}
b_0^{(n+1)} &= b_0^{(n)} + \alpha \frac{\partial \log{L(b_0,b_1)}}{\partial b_0} \Bigr|_{b_0=b_0^{(n)},b_1=b_1^{(n)}} \\
b_1^{(n+1)} &= b_0^{(n)} + \alpha \frac{\partial \log{L(b_0,b_1)}}{\partial b_1} \Bigr|_{b_0=b_0^{(n)},b_1=b_1^{(n)}}
\end{align}
$$

v)
学習率を設定する必要がある勾配法に対して、ニュートン法では$2$階微分の値を用いることで学習率の設定が不要であり、このことより収束させやすいと思われることが読み取れる。

・解説
指数型分布族に関する演習の記載と同様に、ⅱ)のような式変形を行うことで計算間違いを減らすことができることは抑えておくと良いです。

ロジスティック回帰の対数尤度の導出

・問題
GLMの$1$例であるロジスティック回帰は$b_0 + b_1 x_i$などで表す回帰の結果に対して、$0$〜$1$の確率が対応するように下記のような変換を行なって予測値の$p_i$の計算を行う。
$$
\begin{align}
p_i = \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \quad (1)
\end{align}
$$
$(1)$式をベルヌーイ分布を元に計算した対数尤度の計算に代入することで、パラメータ$b_0,b_1$に関する対数尤度の導出を行う。

ここまでの内容に基づいて以下の問題に答えよ。
i) $(1)$式において、$b_0=1, b_1=1$のとき、$x_1=-1, x_2=-\infty, x_3=\infty$に対応する$p_1, p_2, p_3$を求めよ。ただし、$x_2, x_3$は極限を求めるだけで良い。
ⅱ) ロジスティック回帰においてリンク関数は$g(p_i)=b_0+b_1x_i$のように表記するが、$(1)$式を元にリンク関数$g(p_i)$を求めよ。
ⅲ) 下記のようにベルヌーイ分布$Bin(1,p)$の確率関数$P(x|p)$を考える際に、$P(x=0|p), P(x=1|p)$をそれぞれ求めよ。
$$
\begin{align}
P(x|p) = p^{x}(1-p)^{1-x} \quad (2)
\end{align}
$$
iv) $(2)$式において$x=\exp(\log{x})$が成立することを利用して、下記の式が成立することを確認せよ。
$$
\begin{align}
P(x|p) = \exp(x \log{p} + (1-x) \log{(1-p)}) \quad (3)
\end{align}
$$
v) 観測値$x_1,x_2,…,x_n$がパラメータ$p$のベルヌーイ分布$Bin(1,p)$に従ってそれぞれ独立に得られる際に、$(3)$式を参考に同時確率$P(x_1,x_2,…,x_n|p)$を計算せよ。
vi) v)で導出した同時確率$P(x_1,x_2,…,x_n|p)$において$p$に着目すると、$p$の尤度と考えることができる。この尤度$L(p)=P(x_1,x_2,…,x_n|p)$の対数の対数尤度$\log{L(p)}$を求めよ。
vⅱ) ロジスティック回帰ではvi)で導出した対数尤度の式の$p$を$p_i$で、$x$を$y_i$で置き換える。このときの尤度を$L(b_0,b_1)$とするとき、対数尤度$\log{L(b_0,b_1)}$を求めよ。

・解答

i)
$(1)$式を用いることで$p_1, p_2, p_3$は下記のように求めることができる。
$$
\large
\begin{align}
p_1 &= \frac{\exp(1 + 1 \cdot (-1))}{1+\exp(1 + 1 \cdot (-1))} \\
&= \frac{e^0}{1+e^0} \\
&= \frac{1}{2} \\
p_2 &= \lim_{x_2 \to -\infty} \frac{\exp(1 + x_2)}{1+\exp(1 + x_2)} \\
&= \frac{0}{1+0} \\
&= 0 \\
p_3 &= \lim_{x_3 \to \infty} \frac{\exp(1 + x_3)}{1+\exp(1 + x_3)} \\
&= \lim_{x_3 \to \infty} \frac{1}{1/\exp(1 + x_3)+1} \\
&= \frac{1}{0+1} = 1
\end{align}
$$

ⅱ)
下記のように$(1)$式を$b_0+b_1x_i$について解くことでリンク関数を導出することができる。
$$
\large
\begin{align}
p_i &= \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \\
p_i(1+\exp(b_0+b_1x_i)) &= \exp(b_0+b_1x_i) \\
p_i &= (1-p_i)\exp(b_0+b_1x_i) \\
\exp(b_0+b_1x_i) &= \frac{p_i}{1-p_i} \\
b_0+b_1x_i &= \log{ \left( \frac{p_i}{1-p_i} \right) }
\end{align}
$$
上記よりリンク関数$g(p_i)$は$\displaystyle g(p_i) = \log{ \left( \frac{p_i}{1-p_i} \right) }$のように表せることがわかる。

ⅲ)
(2)式に$x=0,1$を代入することで下記のように$P(x=0|p), P(x=1|p)$は求めることができる。
$$
\large
\begin{align}
P(x=0|p) &= p^{0}(1-p)^{1-0} \\
&= (1-p)^1 \\
&= 1-p \\
P(x=1|p) &= p^{1}(1-p)^{1-1} \\
&= p^1 \\
&= p
\end{align}
$$

iv)
$x=\exp(\log{x})$が成立することに基づいて、$(2)$式は下記のように変形することができる。
$$
\large
\begin{align}
P(x|p) &= p^{x}(1-p)^{1-x} \\
&= \exp(\log{(p^{x}(1-p)^{1-x})}) \\
&= \exp(x \log{p} + (1-x) \log{(1-p)})
\end{align}
$$

v)
同時確率$P(x_1,x_2,…,x_n|p)$は下記のように求めることができる。
$$
\large
\begin{align}
P(x_1,x_2,…,x_n|p) &= \prod_{i=1}^{n} P(x_i|p) \\
&= \prod_{i=1}^{n} \exp(x_i \log{p} + (1-x_i) \log{(1-p)})
\end{align}
$$

vi)
対数尤度$\log{L(p)}=\log{P(x_1,x_2,…,x_n|p)}$は下記のように求めることができる。
$$
\large
\begin{align}
\log{L(p)} &= \log{P(x_1,x_2,…,x_n|p)} \\
&= \log{ \prod_{i=1}^{n} P(x_i|p) } \\
&= \log{(P(x_1|p) \times P(x_2|p) \times … \times P(x_n|p))} \\
&= \log{P(x_1|p)} + \log{P(x_2|p)} + … + \log{P(x_n|p)} \\&= \sum_{i=1}^{n} \log{P(x_i|p)} \\
&= \sum_{i=1}^{n} \log{(\exp(x_i \log{p} + (1-x_i) \log{(1-p)}))} \\
&= \sum_{i=1}^{n} (x_i \log{p} + (1-x_i) \log{(1-p)})
\end{align}
$$

vⅱ)
vi)の結果における$p$を$p_i$で、$x$を$y_i$で置き換えると下記のようになる。
$$
\large
\begin{align}
\log{L(p)} = \sum_{i=1}^{n} (y_i \log{p_i} + (1-y_i) \log{(1-p_i)})
\end{align}
$$

上記に$(1)$式を代入すると下記のようになる。
$$
\large
\begin{align}
\log{L(p)} &= \sum_{i=1}^{n} (y_i \log{p_i} + (1-y_i) \log{(1-p_i)}) \\
&= \sum_{i=1}^{n} \left( y_i \log{ \left( \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) } + (1-y_i) \log{ \left( 1 – \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) } \right) \\
&= \sum_{i=1}^{n} \left( y_i \log{ \left( \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) } + (1-y_i) \log{ \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) } \right)
\end{align}
$$

・解説
vⅱ)で導出した対数尤度$\log{L(b_0,b_1)}$をパラメータ$b_0, b_1$に関して偏微分を行うことで、ロジスティック回帰vⅱ)で導出した対数尤度$\log{L(b_0,b_1)}$をパラメータ$b_0, b_1$に関して偏微分を行うことで、ロジスティック回帰のパラメータ推定を行うことができます。尤度や対数尤度に用いる大元の式はベルヌーイ分布から導出されますが、ベルヌーイ分布の式はⅲ)のように$x=0, 1$を代入することで理解しやすいことも知っておくと良いです。$p^{x}(1-p)^{1-x}$のように書かれると難しく見えますが、$x=0,1$しか代入されないことを把握すればそれほど難しくないことがわかると思います。

合成関数の微分とロジスティック回帰のパラメータ推定

・問題
$$
\begin{align}
\log{L(b_0,b_1)} &= \sum_{i=1}^{n} \left( y_i \log{p_i} + (1-y_i) \log{(1-p_i)} \right) \\
&= \sum_{i=1}^{n} \left( y_i \log \left( \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) + (1-y_i) \log{ \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) } \right) \quad (1)
\end{align}
$$
$(1)$式で表した前問vⅱ)で導出したパラメータに関する対数尤度$\log{L(b_0,b_1)}$の勾配の計算はなかなか複雑になる。これらの取り扱いにあたっては、合成関数の微分の公式を元に形式的に計算を行う方が暗算を行うよりもミスを減らしやすいと思われる。

ここまでの内容に基づいて以下の問題に答えよ。
i) 下記の関数の微分を計算せよ。
$$
\begin{align}
f(x) &= \log{x} \\
f(x) &= \sqrt{x} \\
f(x) &= \frac{x}{1+x} \\
f(x) &= \frac{1}{1+x} \\
f(x) &= e^x
\end{align}
$$
ⅱ) 合成関数の公式は下記で与えられる。
$$
\begin{align}
\frac{dy}{dx} = \frac{dy}{du} \cdot \frac{du}{dx}
\end{align}
$$
上記の公式を元に下記の関数の微分を考えよ。
$$
\begin{align}
y &= \log{2x} \\
y &= e^{2x} \\
y &= e^{x^2}
\end{align}
$$
ⅲ) 下記のように$y, u, v, r$を設定する。
$$
\begin{align}
y &= \log{u} \\
u &= \frac{v}{1+v} \\
v &= \exp(r) \\
r &= b_0+b_1x_i
\end{align}
$$
このとき$\displaystyle \frac{dy}{du}, \frac{du}{dv}, \frac{dv}{dr}, \frac{dr}{db_0}$をそれぞれ求めよ。
iv) ⅲ)の結果に基づいて$\displaystyle \frac{dy}{d b_0}$を計算せよ。
v) ⅲ)とiv)の結果と同様に下記の関数を$b_0$について微分せよ。
$$
\begin{align}
y = \log{ \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) }
\end{align}
$$
vi) iv), v)の計算に基づいて、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_0}$を求めよ。
vⅱ) iv)、v)でそれぞれ求めた$\displaystyle \frac{dy}{d b_0}$に対する$\displaystyle \frac{dy}{d b_1}$の関係式が$\displaystyle \frac{dy}{d b_1} = x_i \frac{dy}{d b_0}$のように表せることを利用して、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_1}$を求めよ。

・解答

i)
それぞれ各関数の微分の公式や商の導関数の公式を用いることで下記のように導出することができる。
$$
\large
\begin{align}
\left( \log{x} \right)’ &= \frac{1}{x} \\
\left( \sqrt{x} \right)’ &= \frac{1}{2\sqrt{x}} \\
\left( \frac{x}{1+x} \right)’ &= \frac{x'(1+x)-x(1+x)’}{(1+x)^2} \\
&= \frac{1}{(1+x)^2} \\
\left( \frac{1}{1+x} \right)’ &= \frac{(1)'(1+x)-1 \cdot (1+x)’}{(1+x)^2} \\
&= -\frac{1}{(1+x)^2} \\
\left( e^x \right)’ &= e^x
\end{align}
$$

ⅱ)
$$
\large
\begin{align}
y &= \log{2x} \\
y &= e^{2x} \\
y &= e^{x^2}
\end{align}
$$
上記に対してそれぞれ、$u=2x, u=2x, u=x^2$とおくことで、合成関数の微分の公式を用いると下記のように計算できる。
$$
\large
\begin{align}
\frac{dy}{dx} &= \frac{dy}{du} \cdot \frac{du}{dx} \\
&= \frac{1}{u} \cdot 2 \\
&= \frac{2}{2x} \\
&= \frac{1}{x} \\
\frac{dy}{dx} &= \frac{dy}{du} \cdot \frac{du}{dx} \\
&= e^{2x} \cdot 2 \\
&= 2 e^{2x} \\
\frac{dy}{dx} &= \frac{dy}{du} \cdot \frac{du}{dx} \\
&= e^{2x} \cdot (2x) \\
&= 2x e^{2x}
\end{align}
$$

ⅲ)
i)の結果を参考にすることで、$\displaystyle \frac{dy}{du}, \frac{du}{dv}, \frac{dv}{dr}, \frac{dr}{d b_0}$はそれぞれ下記のように求めることができる。
$$
\large
\begin{align}
\frac{dy}{du} &= \frac{1}{u} \\
\frac{du}{dv} &= \frac{1}{(1+v)^2} \\
\frac{dv}{dr} &= \exp(r) \\
\frac{dr}{d b_0} &= 1
\end{align}
$$

iv)
合成関数の微分の公式にⅲ)の結果を適用することで、$\displaystyle \frac{dy}{d b_0}$は下記のように求めることができる。
$$
\large
\begin{align}
\frac{dy}{d b_0} &= \frac{dy}{du} \cdot \frac{du}{dv} \cdot \frac{dv}{dr} \cdot \frac{dr}{d b_0} \\
&= \frac{1}{u} \cdot \frac{1}{(1+v)^2} \cdot \exp(r) \cdot 1 \\
&= \frac{1+\exp(b_0+b_1x_i)}{\cancel{\exp(b_0+b_1x_i)}} \cdot \frac{1}{(1+\exp(b_0+b_1x_i))^2} \cdot \cancel{\exp(b_0+b_1x_i)} \cdot 1 \\
&= \frac{1}{1+\exp(b_0+b_1x_i)}
\end{align}
$$

v)
下記のように$y, u, v, r$を設定する。
$$
\large
\begin{align}
y &= \log{u} \\
u &= \frac{1}{1+v} \\
v &= \exp(r) \\
r &= b_0+b_1x_i
\end{align}
$$

i)の結果を参考にすることで、$\displaystyle \frac{dy}{du}, \frac{du}{dv}, \frac{dv}{dr}, \frac{dr}{d b_0}$はそれぞれ下記のように求めることができる。
$$
\large
\begin{align}
\frac{dy}{du} &= \frac{1}{u} \\
\frac{du}{dv} &= -\frac{1}{(1+v)^2} \\
\frac{dv}{dr} &= \exp(r) \\
\frac{dr}{d b_0} &= 1
\end{align}
$$
合成関数の微分の公式に上記を適用することで、$\displaystyle \frac{dy}{d b_0}$は下記のように求めることができる。
$$
\large
\begin{align}
\frac{dy}{d b_0} &= \frac{dy}{du} \cdot \frac{du}{dv} \cdot \frac{dv}{dr} \cdot \frac{dr}{d b_0} \\
&= \frac{1}{u} \cdot \left( -\frac{1}{(1+v)^2} \right) \cdot \exp(r) \cdot 1 \\
&= – \frac{1+\exp(b_0+b_1x_i)}{1} \cdot \frac{1}{(1+\exp(b_0+b_1x_i))^2} \cdot \exp(b_0+b_1x_i) \cdot 1 \\
&= – \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)}
\end{align}
$$

iv)の結果を$f_1$、v)の結果を$f_2$とするとき、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_0}$は下記のように表すことができる。
$$
\large
\begin{align}
\frac{d \log{L(b_0,b_1)}}{d b_0} &= \frac{d}{d b_0} \sum_{i=1}^{n} \left( y_i \log \left( \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) + (1-y_i) \log{ \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) } \right) \\
&= \sum_{i=1}^{n} \left( y_i f_1 + (1-y_i)f_2 \right)
\end{align}
$$
$f_1, f_2$にiv)、v)の結果を代入することで、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_0}$は下記のように計算できる。
$$
\large
\begin{align}
\frac{d \log{L(b_0,b_1)}}{d b_0} &= \sum_{i=1}^{n} \left( y_i f_1 + (1-y_i)f_2 \right) \\
&= \sum_{i=1}^{n} \left( y_i \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) + (1-y_i) \left(- \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) \right) \\
&= \sum_{i=1}^{n} \left( y_i \left( \frac{1+\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) – \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) \\
&= \sum_{i=1}^{n} \left( y_i – \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) \\
&= \sum_{i=1}^{n} (y_i – p_i)
\end{align}
$$

vⅱ)
iv)の結果を$f_1$、v)の結果を$f_2$とするとき、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_1}$は下記のように表すことができる。
$$
\large
\begin{align}
\frac{d \log{L(b_0,b_1)}}{d b_1} &= \frac{d}{d b_0} \sum_{i=1}^{n} \left( y_i \log \left( \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i)} \right) + (1-y_i) \log{ \left( \frac{1}{1+\exp(b_0+b_1x_i)} \right) } \right) \\
&= \sum_{i=1}^{n} \left( y_i f_1 x_i + (1-y_i)f_2 x_i \right)
\end{align}
$$
$f_1, f_2$にiv)、v)の結果を代入しvi)と同様の計算を行うことで、$\displaystyle \frac{d \log{L(b_0,b_1)}}{d b_0}$は下記のように計算できる。
$$
\large
\begin{align}
\frac{d \log{L(b_0,b_1)}}{d b_1} &= \sum_{i=1}^{n} \left( y_i f_1 x_i + (1-y_i)f_2 x_i \right) \\
&= \sum_{i=1}^{n} \left( y_i – \frac{\exp(b_0+b_1x_i)}{1+\exp(b_0+b_1x_i) } \right) x_i \\
&= \sum_{i=1}^{n} (y_i – p_i) x_i
\end{align}
$$

・解説
vi)、vⅱ)で導出した勾配を元に勾配降下法などを用いることで、ロジスティック回帰のパラメータ推定を行うことができます。ロジスティック回帰のパラメータ推定の導出は途中計算の式がポアソン回帰に比べ複雑であることは抑えておくと良いと思います。

ポアソン回帰の対数尤度と勾配の導出

・問題
$y_i$が離散値の観測値$(x_i, y_i)$が与えられ、$x_i$より$y_i$を予測することを考える。この際によく用いられるのがポアソン分布に基づくGLMのポアソン回帰である。
$$
\begin{align}
P(y_i|\lambda_i) &= \frac{\lambda^{y_i} e^{\lambda_i}}{y_i!} \\
\log{\lambda_i} &= \beta_0 + \beta_1 x_i
\end{align}
$$
ポアソン回帰では上記のような式に基づいて、$\lambda_i$を計算し、この値を予測値に考える。このとき回帰式のパラメータ$\beta_0, \beta_1$の推定にあたっては、ポアソン分布$P(y_i|\lambda_i)$の確率を元に尤度を考え、最尤法を用いる。

ここまでの内容に基づいて、以下の問いに答えよ。
i) $x = \exp(\log(x))$が成立することを利用して、下記が成立することを確認せよ。
$$
\begin{align}
\frac{\lambda^{y_i} e^{\lambda_i}}{y_i!} = \exp( y_i \log{\lambda_i} – \lambda_i – \log{y_i!})
\end{align}
$$
ⅱ) 同時確率$P(y_1,…,y_n|\lambda_1,…,\lambda_n)$を計算せよ。
ⅲ) ⅱ)で計算した同時確率を元にパラメータの推定にあたっての尤度が定義される。ここで$\lambda_i = \beta_0 + \beta_1x_i$をⅱ)の結果に代入することで尤度関数$L(\beta_0,\beta_1)$を計算せよ。
iv) ⅲ)の結果に基づいて対数尤度$\log{L(\beta_0,\beta_1)}$を計算せよ。
v) 対数尤度$\log{L(\beta_0,\beta_1)}$を$\beta_0$で偏微分せよ。
vi) 対数尤度$\log{L(\beta_0,\beta_1)}$を$\beta_1$で偏微分せよ。
vⅱ) v)、vi)の結果を元に勾配法の実行を行うにあたって、注意すべき点に関して答えよ。

・解答
i)
下記のように導出を行うことができる。
$$
\large
\begin{align}
\frac{\lambda^{y_i} e^{\lambda_i}}{y_i!} &= \exp \left( \log \left( \frac{ \lambda_i^{y_i} e^{-\lambda_i}}{y_i!} \right) \right) \\
&= \exp( y_i \log{\lambda_i} – \lambda_i – \log{y_i!})
\end{align}
$$

ⅱ)
同時確率$P(y_1,…,y_n|\lambda_1,…,\lambda_n)$は下記のように計算できる。
$$
\large
\begin{align}
P(y_1,…,y_n|\lambda_1,…,\lambda_n) &= \prod_{i=1}^{n} P(y_i|\lambda_i) \\
&= \prod_{i=1}^{n} \exp( y_i \log{\lambda_i} – \lambda_i – \log{y_i!}) \\
&= \exp( y_1 \log{\lambda_1} – \lambda_1 – \log{y_1!}) … \exp( y_n \log{\lambda_n} – \lambda_n – \log{y_n!}) \\
&= \exp( y_1 \log{\lambda_1} – \lambda_1 – \log{y_1!} + … + y_n \log{\lambda_n} – \lambda_n – \log{y_n!}) \\
&= \exp \left( \sum_{i=1}^{n} (y_i \log{\lambda_i} – \lambda_i – \log{y_i!}) \right)
\end{align}
$$

ⅲ)
ⅱ)の結果に$\lambda_i = \beta_0 + \beta_1x_i$をⅱ)を代入することで尤度関数$L(\beta_0,\beta_1)$は下記のように計算できる。
$$
\large
\begin{align}
L(\beta_0,\beta_1) &= \exp \left( \sum_{i=1}^{n} (y_i \log{\lambda_i} – \lambda_i – \log{y_i!}) \right) \\
&= \exp \left( \sum_{i=1}^{n} (y_i \log{(\exp(\beta_0 + \beta_1x_i))} – \exp(\beta_0 + \beta_1x_i) – \log{y_i!}) \right) \\
&= \exp \left( \sum_{i=1}^{n} (y_i (\beta_0 + \beta_1x_i) – \exp(\beta_0 + \beta_1x_i) – \log{y_i!}) \right)
\end{align}
$$

iv)
対数尤度$\log{L(\beta_0,\beta_1)}$は下記のようになる。
$$
\large
\begin{align}
\log{L(\beta_0,\beta_1)} &= \log{ \left( \exp \left( \sum_{i=1}^{n} (y_i (\beta_0 + \beta_1x_i) – \exp(\beta_0 + \beta_1x_i) – \log{y_i!}) \right) \right) } \\
&= \sum_{i=1}^{n} (y_i (\beta_0 + \beta_1x_i) – \exp(\beta_0 + \beta_1x_i) – \log{y_i!})
\end{align}
$$

v)
下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial \log{L(\beta_0,\beta_1)}}{\partial \beta_0} &= \frac{\partial}{\partial \beta_0} \sum_{i=1}^{n} (y_i (\beta_0 + \beta_1x_i) – \exp(\beta_0 + \beta_1x_i) – \log{y_i!}) \\
&= \sum_{i=1}^{n} (y_i – \exp(\beta_0 + \beta_1x_i))
\end{align}
$$

vi)
下記のように計算できる。
$$
\large
\begin{align}
\frac{\partial \log{L(\beta_0,\beta_1)}}{\partial \beta_1} &= \frac{\partial}{\partial \beta_1} \sum_{i=1}^{n} (y_i (\beta_0 + \beta_1x_i) – \exp(\beta_0 + \beta_1x_i) – \log{y_i!}) \\
&= \sum_{i=1}^{n} (y_i x_i – \exp(\beta_0 + \beta_1x_i)x_i) \\
&= \sum_{i=1}^{n} (y_i – \exp(\beta_0 + \beta_1x_i)) x_i
\end{align}
$$

vⅱ)
式に指数関数が含まれていることから、パラメータの値によってはオーバーフローを起こす可能性があることに注意が必要である。

・解説
下記でPythonを用いた具体的なプログラミングを取り扱ったので、合わせて確認するとわかりやすいと思います。
Pythonを用いたポアソン回帰のプログラミング

参考書籍

・Pythonで学ぶ統計モデリング（運営者作成）
https://www.amazon.co.jp/dp/B08FYMTYBW/