院生卒業してからやっと最尤推定を理解した話

はじめに

こんにちは！
GA technologiesのAI Strategy Center（AISC）に所属している王 皓と申します。
現在、AISCでは有志で集まって数学の勉強会を開いております。週一頻度で実施し、決めた教材を輪読していき、疑問点を議論したり問題を解いたりしています。私は大学時代が財務系、院生時代は経済学専門と、基本的には文系でしたから、数学はずっと苦手でした。（ちなみに同じチーム内の三田さんも数学が苦手と宣言していますが、それは嘘です）。データサイエンスの専門性を高めるために不可欠だと思って参加しました。会社員になってからになってから改めて講義のような形で数学を学べるのは懐かしかったです。
最近はその勉強会の参加者で「数学についての情報を発信をしよう！」となっており、今のところ以下の記事がでております。よろしければこちらもご覧いただければ幸いです。

「最尤推定」は恐らく理系の方なら誰でも聞いたことがある一般概念ですが、自分も学生時代に2度も最尤法の授業を受けました。不才な人ですから、総積に対して微分を取ってその極値を求める事しか覚えていなかったです。なぜそうやっているのか？どの目的でやっているのか？などの事は全く理解できませんでした。

一生最尤推定を理解できない可能性もありましたが、今度のチーム育成制度に恵まれてようやく理解しました。

以下の内容は福中先生が作った最尤推定法の資料を参考にして作成します。

尤度（likelihood）とは？

まず尤度の概念を理解しましょう！
尤度は想定するパラメーターがある値をとる場合に観測している事柄や事象が起こりうる確率のことで、簡単にいうと、データ全体の得られやすさを表す確率です。

例えば：
ホームランを打つ確率が30%の野球選手が5回バットを振ったところ、下の表のような結果になったとします。このときの『このデータ全体』の得られる確率はいくらになるでしょうか？

（ホームランを打った時が1、外れは0）
試行    |  データ  | 確率
-------|---------|------
1回目   |   0     |  0.7
2回目   |   0     |  0.7
3回目   |   1     |  0.3
4回目   |   0     |  0.7
5回目   |   1     |  0.3

各データの発生確率は上表の通りなので、データ全体の得られる確率はその総積で求めることができます。

0.7 * 0.7 * 0.3 * 0.7 * 0.3 = 0.03087 (データ全体の得られる確率) 

さらに$${X}$$の文字式で表現してみましょう！

（ホームランを打った時が1、外れは0）
試行    |  データ  | 確率
-------|---------|------
1回目   |   0     |  1-x
2回目   |   0     |  1-x
3回目   |   1     |   x
4回目   |   0     |  1-x
5回目   |   1     |   x

(1- X) * (1- X) * X * (1- X) * X = (1- X)^3 * X^2 = L(X)

この$${L(X)}$$は尤度（『このデータ全体』の得られる確率）を求める尤度関数です。

最尤推定法

上では「外れ、外れ、ホームラン、外れ、ホームラン」の結果を得られる確率を関数$${L(X)}$$で表現しました。ただ$${X}$$の値によって、$${L(X)}$$も変わります。見やすいように、$${L(X)}$$をグラフで表現しましょう！pythonを利用し、描画することができます。

import matplotlib.pyplot as plt
import numpy as np
import math

def Likelihood(a):
    s = a**2 * (1-a)**3
    return s

x = np.arange(0, 1, 0.01)
y = Likelihood(x)

plt.plot(x, y)
plt.show()

その描画の結果はこちらになります。

関数L(x)

図のとおり、$${X}$$が0から1まで変化すると、関数$${L(X)}$$は最初徐々に増加し、0.4を超えたあたりで急速に減少します。山のような形をしています。間違いなく$${L(X)}$$が最大化する点は山のピークになります。この尤度を最大化するような$${X}$$を探索する推定方法は最尤推定法と呼び、この時の$${X}$$が最尤推定量になります。
最尤推定量の求め方に関しては以下の３つの方法があります：

1. 尤度関数を直接微分する方法

2. 尤度関数を対数変換してから微分する方法

3. 数値計算

尤度関数は簡単なら、直接微分法でも解けますが、複雑な関数ならpythonなどを利用する事を推奨します。例えば今度の尤度関数に関して以下のコードで解を求められます。

# 関数定義の際、関数全体をマイナス変換しておくことに注意

from scipy import optimize
def Likelihood(a):
    s = -(a**2 * (1-a)**3)
    return s
