PRML 第1章の「最尤推定によるパラメータフィッティング」の解説

何の話かというと

「PRML 第1章の多項式フィッティングの例を再現」では、二乗平均平方根誤差を最小にするという条件でフィッティングを行いました。しかしながら、フィッティングを行う条件は他にもいろいろ考えられます。「有用な結果（つまり、未知のデータをよりよく推定する関数）を得るには、どのような条件でフィッティングするのがよいか」という非常に難しい問題があります。

ここでは、「どれがよいか」という評価方法は別にして、数学的に扱いやすい（数学的な性質の分析が行いやすい）フィッティング方法を2つ紹介します。

尤度関数を用いる方法（最尤推定）
パラメータの事後確率を求める方法（ベイズ推定）

これらは、どちらも「ある確率」を計算してパラメータを決定します。ただし、「ある確率」として何を使用するかが異なります。「どのような確率を使ってどのように推定するのか」も、採用するモデルの一部と考えてください。

使用するトレーニングセットとフィッティング関数は、「PRML 第1章の多項式フィッティングの例を再現」と同じ例を用います。念のために再掲しておきます。

――――
トレーニングセット $\{(x_n, t_n)\}_{n=1}^{N}$ は、次のように生成します。

説明変数 $x_n$ に対して、下記の関数値（正弦波）に正規分布（平均 0、標準偏差 0.3）のノイズを加えたものを目的変数 $t_n$ として生成します。 $x$ の値は $0 \le x \le 1$ の範囲を10個のデータポイントに等分します。

　 $y=\sin(2\pi x)\,\,\,\, (0 \le x \le 1)$

このようにして得られたトレーニングセットをM次多項式でフィッティングします。

　 $y(x,\mathbf{w}) = \sum_{m=0}^{M} w_m x^m$

　 $\mathbf{w} := \begin{pmatrix}w_0 \\ \vdots \\ w_M\end{pmatrix}$
――――

ここでは、フィッティング対象のM次多項式を明示的に $y(x,\mathbf{w})$ と表現しています。また、この後の説明のために、平均 $\mu$ 、分散 $\sigma^2$ の正規分布を次の記号で表します。

　 ${\mathscr N}(x|\mu,\sigma^2) := \frac{1}{\sqrt{2\pi \sigma^2}}\exp\left\{-\frac{1}{2\sigma^2}(x-\mu)^2\right\}$

最尤推定

今回得られたトレーニングセットを与える情報源は、「何らかの関数 $y(x)$ に従って値が決まっているが、それぞれのデータには一定の正規分布の誤差が入り込む」という性質を持っていることが予めわかっているものとします。（ビジネス上のデータを解析する際は、そのデータを収集した背景から、このような付加情報が得られることもあるでしょう。）

そこで、この誤差の分散を $\beta^{-1}$ と仮定すると、パラメータ $\mathbf w$ が決定された暁には、説明変数 $x$ に対して、目的変数の値が $t$ になる確率は次式で与えられます。（正確には「確率密度」ですが、ここでは簡単のために「確率」と呼んでおきます。）

　 $p(t | x,{\mathbf w},\beta) = {\mathscr N}\left(t\,\middle|\,y(x,{\mathbf w}),\beta^{-1}\right)$

ということは、トレーニングセットに含まれる $N$ 個のデータが得られる確率は、次式になります。

　 $p({\mathbf t} | {\mathbf x},{\mathbf w},\beta) = \prod_{n=1}^N{\mathscr N}\left(t_n\,\middle|\,y(x_n,{\mathbf w}),\beta^{-1}\right)\\\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=\left(\frac{\beta}{2\pi}\right)^{\frac{N}{2}}\exp\left\{-\frac{\beta}{2}\sum_{n=1}^N\left(y(x_n,{\mathbf w})-t_n\right)^2\right\}\\\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,=\left(\frac{\beta}{2\pi}\right)^{\frac{N}{2}}\exp\left\{-\beta E_{\rm D}({\mathbf w})\right\}$

　 $E_{\rm D}({\mathbf w}) := \frac{1}{2}\sum_{n=1}^N\left(y(x_n,{\mathbf w})-t_n\right)^2$

　 $\mathbf{t} := \begin{pmatrix}t_1 \\ \vdots \\ t_N\end{pmatrix}\,\,\,\,\mathbf{x} := \begin{pmatrix}x_1 \\ \vdots \\ x_N\end{pmatrix}$

