【JDLA E資格】最尤推定法

#はじめに
JDLA E資格試験で出題される最尤推定法について解説した記事です。
E資格試験の応用数学パートでは、最尤推定法を用いて、確率分布のパラメータを推定する問題が出題されます。
特に、正規分布やベルヌーイ分布のパラメータ推定は頻出です。
また、最尤推定法も機械学習・深層学習全般において使用されるため、本稿の内容を理解しておくのは必須です。

なお、他パートの具体的な解説については、下記をご覧ください。
[E資格試験に関する私の投稿記事リスト][link-1]
[link-1]:https://qiita.com/fridericusgauss/items/5a97f2645cdcefe15ce0

###目次

1. 尤度
2. 最尤推定法
3. 最尤推定法の適用
4. おわりに

###数学表記
$\boldsymbol{0}$は零ベクトル、$\mathbb{R}$は実数集合です。$\Omega$は確率変数$X$が取り得る標本空間、$\Theta$は確率分布のパラメータ$\theta$が取り得るパラメータ空間を表します。
また、$f(x; \theta)$と表記したとき、$;$の前の$x$を関数$f$の変数、$;$の後ろの$\theta$を関数$f$のパラメータとして扱います。
確率や確率分布の諸定義については、下記をご覧ください。
[期待値・分散][link-2]
[確率分布][link-3]
[link-2]:https://qiita.com/fridericusgauss/items/ab36d7cf91093d284ba0
[link-3]:https://qiita.com/fridericusgauss/items/14bbcdc0423fc39067c2

#尤度
###尤度
最尤推定法は、__サンプリングされたデータ（標本）を観測した後、その標本は元々どのようなパラメータを有する確率分布から生成されたのか？__という問いについて答えるための一つの方法です。
つまり、「母集団を直接調査するのを諦めて、代わりに母集団から抽出した標本を調査し、その標本を用いて母集団の統計量を推定する」という標本理論の考えに従っています。
[標本理論][link-3]

最尤推定法では尤度を用います。__尤度（ゆうど）は、尤もらしい度合（あり得る度合）を表す概念__です。確率もある得る度合を表すため、尤度の一種として考えられます。尤度は確率よりも少し広い概念で、立場が少し異なりますが、確率と同じ扱いをしても計算上問題ないことが多いです。

実用的な場面で強調される主な違いは下記の点です。

• 確率密度関数はパラメータが与えられた（固定した）とき、確率変数を変数とした（動かした）ときのあり得る度合を表す。
• 尤度は確率変数がデータから与えられた（固定した）とき、パラメータを変数とした（動かした）ときのあり得る度合を表す。

この違いに、__標本（データ集合）を用いて母集団が従う確率分布のパラメータを推定する__考えが反映されています。
尤度関数は、尤度を関数として数式表現したもので、本質的には尤度と同義です。

なお、本稿では、尤度の説明をこの程度に留めておき、最尤推定法の解説に注力します。
例えば、尤度の視覚的イメージは下記などを参考にしてください。
[【統計学】尤度って何？をグラフィカルに説明してみる。][link-4]
[link-4]:https://qiita.com/kenmatsu4/items/b28d1b3b3d291d0cc698

###尤度関数
確率分布のパラメータを$\theta \in \Theta$、それに従う確率変数を$X\in \Omega$としたとき、その確率密度関数（確率質量関数）を$f(X=x;\theta)=f(x;\theta)$と表します。
また、データ$x$を与えたとき、それに対応する尤度関数を$L(\theta; X=x)=L(\theta; x)$と表し、式(1)とします。

L(\theta; x)=f(x;\theta)
\tag{1}

データ集合$D=\{x_1,x_2,\cdots,x_N\}$を与えたとき、その尤度関数$L(\theta;D)$は式(2)で表されます。

\begin{array}{ll}
L(\theta; D) &= L(\theta; x_1,x_2,\cdots,x_N)\\
 &= \prod_{i=1}^{N} L(\theta; x_{i}) \\
 &=\prod_{i=1}^{N} f(x_{i};\theta)
\end{array}
\tag{2}

式(2)は、各データを同時に与えたときの確率分布（同時分布）なので、各データに対する尤度関数（確率分布）の積としています。

式(2)の負の対数尤度関数$L_{\log}(\theta; D)$は式(3)で表されます。

\begin{array}{ll}
L_{\log}(\theta; D)
 &= -\ln L(\theta; D) \\
