前回は普通のニューラルネットワークをtensorflowやchainerなどのライブラリを使わずに実装しました。ただ、その類の実装はネット上にすでにたくさんあり、それではあまり面白みがありませんので、前回のはあくまで準備として今回は混合密度ネットワークを実装しました。前回実装したニューラルネットワークのコードを利用するので、前回の記事も適宜参考にしてみて下さい。

混合密度ネットワーク

従来モデルの問題点

回帰問題を解くときに、二乗和誤差関数を使っていたのではうまく対処できない場合があります。下の図（PRML図5.19の再現）は、青点をトレーニングデータとして、二乗和誤差関数を用いてニューラルネットワークを学習させたときの結果を表しています。

$x=0.5$のときは$y$の値が0.2、0.5、0.8のどれでもよさそうです。しかし、その中で一つを選ばないといけないので0.5になり、他の２つにフィッティングしません。ニューラルネットワークを学習させたのにあまりフィッティングできていないのは、コスト関数が単峰性のガウス分布でモデルされているからです。入力$x$とターゲット$t$の関係を、ガウス分布で以下のように表し、

p(t|x) = \mathcal{N}(t|\mu(x),\sigma^2),

これに基づいてフィッティング、関数$\mu(x)$の推定、を行っていました。そして、推定した$\mu(x)$が上の図の赤色となっていて、ダメダメな結果となってしまいました。

混合ガウス分布

混合密度ネットワークではコスト関数に多峰性の混合ガウス分布でモデルされたものを用いることで、上の問題を解決します。K個のガウス分布に一つ一つ混合係数$\pi_k$を用意して、入力$x$とターゲット$t$の関係を

p(t|x) = \sum_{k=1}^K\pi_k(x)\mathcal{N}(t|\mu_k(x),\sigma^2_k(x))

とモデルします。ただし、混合係数の和は1とします($\sum_k\pi_k=1$)。上のトレーニングデータの例だと、$x=0.5$では$y=0.2,0.5,0.8$のあたりに３つのデータ点の集合があります。よって、$K=3$とするのが良さそうです。3個のガウス分布がそれぞれ別のデータ点の集合にフィッティングされることで、多峰性のあるデータでも対処できます。
よって、トレーニングデータの組$\{x_n,t_n\}_{n=1}^N$が与えられているとき、最適化するコスト関数は

\begin{align}
E(w) &= \sum_{n=1}^N E_n\\
&= -\sum_{n=1}^N \ln p(t_n|x_n)\\
&= -\sum_{n=1}^N \ln \left\{\sum_{k=1}^K\pi_k(x_n,w)\mathcal{N}\left(t_n|\mu_k(x_n,w),\sigma^2_k(x_n,w)\right)\right\}
\end{align}

となります。そして、上の式を最小化するようなパラメータ$w$を推定します。関数$\pi_k,\mu_k,\sigma_k$にはニューラルネットワークを使います。パラメータ$w$のネットワークに入力$x$を入れると、混合係数、平均、分散が出力され、$t$の確率分布が算出されます。

ネットワーク構造

今回は上の図のようなデータに対して回帰を行うので、入力${\bf x}$は１次元となります。中間${\bf z}$のノードの数はPRMLにならって5つ、出力${\bf y}$のノード数は９つとします。第一層のパラメータを$W^{(1)},{\bf b}^{(1)}$、活性化関数を$f(\cdot)$、第二層のパラメータを$W^{(2)},{\bf b}^{(2)}$として、順伝播は次のようになります。

\begin{align}
{\bf z} &= f(W^{(1)}{\bf x} + {\bf b}^{(1)})\\
{\bf y} &= W^{(2)}{\bf z} + {\bf b}^{(2)}
\end{align}

この出力9つは混合ガウス分布のパラメータ$\pi_k,\mu_k,\sigma_k$の活性です。
まず、出力${\bf y}$は上から３つずつ$\pi,\mu,\sigma$に対応している活性とします。

{\bf y} =
\begin{bmatrix}
a_{\pi_1}\\
a_{\pi_2}\\
a_{\pi_3}\\
a_{\mu_1}\\
a_{\mu_2}\\
a_{\mu_3}\\
a_{\sigma_1}\\
a_{\sigma_2}\\
a_{\sigma_3}
\end{bmatrix}

混合係数$\pi_k$は$\sum_k\pi_k=1$という制約条件があるので、ソフトマックス関数により混合係数を出力します。

\pi_k = {\exp(a_{\pi_k})\over\sum_{l=1}^K\exp(a_{\pi_l})}

