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

多次元正規分布のパラメータを勾配法で求めることができないか,やってみました. 今回扱うモデルは,共分散が無く,分散が等しいモデルです. 共分散が無く,分散が等しい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

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

広告

フィボナッチ数列

ここによると,動的計画法の基礎は以下の三つに分かれるらしい. ビデオの中にも出て来たアルゴリズムを試してみた. フィボナッチ数列のプログラムを書くのは何年ぶりだろう.
例のごとく,ここにプログラムは置いている.source

# lazy
def fib(n):
    if n <= 2: f = 1
    else: f = fib(n-1) + fib(n-2)
    return f

# memorization
memo = {}
def fib1(n):
    if n in memo: return memo[n]
    if n <= 2: f = 1
    else: f = fib(n-1) + fib(n-2)
    memo[n] = f
    return f

# buttom up    
def fib2(n):
    fib = {}
    for k in range(1, n+1):
        if k <= 2: f = 1
        else: f = fib[k-1] + fib[k-2]
        fib[k] = f
    return fib[n]

if __name__ == '__main__':
    from time import time
    N = 35
    start = time()
    print fib(N)
    print time() - start

    start = time()
    print fib1(N)
    print time() - start

    start = time()
    print fib2(N)
    print time() - start

簡単な数値計算ライブラリを公開しました

簡単な数値計算ライブラリを公開しました.
https://github.com/nkt1546789/softcream/blob/master/calculus.py

1変数関数に関しては,その偏導関数を求めることができます. n変数関数に関しては,その勾配やヘッシアンを求めることができます. これらは全て関数で返ってくるので,例えば,

H = Hessian(f)
H(x)

などとしてxにおけるヘッシアンの値を求めることができます.Hはあくまでも関数です. 全てこのような形をとってあります. また,微分を簡単に近似計算しているので,それほど正確ではありませんが,実用上は問題ないかと思われます. この辺りはまた改良していきたいと思います.

次に,勾配法やニュートン法も実装してあります. 勾配法は最急降下法のみで,極小値を求めることができます. ニュートン法は極値を求めることができます.極小か極大かはわからないので注意.

勾配法,ニュートン法において明示的に関数の勾配や,ヘッシアンを与える必要がないのがこのライブラリの特徴です.