&= -\ln \prod_{i=1}^{N} L(\theta; x_{i})\\
&= - \sum_{i=1}^{N} \ln L(\theta; x_{i})\\
&= - \sum_{i=1}^{N} \ln f(x_{i};\theta)
\end{array}
\tag{3}

ただし、$\ln x$は自然対数$\log_{\mathrm{e}} x$です。

実際の最尤推定法では、負の対数尤度関数を用いるケースがほとんどです。
その理由は下記の通りです。

• 尤度関数の最大化を行いたいが、一般の最適化問題は最小化問題であるため、その形式に整合するために、$-1$をかけた。
• 式(3)の通り、尤度関数は確率の積の形式で表されるが、計算の都合上、和の形式のほうが容易（例えば指数関数）なので、自然対数をとった。
• $[0,1]$以内の小数点である確率を多数乗算すると、数値計算上、浮動小数点演算の精度を失う可能性があるため、自然対数をとった。

なお、$x$が単調増加するとき、負の自然対数$-\ln x$は単調減少する関係です。

#最尤推定法
###方針
最尤推定法とは、__標本（データ集合）から得られる尤度に基づいて、母集団が従う確率分布のパラメータを推定する方法__です。
これは、標本理論で説明したように、母集団が従う確率分布のパラメータは不明ですが、そこから得た一部の標本から、最も上手く説明できるパラメータを求めることになります。
「最も上手く説明できる」というのが、「標本から判断したときの尤もらしい度合（尤度）が最大となる」と捉えています。

つまり、最尤推定法とは、尤度の最大化であり、それを達成するパラメータを求めます。
形式的には、式(4)に示すように、尤度関数$L(\theta; D)$を最大化する$\hat{\theta}$を求めます。

\hat{\theta}=\arg \max_{\theta\in \Theta} L(\theta; D)
\tag{4}

負の対数尤度関数$L_{\log}(\theta; D)$を用いると、$L_{\log}(\theta; D)$を最小化する$\hat{\theta}$を求める問題となるため、式(5)で表されます。

\hat{\theta}=\arg \min_{\theta\in \Theta} L_{\log}(\theta; D)
\tag{5}

この過程で求めた$\hat{\theta}$を__最尤推定量__と呼びます。

さらに、この考え方は、多くの機械学習においても共通します。
機械学習では、タスクを達成するための数理モデルを仮定し、与えたデータ集合にモデルが適合するように、モデルのパラメータの値を決定します。
その方法として、データとモデルとの間の誤差を表現した__誤差関数__を定義し、その誤差関数を最小化するようなパラメータを求めます。
つまり、機械学習における学習とは、誤差の最小化であり、それを達成するパラメータを求めます。
このため、他の機械学習においても、いくつかの尤度関数を誤差関数として利用する場合があります。

###計算方法
式(5)は最適化変数を$\theta$としたとき、目的関数$L_{\log}(\theta; D)$を最小化する最適化問題です。
このため、基本的には$L_{\log}(\theta; D)$を$\theta$について偏微分したものが$0$となる方程式を満たす$\theta$が最尤推定量$\hat{\theta}$となります。
この方程式を__尤度方程式__と呼び、式(6)で表されます。

\frac{\partial L_{\log}(\theta; D)}{\partial \theta}=0
\tag{6}

また、最適化変数がスカラ$\theta$ではなく、ベクトル$\boldsymbol{\theta}=(\theta_{1},\theta_{2},\cdots,\theta_{d})^{\mathrm{T}}$の場合、尤度方程式は下記のように表されます。

\frac{\partial L_{\log}(\boldsymbol{\theta}; D)}{\partial \boldsymbol{\theta}}=\boldsymbol{0}

等式の制約条件が課されているときはラグランジュの未定乗数法などで求めます。
しかし、解析的に解けない場合、あるいは、目的関数$L_{\log}(\theta; D)$の景観が明らかでない場合、尤度方程式（式(6)）を直接解くことは一般に困難です。
このとき、最適化アルゴリズムを用いて最適化問題を解くことで、数値計算的に最尤推定量を求めます。
例えば、ニューラルネットワークの場合、モデルパラメータが多く存在し、目的関数の景観が明らかでないため、誤差逆伝播法や最適化アルゴリズムを用いることになります。

#最尤推定法の適用
###ベルヌーイ分布への適用
最尤推定法を用いて、ベルヌーイ分布のパラメータをデータ集合から推定してみます。
ベルヌーイ分布の確率質量関数$f(x)$は式(7)で表されます。

f(x; p)=p^{x}(1-p)^{1-x}
\tag{7}