平均は非線形変換を用いずに$\mu_k = a_{\mu_k}$とします。
標準偏差$\sigma$は０以上なので指数関数を使って$\sigma_k = \exp(a_{\sigma_k})$。
これで、混合ガウス分布を計算するのに必要な関数$\pi_k,\mu_k,\sigma_k$が求まりました。

勾配

前回のニューラルネットワークの実装で書いたようにニューラルネットワークを学習させるのにはコスト関数をニューラルネットワークの出力で微分した勾配が必要です。それがあれば、あとは誤差逆伝播を使ってコスト関数のパラメータ$w$についての微分も計算できます。
そのコスト関数のニューラルネットワークの出力での勾配は、

\gamma_{nk}(t_n|x_n) = {\pi_k\mathcal{N}(t_n|\mu_k(x_n,w),\sigma_k^2(x_n,w))\over\sum_l\pi_l\mathcal{N}(t_n|\mu_l(x_n,w),\sigma_l^2(x_n,w))}

を用いて、

\begin{align}
{\partial E_n\over\partial a_{\pi_k}} &= \pi_k - \gamma_{nk}\\
{\partial E_n\over\partial a_{\mu_k}} &= \gamma_{nk}{\mu_k - t_n\over\sigma^2_k}\\
{\partial E_n\over\partial a_{\sigma_k}} &= \gamma_{nk}\left(1 - {||t_n - \mu_k||^2\over \sigma^2_k}\right)
\end{align}

となります。式の導出はPRMLでは演習問題になっているのでその回答、もしくは混合密度ネットワークの論文を参照してください。

実装

コスト関数

上で書いた数式をコードにすると、次のようになります。

# コスト関数のクラス
class GaussianMixture(object):
    """Negative Log Likelihood of Gaussian Mixture model"""
    def __init__(self, n_components):
        # ガウス分布の個数、今までの例だと3
        self.n_components = n_components

    # コスト関数の値を計算するメソッド
    def __call__(self, X, targets):

        # ネットワークの出力Xを活性化関数で変換して、標準偏差、混合係数、平均を計算
        sigma, weight, mu = self.activate(X)

        # ガウス関数N(t|mu,sigma^2)の値を計算
        gauss = self.gauss(mu, sigma, targets)

        # 混合ガウス分布の負の対数尤度E(w):PRML(式5.153)
        return -np.sum(np.log(np.sum(weight * gauss, axis=1)))

    # 活性化関数で変換
    def activate(self, X):
        assert np.size(X, 1) == 3 * self.n_components

        # Xを標準偏差、混合係数、平均に対応するところに分割
        X_sigma, X_weight, X_mu = np.split(X, [self.n_components, 2 * self.n_components], axis=1)

        # 標準偏差を活性化関数で変換
        sigma = np.exp(X_sigma)

        # 混合係数を活性化関数で変換、桁があふれないように最大値で引く
        weight = np.exp(X_weight - np.max(X_weight, 1, keepdims=True))
        weight /= np.sum(weight, axis=1, keepdims=True)

        return sigma, weight, X_mu

    # ガウス関数N(target|mu,sigma^2)を計算
    def gauss(self, mu, sigma, targets):
        return np.exp(-0.5 * (mu - targets) ** 2 / np.square(sigma)) / np.sqrt(2 * np.pi * np.square(sigma))

    # コスト関数を活性で微分
    def delta(self, X, targets):
        sigma, weight, mu = self.activate(X)
        var = np.square(sigma)
        gamma = weight * self.gauss(mu, sigma, targets)
        gamma /= np.sum(gamma, axis=1, keepdims=True)

        # それぞれの微分を計算
        delta_mu = gamma * (mu - targets) / var
        delta_sigma = gamma * (1 - (mu - targets) ** 2 / var)
        delta_weight = weight - gamma

        # 連結させてから返す
        delta = np.hstack([delta_sigma, delta_weight, delta_mu])
        return delta

学習データの作成

# 関数funcの逆関数に従うような学習データを作成する
def create_toy_dataset(func, n=300):
    t = np.random.uniform(size=(n, 1))
    x = func(t) + np.random.uniform(-0.05, 0.05, size=(n, 1))
    return x, t

ミニバッチ処理

学習がうまくいかないことがあったので急遽使用した手法の一つ。
学習データ全てを使って勾配を計算するのではなく、学習データからいくつかを取ってきて勾配計算を行う。こうすることで、局所解から抜け出せるかもしれません。

# 学習データの組x,tからランダムにペアをn個選択する
def sample(x, t, n=None):
    assert len(x) == len(t)
    N = len(x)
    if n is None:
        n = N
    indices = np.random.choice(N, n, replace=False)
    return x[indices], t[indices]