正規分布のパラメータを勾配法で求めてみる

多次元正規分布のパラメータを勾配法で求めることができないか,やってみました. 今回扱うモデルは,共分散が無く,分散が等しいモデルです. 共分散が無く,分散が等しいd次元正規分布の確率密度関数は以下のようになります.

\displaystyle  p(x; \mu, \sigma^2) = \frac{1}{(2\pi\sigma^2)^{\frac{d}{2}}}  \exp \left(-\frac{(x-\mu)^T(x-\mu)}{2\sigma^2} \right)

この分布からn個のデータが発生したとしてパラメータ\mu, \sigma^2を最尤推定を用いて求めてみます. 対数尤度関数は以下のようになります.

\displaystyle  \log L(\theta) = \sum^n_{i=1} \log p(x_i; \mu, \sigma^2)  = -\frac{nd}{2}\log 2\pi\sigma^2 - \frac{1}{2\sigma^2} \sum^n_{i=1} (xi-\mu)^T(xi-\mu)

これは解析的に解くことができますが,今回は勾配法を用いて解けるか,やってみました.

で作成した,明示的に勾配を与えなくても良い最急降下法のアルゴリズムを用います.

Source code

import numpy as np
from utils import rand_vcm, gradient_descent, Newtons_method

if __name__ == '__main__':
    d = 2
    n = 1000
    mu = np.random.normal(size=d)
    sigma = np.random.rand()
    x = np.random.multivariate_normal(mu, np.eye(d)*sigma**2, n)

    print "真のパラメータ"
    print mu
    print sigma**2

    print "最尤推定で解析的に求めた解"
    mu_mle = x.mean(axis=0)
    print mu_mle
    print np.sum([(xx-mu_mle).T.dot(xx-mu_mle) for xx in x])/(n*d)
    print 
    def logL(t):
        mu = t[:d]
        sigma = t[d]
        y = np.sum([(xx-mu_mle).T.dot(xx-mu_mle) for xx in x])
        return -0.5*(n*d*np.log(2*np.pi*sigma**2)+(y/sigma**2))

    t0 = np.random.random_sample(d+1)
    f = lambda t: -logL(t)
    t = gradient_descent(f, d+1, dx=1./10**8, h0=0.03)
    print "勾配法による解(2次元目までは平均,3次元目は分散)"
    print t

Result

真のパラメータ
[ 0.1698449   0.27160279]
0.0282500918966
最尤推定で解析的に求めた解
[ 0.16839332  0.26662506]
0.0294700109775

勾配法による解(2次元目までは平均,3次元目は分散)
[  6.18238880e-01   2.63384354e-01  -9.47000342e+07]

Summary

何回か試してみましたが,うまいこといかないですね. ステップ幅や初期値の問題でしょうか? これを解決していきたいです.

コメントを残す

以下に詳細を記入するか、アイコンをクリックしてログインしてください。

WordPress.com ロゴ

WordPress.com アカウントを使ってコメントしています。 ログアウト / 変更 )

Twitter 画像

Twitter アカウントを使ってコメントしています。 ログアウト / 変更 )

Facebook の写真

Facebook アカウントを使ってコメントしています。 ログアウト / 変更 )

Google+ フォト

Google+ アカウントを使ってコメントしています。 ログアウト / 変更 )

%s と連携中