ただし、標本空間は$\Omega=\{0,1\}$、$p\in [0,1]$は確率パラメータであり、式(8)の関係を満たしています。

\left\{
\begin{array}{ll}
f(0)=1-p \\
f(1)=p
\end{array}
\right.
\tag{8}

以下では、データ集合$D=\{x_1,x_2,\cdots,x_N\}$から、ベルヌーイ分布の確率パラメータ$p$を推定します。
式(2)から、尤度関数$L(p;D)$は式(9)で表されます。

\begin{array}{ll}
L(p; D) &= \prod_{i=1}^{N} L(p; x_{i})\\ 
&=\prod_{i=1}^{N} p^{x_{i}}(1-p)^{1-x_{i}}
\end{array}
\tag{9}

式(9)の負の対数尤度関数$L_{\log}(p; D)=L_{\log}(p)$は式(10)で表されます。

\begin{array}{ll}
L_{\log}(p) &= -\ln L(p; D)\\
&= - \sum_{i=1}^{N} \ln (p^{x_{i}}(1-p)^{1-x_{i}})\\
&= - \sum_{i=1}^{N} (\ln p^{x_{i}} + \ln (1-p)^{1-x_{i}})\\ 
&= - \sum_{i=1}^{N} (x_{i}\ln p + (1-x_{i}) \ln (1-p))
\end{array}
\tag{10}

この尤度方程式は式(11)で表されます。

\frac{d L_{\log}(p)}{d p}=0
\tag{11}

式(11)は下記のように変形できます。

\begin{array}{ll}
\frac{d L_{\log}(p)}{d p}
&= - \sum_{i=1}^{N} \left(\frac{x_{i}}{p} - \frac{1-x_{i}}{1-p}\right)\\
&= - \frac{1}{p(1-p)} \sum_{i=1}^{N} \left((1-p)x_{i} - p(1-x_{i})\right)\\
&= - \frac{1}{p(1-p)} \sum_{i=1}^{N} \left(x_{i}-px_{i} - p + px_{i}\right)\\
&= - \frac{1}{p(1-p)} \sum_{i=1}^{N} (x_{i}- p)\\
&= - \frac{1}{p(1-p)} \left(\sum_{i=1}^{N} x_{i}- \sum_{i=1}^{N}p\right)\\
&= -\frac{1}{p(1-p)}\left( \sum_{i=1}^{N} x_{i} - pN\right) = 0
\end{array}

よって、下記のように最尤推定量$\hat{p}$が求められます。

\begin{array}{ll}
\sum_{i=1}^{N} x_{i} - \hat{p}N = 0 
&\Leftrightarrow& \hat{p} = \frac{1}{N}\sum_{i=1}^{N} x_{i}
\end{array}

なお、これは__最尤推定量$\hat{p}$が期待値$\mathbb{E}[X=x]$に一致する__ことを示しています。
ベルヌーイ試行（コイントス）のデータが多く存在するとき、最尤推定量$\hat{p}$は母平均$1/2$に近づきます。
[ベルヌーイ分布の期待値][link-3]

###正規分布への適用
最尤推定法を用いて、正規分布のパラメータをデータ集合から推定してみます。
正規分布の確率密度関数$f(x)$は式(12)で表されます。

f(x; \mu ,\sigma^2)=\frac{1}{\sqrt {2\pi \sigma^2}}\exp{\biggl (}-{\frac {(x-\mu)^2}{2\sigma^{2} }}{\biggr )}
\tag{12}

ただし、標本空間は$\Omega\subseteq\mathbb{R}$、$\mu\in\mathbb{R},\sigma^2>0$は平均パラメータ、分散パラメータです。

以下では、データ集合$D=\{x_1,x_2,\cdots,x_N\}$から、ベルヌーイ分布の平均パラメータ$\mu$と、分散パラメータ$\sigma^2$を推定します。
式(2)から、尤度関数$L(\mu ,\sigma^2; D)$は式(13)で表されます。

\begin{array}{ll}
L(\mu ,\sigma^2; D) 
&= \prod_{i=1}^{N} L(\mu ,\sigma^2; x_{i})\\ 
&= \prod_{i=1}^{N} \frac{1}{\sqrt {2\pi \sigma^2}}\exp{\left(-{\frac{(x_{i}-\mu)^2}{2\sigma^{2}}}\right)}
\end{array}
\tag{13}

式(13)の負の対数尤度関数$L_{\log}(\mu ,\sigma^2; D)=L_{\log}(\mu ,\sigma^2)$は式(14)で表されます。

\begin{array}{ll}
L_{\log}(\mu ,\sigma^2) 
&= -\ln L(\mu ,\sigma^2; D)\\
&= - \sum_{i=1}^{N} \ln \biggl(\frac{1}{\sqrt{2\pi \sigma^2}}\exp{\left(-\frac{(x_{i}-\mu)^2}{2\sigma^2}\right)}\biggr)\\
&= - \sum_{i=1}^{N} \biggl(\ln \frac{1}{\sqrt{2\pi \sigma^2}} + \ln \exp{\left(-\frac{(x_{i}-\mu)^2}{2\sigma^2}\right)}\biggr)\\ 
&= - \sum_{i=1}^{N} \biggl(\ln \frac{1}{\sqrt{2\pi \sigma^2}} -\frac{1}{2\sigma^2}(x_{i}-\mu)^2\biggr)\\ 
&= \sum_{i=1}^{N} \biggl(\frac{1}{2}\ln (2\pi \sigma^2) +\frac{1}{2\sigma^2}(x_{i}-\mu)^2\biggr)\\
&= \frac{1}{2}\sum_{i=1}^{N} \biggl(\ln (2\pi \sigma^2) +\frac{1}{\sigma^2}(x_{i}-\mu)^2\biggr)
\end{array}
\tag{14}

$\mu$に関する尤度方程式は式(15)で表されます。

\frac{\partial L_{\log}(\mu,\sigma^2)}{\partial \mu}=0
\tag{15}

式(15)は下記のように変形できます。

\begin{array}{ll}
\frac{\partial L_{\log}(\mu,\sigma^2)}{\partial \mu}
&= -\frac{1}{\sigma^2}\sum_{i=1}^{N} (x_{i}-\mu)\\
&= -\frac{1}{\sigma^2}\left(\sum_{i=1}^{N} x_{i}- \sum_{i=1}^{N}\mu\right)\\
&= -\frac{1}{\sigma^2}\left(\sum_{i=1}^{N} x_{i} - N\mu\right) = 0
\end{array}

よって、下記のように最尤推定量$\hat{\mu}$が求められます。

\begin{array}{ll}
&\sum_{i=1}^{N} x_{i} - N\hat{\mu} = 0 \\
\Leftrightarrow& \hat{\mu} = \frac{1}{N}\sum_{i=1}^{N} x_{i}
\end{array}

なお、これは__最尤推定量$\hat{\mu}$が期待値$\mathbb{E}[X=x]$に一致する__ことを示しています。

また、$\sigma^2$に関する尤度方程式は式(16)で表されます。

\frac{\partial L_{\log}(\mu,\sigma^2)}{\partial \sigma^2}=0
\tag{16}

式(16)は下記のように変形できます。

\begin{array}{ll}
\frac{\partial L_{\log}(\mu,\sigma^2)}{\partial \sigma^2}
&= \frac{1}{2}\sum_{i=1}^{N} \biggl(2\pi\frac{1}{2\pi \sigma^2} -\frac{1}{(\sigma^2)^2}(x_{i}-\mu)^2\biggr)\\
&= \frac{1}{2(\sigma^2)^2}\sum_{i=1}^{N} \biggl(\sigma^2 -(x_{i}-\mu)^2\biggr)\\
&= \frac{1}{2(\sigma^2)^2}\left(\sum_{i=1}^{N} \sigma^2 - \sum_{i=1}^{N}(x_{i}-\mu)^2\right)\\
&= \frac{1}{2(\sigma^2)^2}\left(\sigma^2 N - \sum_{i=1}^{N}(x_{i}-\mu)^2\right)\\
&= 0
\end{array}

よって、下記のように最尤推定量${\hat{\sigma}}^2$が求められます。

\begin{array}{ll}
&{\hat{\sigma}}^2 N - \sum_{i=1}^{N}(x_{i}-\hat{\mu})^2 = 0 \\
\Leftrightarrow& {\hat{\sigma}}^2 = \frac{1}{N}\sum_{i=1}^{N} (x_{i}-\hat{\mu})^2
\end{array}

なお、これは__最尤推定量${\hat{\sigma}}^2$が分散$\mathbb{V}[X=x]$に一致する__ことを示しています。

#おわりに
E資格向けの最尤推定法について解説しました。
なお、上記は、2021年2月時点における内容であることにご注意ください。

[E資格試験に関する私の投稿記事リスト][link-1]