Python でPRML7.2章 RVMを実装。

いつ使うかは分からんが、RVMを実装しておこう。
7.2章の数式をそのまま書いた。動作確認は1次元のガウス基底で行った。

def rvm(t, PHI, alphas0, beta0, maxloop=10):
    """
    t: observed values.
    X: designe matrix
    alphas0: initial weights for each parameter.
    beta0: an initial precision value.
    
    returns
    alpha: prior weights
    beta: prior precision
    m: prior means
    V: prior covariance matrix
    
    See PRML 7.2 for details.
    
    ** to many loop causes 0 devision for some reason.
    """
    beta = beta0
    alphas = alphas0
    
    for i in range(maxloop):
        A = np.diag(alphas)
        V = np.linalg.inv(A+beta*np.dot(PHI.T, PHI))
        m = beta*np.dot(V, np.dot(PHI.T, t))
        gammas = 1 - alphas * np.diag(V)
        
        #new alphas and beta
        alphas = gammas / m**2 
        beta =  (N - np.sum(gammas)) / np.linalg.norm(t - np.dot(PHI, m))
        
    return alphas, beta, m, V

10個のガウス基底でフィッティングをした。
今回書いたコードではループの回数が多すぎるとzero divisionでエラーとなった。でもループは10回くらいで十分フィットしているように見えた。