$N$ 個の観測点 $\{x_n\}_{n=1}^{N}$ を固定して何度も観測を繰り返すと、入り込む誤差の影響で、毎回異なる $\{t_n\}_{n=1}^{N}$ が得られるわけですが、実際に観測で得られたデータは、どのぐらいの確率で発生するものかを計算していると考えてください。（サイコロを2個振ったら、毎回異なる結果が出ますが、仮に「1」「1」のゾロ目が出たとして、これがどのぐらい確率の高い事象だったのかを後から振り返って考えているようなものです。）

また、興味深い点として、「PRML 第1章の多項式フィッティングの例を再現」で登場した誤差関数 $E_{\rm D}({\mathbf w})$ と同じものが上記の式に含まれています。

そして、ここで重大な想定（仮説）を置きます。

　「実際に観測されたデータは、発生する確率が一番高いものに違いない」

つまり、前述の $p({\mathbf t} | {\mathbf x},{\mathbf w},\beta)$ の値が最も大きくなるように、パラメータ $({\mathbf w}, \beta)$ を決定することにします。

　「そんな仮説、ほんまに正しいんかいな？」

と思うかもしれませんが、これが本当に正しいかどうかは、ここでは問題ではありません。前述のように「どのような方法でパラメータを決定するか」は、採用するモデルの一部であり、ここでは、「そういう仮説の下にパラメータを決定するモデル」を採用したというだけのことです。このモデルが実際に有用かどうか（未知のデータを正しく推定するかどうか）は、得られた結果を実際に利用した結果から判断することにしておきます。

パラメータ $({\mathbf w}, \beta)$ を決定する際は、先ほどの確率 $p({\mathbf t} | {\mathbf x},{\mathbf w},\beta)$ を $({\mathbf w}, \beta)$ の関数と見なすことになり、これを尤度関数と呼びます。そして、計算を簡単にするため、尤度関数そのものの代わりに、次の対数尤度関数を最大にします。すなわち、対数尤度関数の $({\mathbf w}, \beta)$ による偏微係数が 0 になる条件を求めます。

　 $\ln p({\mathbf t} | {\mathbf x},{\mathbf w},\beta)=-\beta E_{\rm D}({\mathbf w})+\frac{N}{2}\ln\beta-\frac{N}{2}\ln(2\pi)$

まず、 ${\mathbf w}$ に対する依存性を考えると、ちょうど $E_{\rm D}({\mathbf w})$ を通してのみ依存することが分かります。つまり、

　 $\frac{\partial (\ln p)}{\partial {\mathbf w}} = 0\,\,\,\,\Leftrightarrow\,\,\,\,\frac{\partial E_{\rm D}}{\partial {\mathbf w}} = 0$

より、 ${\mathbf w}$ は、 $E_{\rm D}$ を最小化するという条件（つまり、先に計算した、二乗平均平方根誤差を最小化するという条件）と同じ結果になります。

一方、

　 $\frac{\partial (\ln p)}{\partial \beta} = 0$

からは、下記の関係が得られます。これに先の条件で決定した ${\mathbf w}$ を代入すると $\beta^{-1}$ （データに含まれる未知の誤差の分散）が決まります。

　 $\beta^{-1}=\frac{1}{N}\sum_{n=1}^N\left\{y(x_n,{\mathbf w})-t_n\right\}^2$

これは、「トレーニングセットの分散」を「未知の誤差の分散」として採用することを表します。

ちなみに、ここでは $E_{\rm D}({\mathbf w})$ を最小化するという条件が結果として得られたわけですが、何か神秘的なものを感じたでしょうか？　実は、これは神秘的でも何でもなくて、未知の誤差の分布として正規分布を選んだからに過ぎません。言い換えると、尤度関数が誤差関数 $E_{\rm D}({\mathbf w})$ を通じて ${\mathbf w}$ に依存するような確率分布を探せば、結果的に正規分布が選ばれることになります。

それでは、この問題を実際に数値計算で解いてみましょう。多項式の推定は、最小二乗法と同じ結果になりますが、今回は、データに含まれる「ノイズ」の分散 $\sigma^2$ も同時に推定されます。つまり、推定結果は、多項式から得られる値に標準偏差 $\sigma$ の幅をもたせた形になります。次は、この幅をグラフに表示するコードになります。

# -*- coding: utf-8 -*-
#
# 最尤推定による回帰分析
#
# 2015/05/19 ver1.0
#

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame

from numpy.random import normal

#------------#
# Parameters #
#------------#
N=10            # サンプルを取得する位置 x の個数