x_min = optimize.brent(Likelihood, brack=(0,0.5,1))
print(x_min)
0.3999999941152572

回帰分析への応用

ようやく本番に入りますね。今までの内容は学校でも勉強しましたが、最尤推定法を回帰分析へどう応用するのかはわかりませんでした。
データ解析の目的の一つは母数の推定になります。

データ解析の目的の一つ

今まで勉強した母数の推定法は最小二乗法がmainになり、確率ベースではありませんでした。なので確率ベースでどう母数を推定するのか興味があります。
実はなんと、1つ1つのデータは『確率的』に発生したものだと考えます。例えば上の例は下の表みたいにとある確率で得られたとも考えられます。

x   |  y  | 確率
----|-----|------
3   |  2  |   ?
5   |  7  |   ?
8   |  6  |   ?
9   |  8  |   ?

ただ、各データの得られる確率分布は未知なものだから、尤度関数も作れないです。やっぱり行き詰まりじゃないかと思います。その時に人の知恵に感動しました。
確かに各データの確率分布はわからないが、でもモデルには通常 『誤差』が含まれるので、回帰モデルを正確に記述すると以下のように書けます。ここで$${e}$$は誤差を表しています。

y = b + ax + e

先ほどのデータの各行に対してこのモデルを当てはめ、誤差$${e}$$について解くと以下のようになります。

x   |  y  |  e
----|-----|------
3   |  2  |  2 - b - 3a
5   |  7  |  7 - b - 5a
8   |  6  |  6 - b - 8a
9   |  8  |  8 - b - 9a

ここでは、$${X}$$と$${Y}$$のデータの発生確率はわからないけど、$${e}$$の得られやすさを表す確率はわかるんです！なぜかと言うと、古典的な回帰分析では「誤差$${e}$$はお互いに独立で$${N(0, \sigma^2)}$$の正規分布に従う」という仮定を置いています。ただし今回は簡単のため、$${N(0, 1)}$$の標準正規分布に従うことにします。（ここで感動した笑、ああ確かにそういうことだ）

$$
f(e) = \frac{1}{\sqrt{2\pi}} \exp\left(-\frac{e^2}{2}\right)
$$

それで、各データにおいて、母数$${a}$$と$${b}$$が与えられた場合の各データに関する確率（密度）は以下の表の$${f(e)}$$のようになります。

これでやっと尤度関数を作れます！尤度関数を作れば、残りは解ける問題だけです。人力では解きにくいので、pythonのコードを共有します。

# 多変数関数の最適化にはscipy.optimize.fminを用いる（ネルダー・ミード法）

from scipy import optimize
import numpy as np

def Likelihood(X):
    s = -(
	(np.exp(-0.5*(2-X[1]-3*X[0])**2)/np.sqrt(2*np.pi))*
	(np.exp(-0.5*(7-X[1]-5*X[0])**2)/np.sqrt(2*np.pi))*
	(np.exp(-0.5*(6-X[1]-8*X[0])**2)/np.sqrt(2*np.pi))*
	(np.exp(-0.5*(8-X[1]-9*X[0])**2)/np.sqrt(2*np.pi))
	)
    return s

x_min = optimize.fmin(Likelihood, [1,1])
print(x_min)

解は$${[a=0.758, b=1.011]}$$になります。つまり$${y = 0.758 + 1.011x}$$は今度のモデルになりますね。

誤差項の分布

今回は誤差項$${e}$$が正規分布に従うという仮定のもとで尤度関数を作り、解を求めましたが、統計モデリングにおいて、誤差項の分布は問題の性質やデータの特性に応じて異なります。

例えば、線形回帰モデルでは、測定誤差や観測ノイズなどが正規分布に従うと仮定されることが多いです。ロジスティック回帰モデルでは、誤差項がベルヌーイ分布に従うと仮定されることが多いです。顧客がサービスを利用する間隔やシステムのリクエストが到着する間隔など、イベントの到着間隔をモデル化する際には、指数分布が使用されますので、誤差項が指数分布に従うと仮定されることが多いです

終わり

以上の内容を踏まえて、最尤推定法について理解することができました。数学の理解が苦手だった私たちも、チームの協力と学習の機会によって、この重要な概念を克服することができました。最尤推定法は、データ分析や統計モデリングにおいて不可欠な手法であり、その理解は私たちの専門性を高める上で重要な一歩です。今後も、数学の理解を深め、データサイエンスの領域での成長を目指していきたいと考えています。