尤度・最尤推定

コインの例

最尤推定法(maximum likelihood estimation)は確率分布 $f_D$ のパラメーター $\theta$ を、その母集団から得られたサンプルから推定する方法。

たとえば、表の出る確率が必ずしも1/2でないコインを50回投げた時、表が15回出たとすると、このコインの表が出る確率はいくらか、といった問題。20～30回くらいの間だと普通のコインで表／裏の確率は1/2かなと思うが、50回投げて表が15回しかでないとちょっと疑いたくなる。

ここで、表が出る確率を $p$ として、 $n$ 回の試行中表が $k$ 回出る確率は以下の通り。

(1) $\begin{equation*} P(k:n\,|\,p) = \binom{n}{k} p^k (1-p)^{n-k} \end{equation*}$

先の問題の結果について、母集団の確率をいくつか仮定して計算してみると

(2) $\begin{equation*} \begin{align} p &= 1/2 \quad \rightarrow \quad P(15:50\,|\,0.5) \approx 0.002 \\ p &= 1/3 \quad \rightarrow \quad P(15:50\,|\,0.3) \approx 0.107 \\ p &= 1/4 \quad \rightarrow \quad P(15:50\,|\,0.1) \approx 0.089 \end{align} \end{equation*}$

得られたサンプルの下では、少なくともぼ表が出る確率を1/2や1/4と考えるよりは、1/3と考えた方がこのサンプルのパターンが発生する確率が高い。

そこで、サンプルのパターンに対して母集団の確率を変化させてみると、以下のようになる。

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import binom

n = 50
k = 15
p = np.linspace(0, 1, 100)

fig, ax = plt.subplots(figsize=(4.8, 3.6))

ax.plot(p, binom.pmf(k, n, p), c='navy')
ax.plot((k/n, k/n), (0, binom.pmf(k, n, k/n)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)
ax.set_ylim(0, 0.15)
ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from scipy.stats import binom

n = 50

k = 15

p = np.linspace(0, 1, 100)

fig, ax = plt.subplots(figsize=(4.8, 3.6))

ax.plot(p, binom.pmf(k, n, p), c='navy')

ax.plot((k/n, k/n), (0, binom.pmf(k, n, k/n)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)

ax.set_ylim(0, 0.15)

ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

グラフ上は、母集団確率が15/50=0.3のときにサンプルパターンの発生確率が最も高くなっている。

解析的に確率が最大となる母集団確率を求めるために、式(1)を $p$ で微分してゼロとなる点を求める。

(3) $\begin{gather*} \frac{d P(k:n\,|\,p)}{dp} = \binom{n}{k}\left(kp^{k-1} (1-p)^{n-k} - p^k (n-k)(1-p)^{n-k-1} \right) = 0 \\ p^{k-1} (1-p)^{n-k-1} \left( k (1-p) - (n-k) p \right) = 0 \\ p^{k-1} (1-p)^{n-k-1} \left( k - np \right) = 0 \end{gather*}$

上式の解は $p = 0, 1, k/n$ となるが、 $p=0, 1$ の場合は式(1)がゼロとなることから、確率が最大となるのは $p=k/n$ であることがわかる。

尤度関数と最尤法

母集団の確率分布を $f$ 、母集団のパラメータを $\theta$ 、得られたサンプルセットを $x_1, \ldots, x_n$ とすると、パラメーターが $\theta$ である確率は

(4) $\begin{equation*} Pr(x_1, \ldots, x_n \,|\, \theta) = f(x_1, \ldots, x_n \,|\, \theta) \end{equation*}$

このとき、パラメーター $\theta$ に関する尤度関数は以下のように定義される。

(5) $\begin{equation*} L(\theta) = f(x_1, \ldots, x_n \,|\, \theta) \end{equation*}$

通常、尤度関数が最大となるパラメーターを求めるためには対数尤度関数が用いられる。

(6) $\begin{equation*} \frac{d}{d \theta} \ln L(\theta) = 0 \end{equation*}$

サンプルの質と尤度

先のコインの問題で、試行数が10回の場合と100回の場合を比べてみる。それぞれについて、表の確率が0.3、0.5、0.7に対応する回数は試行回数が10回の場合は3回、5回、7回、試行回数100回に対しては30回、50回、70回である。それぞれの表の階数に対して、母集団確率の尤度関数を描くと以下のようになる。

import numpy as np
import matplotlib.pyplot as plt
from functools import reduce

def comb(n, k):
    numer = reduce(lambda x,y:x*y, range(n - k + 1, n))
    denom = reduce(lambda x,y:x*y, range(1, k + 1))
    return numer // denom

def binom(n, k, p):
    return comb(n, k) * p**k * (1 - p)**(n - k)

n = [10, 100]
k = [np.array([3, 5, 7]), np.array([30, 50, 70])]

p = np.linspace(0, 1, 100)

fig, axs = plt.subplots(1, 2, figsize=(9.6, 3.6))

for ax, n, k in zip(axs, n, k):
    ax.set_title("n={}".format(n))

    for t in k:
        p_mle = t / n
        prob = binom(n, t, p)

        ax.plot(p, prob, c='navy')
        ax.plot((p_mle, p_mle), (0, binom(n, t, p_mle)), c='gray', linestyle='dashed')

        ax.set_xlim(0, 1)
        ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

import numpy as np

import matplotlib.pyplot as plt

from functools import reduce

def comb(n, k):

numer = reduce(lambda x,y:x*y, range(n - k + 1, n))

denom = reduce(lambda x,y:x*y, range(1, k + 1))

return numer // denom

def binom(n, k, p):

return comb(n, k) * p**k * (1 - p)**(n - k)

n = [10, 100]

k = [np.array([3, 5, 7]), np.array([30, 50, 70])]

p = np.linspace(0, 1, 100)

fig, axs = plt.subplots(1, 2, figsize=(9.6, 3.6))

for ax, n, k in zip(axs, n, k):

ax.set_title("n={}".format(n))

for t in k:

p_mle = t / n

prob = binom(n, t, p)

ax.plot(p, prob, c='navy')

ax.plot((p_mle, p_mle), (0, binom(n, t, p_mle)), c='gray', linestyle='dashed')

ax.set_xlim(0, 1)

ax.set_xticks(np.arange(0, 1.1, 0.1))

plt.show()

試行回数が100回の場合には分布が尖っており、最尤推定されたパラメーターの確度が高いが、10回の場合には分布がなだらかで、最尤推定された値以外のパラメーターの確率も高くなっている。

最尤推定の場合、サンプルの規模が推定値の信頼度に関わってくる。

TauStation

尤度・最尤推定

コインの例

尤度関数と最尤法

サンプルの質と尤度

コメントを残すコメントをキャンセル

コインの例

尤度関数と最尤法

サンプルの質と尤度

コメントを残す コメントをキャンセル

コメントを残すコメントをキャンセル