# データセット {x_n,y_n} (n=1...N) を用意
def create_dataset(num):
    dataset = DataFrame(columns=['x','y'])
    for i in range(num):
        x = float(i)/float(num-1)
        y = np.sin(2*np.pi*x) + normal(scale=0.3)
        dataset = dataset.append(Series([x,y], index=['x','y']),
                                 ignore_index=True)
    return dataset

# 最大対数尤度（Maximum log likelihood）を計算
def log_likelihood(dataset, f):
    dev = 0.0
    n = float(len(dataset))
    for index, line in dataset.iterrows():
        x, y = line.x, line.y
        dev += (y - f(x))**2
    err = dev * 0.5
    beta = n / dev
    lp = -beta*err + 0.5*n*np.log(0.5*beta/np.pi)
    return lp

# 最尤推定で解を求める（解法は最小二乗法と同じ）
def resolve(dataset, m):
    t = dataset.y
    phi = DataFrame()
    for i in range(0,m+1):
        p = dataset.x**i
        p.name="x**%d" % i
        phi = pd.concat([phi,p], axis=1)
    tmp = np.linalg.inv(np.dot(phi.T, phi))
    ws = np.dot(np.dot(tmp, phi.T), t)

    def f(x):
        y = 0.0
        for i, w in enumerate(ws):
            y += w * (x ** i)
        return y

    sigma2 = 0.0
    for index, line in dataset.iterrows():
        sigma2 += (f(line.x)-line.y)**2
    sigma2 /= len(dataset)

    return (f, ws, np.sqrt(sigma2))

# Main
if __name__ == '__main__':
    train_set = create_dataset(N)
    test_set = create_dataset(N)
    df_ws = DataFrame()

    # 多項式近似の曲線を求めて表示
    fig = plt.figure()
    for c, m in enumerate([0,1,3,9]): # 多項式の次数
        f, ws, sigma = resolve(train_set, m)
        df_ws = df_ws.append(Series(ws,name="M=%d" % m))

        subplot = fig.add_subplot(2,2,c+1)
        subplot.set_xlim(-0.05,1.05)
        subplot.set_ylim(-1.5,1.5)
        subplot.set_title("M=%d" % m)

        # トレーニングセットを表示
        subplot.scatter(train_set.x, train_set.y, marker='o', color='blue')

        # 真の曲線を表示
        linex = np.arange(0,1.01,0.01)
        liney = np.sin(2*np.pi*linex)
        subplot.plot(linex, liney, color='green')

        # 多項式近似の曲線を表示
        linex = np.arange(0,1.01,0.01)
        liney = f(linex)
        label = "Sigma=%.2f" % sigma
        subplot.plot(linex, liney, color='red', label=label)
        subplot.plot(linex, liney+sigma, color='red', linestyle='--')
        subplot.plot(linex, liney-sigma, color='red', linestyle='--')
        subplot.legend(loc=1)

    fig.show()

    # 多項式近似に対する最大対数尤度を計算
    df = DataFrame()
    train_mlh = []
    test_mlh = []
    for m in range(0,9): # 多項式の次数
        f, ws, sigma = resolve(train_set, m)
        train_mlh.append(log_likelihood(train_set, f))
        test_mlh.append(log_likelihood(test_set, f))
    df = pd.concat([df,
                    DataFrame(train_mlh, columns=['Training']),
                    DataFrame(test_mlh, columns=['Test'])],
                    axis=1)
    df.plot(title='Maximum log likelihood for N=%d' % N)
    plt.show()

実行結果は次のようになります。

多項式の次数をあげるに従って、推定される標準偏差が小さくなっていきます。 $M=9$ の場合は $\sigma=0$ となり、与えられたデータがそのまま予測値になります。いわゆるオーバーフィッティングの状態です。

最尤推定のモデル評価

最小二乗法では、トレーニングセットとテストセットに対する二乗平均平方根誤差 $E_{\rm RMS}$ の変化からオーバーフィッティングを判定しました。つまり、トレーニングセットに対する $E_{\rm RMS}$ は、次数を上げればいくらでも減少しますが、テストセットに対する $E_{\rm RMS}$ は、ある程度以上には減らなくなりました。

これと同様の判定を最尤推定で行う場合は、「最大尤度」の変化を見ます。最尤推定は、(1)で計算した、下記の「トレーニングセットに含まれる N 個のデータが得られる確率」を最大化するようにパラメータを決定しました。

　 $p({\mathbf t} | {\mathbf x},{\mathbf w},\beta) = \left(\frac{\beta}{2\pi}\right)^{\frac{N}{2}}\exp\left\{-\beta E_{\rm D}({\mathbf w})\right\}$

　 $E_{\rm D}({\mathbf w}) := \frac{1}{2}\sum_{n=1}^N\left(y(x_n,{\mathbf w})-t_n\right)^2$

具体的に計算されるこの最大値の値が大きいほど、トレーニングセットが得られる確率、つまり、トレーニングセットに対する「当てはまり具合」は高いことになります。実際に次のコードで、トレーニングセットとテストセットに対する最大尤度の変化を見てみます。先ほどのコードは、最大尤度の変化を次のように表示します。

$M=3$ を越えるとテストセットに対する最大尤度は増加しなくなるので、ここからオーバーフィッティングが発生していると判断できます。

単純化した例で数値計算

最尤推定をよりよく理解するために、もう少し単純化した例を考えます。先の例では、観測点 $x$ によって正規分布の中心 $y(x,{\mathbf w})$ が異なっていましたが、ここでは、単一の観測点から複数のデータを取得して、それらを元にその点における正規分布の中心（平均）と分散を推定することにします。

まず、未知の正規分布を ${\mathscr N}\left(t\,\middle|\,\mu,\beta^{-1}\right)$ とすると、トレーニングセット $\{t_n\}_{n=1}^N$ が得られる確率（尤度関数）は次式になります。

　 $p({\mathbf t} | \mu,\beta) = \prod_{n=1}^N{\mathscr N}\left(t_n\,\middle|\,\mu,\beta^{-1}\right)$

　 $\ln p({\mathbf t} | \mu,\beta) = -\frac{\beta}{2}\sum_{n=1}^N(x_n-\mu)^2+\frac{N}{2}\ln\beta-\frac{N}{2}\ln(2\pi)$

対数尤度関数の $\mu,\beta$ による偏微分係数が 0 になるという条件からは、次式が得られます。

　 $\mu_{\rm ML} = \frac{1}{N}\sum_{n=1}^N t_n$

　 $\beta_{\rm ML}^{-1} = \frac{1}{N}\sum_{n=1}^N\left(t_n-\mu_{\rm ML}\right)^2$

つまり、トレーニングセットの平均／分散を未知の正規分布の平均／分散の推定量として採用することになります。 ${\rm ML}$ という添字は最尤推定（Maximum Likelihood）による推定量という事を表します。

それでは、これらの手法で推定される正規分布は、実際の正規分布とどのぐらい一致するのでしょうか？　容易に想像されるように、トレーニングセットのサンプルが少ない場合はあまり一致せずに、サンプルが増えるほどよく一致すると考えられます。サンプル数を変化させながら、実際に数値計算で試してみましょう。

# -*- coding: utf-8 -*-
#
# 最尤推定による正規分布の推定
#
# 2015/04/23 ver1.0
#

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame

from numpy.random import normal
from scipy.stats import norm

if __name__ == '__main__':
    fig = plt.figure()
    for c, datapoints in enumerate([2,4,10,100]): # サンプル数
        ds = normal(loc=0, scale=1, size=datapoints)
        mu = np.mean(ds)                # 平均の推定値
        sigma = np.sqrt(np.var(ds))     # 標準偏差の推定値

        subplot = fig.add_subplot(2,2,c+1)
        subplot.set_title("Datapoints=%d" % datapoints)
        # 真の曲線を表示
        linex = np.arange(-10,10.1,0.1)
        orig = norm(loc=0, scale=1)
        subplot.plot(linex, orig.pdf(linex), color='green')
        # 推定した曲線を表示
        est = norm(loc=mu, scale=np.sqrt(sigma))
        label = "Sigma=%.2f" % sigma
        subplot.plot(linex, est.pdf(linex), color='red', label=label)
        subplot.legend(loc=1)
        # サンプルの表示
        subplot.scatter(ds, orig.pdf(ds), marker='o', color='blue')
        subplot.set_xlim(-4,4)
        subplot.set_ylim(0)
    fig.show()

これを実行すると、次の結果が得られます。緑のカーブが元の正規分布で、平均 0、分散 1 になっています。青い点が得られたサンプルで、赤いカーブが前述の手法で推定された正規分布です。

予想通り、サンプルが少ない場合はうまくフィットしていませんが、特にカーブの幅が極端に狭くなっている（分散が小さく推定されている）点に注意してください。一般に、有限個のサンプルでは出現確率が低い端の方のデータが十分に得られないため、裾の広がりが小さく見積もられてしまうことになります。仮に $N$ 個のサンプルを取得して推定するということを何度も繰り返して、「推定値の平均値」を計算した場合、次の関係が成り立ちます。

　 ${\mathbb E}\left[\mu_{\rm ML}\right] = \mu$