PRML4.2.1 pythonで実装してみた。

PRML４．２章を実装してみた。

とりあえず可視化は成功。
共分散行列が同じでない場合とか
事前確率が違う場合とかも試せてわかってイメージが湧いた。

meshgrid()を使って２次元パラメータの可視化を行ったけど、
作ったメッシュで２次元ガウス分布をもっと効率的にできんかね？（６５，６６行のメソッド）
少し時間が掛かってしまう。
まぁよしとしよう。

matplotlibで等高線に凡例をつける場合はググると次のようにすることがわかった。

    lines=[]
    lines.append(plt.contour(X1, X2, px_C1, colors="r"))
    lines.append(plt.contour(X1, X2, px_C2, colors="b"))
    actual_lines = [cs.collections[0] for cs in lines] 
    plt.legend(actual_lines, ("px_C1", "px_C2")) 

以下ソースを載っけます。

import numpy as np
from matplotlib import pyplot as plt



class Gauss2D:
    def __init__(self, m, sig):
        self.m = m
        self.d = 2
        self.sig = sig
        self.sig_i = np.linalg.inv(sig)
        
        det = np.linalg.det(sig)**0.5
        self.c = 1/(2*np.pi*det)
        
        
    
    def meshvalue(self, X1, X2):
        """
        get values by meshgrid X1, X2
        return the array of which size is the same with X1 and X2
        SLOW!!!
        """
        if X1.shape != X2.shape:
            raise TypeError("input array shape are not the same!")
        
        Z = np.empty(X1.shape)
        
        for i in range(X1.shape[0]):
            for j in range(X1.shape[1]):
                tmp = np.array([X1[i,j], X2[i,j]]) - self.m
                Z[i,j] = self.c * np.exp(-0.5 * np.dot(np.dot(tmp, self.sig_i), tmp))
        
        return Z
            

    

if __name__ == "__main__":
    fig = plt.figure(figsize=(6,16))
    
    
        
    x1 = np.arange(-1, 1.0, 0.01)
    x2 = np.arange(-1, 1.0, 0.01)    
    X1, X2 = np.meshgrid(x1,x2)
       
    pC1 = 0.5
    pC2 = 1 - pC1    
    
    """define px_C1 and px_C2 as gaussian"""
    m1 = np.array([-0.5, -0.5]) #means
    m2 = np.array([0.5, 0.5])
    sig1 = np.array([[0.1,0], [0,0.1]]) #covarance
    sig2 = np.array([[0.1,0], [0,0.1]])
    g1 = Gauss2D(m1, sig1)
    g2 = Gauss2D(m2, sig2)
    px_C1 = g1.meshvalue(X1,X2)
    px_C2 = g2.meshvalue(X1,X2)
    
    """plot px_C1 and px_C2"""
    fig.add_subplot(211)
    lines=[]
    lines.append(plt.contour(X1, X2, px_C1, colors="r"))
    lines.append(plt.contour(X1, X2, px_C2, colors="b"))
    actual_lines = [cs.collections[0] for cs in lines] 
    plt.legend(actual_lines, ("px_C1", "px_C2")) 

    """calc and plot pC1_x"""
    a = np.log(px_C1*pC1/(px_C2*pC2))
    pC1_x = 1/(1+np.exp(-a))            #sigmoid   
    
    fig.add_subplot(212)        
    plt.imshow(pC1_x[::-1], extent=[-1,1,-1,1], interpolation='nearest')
    plt.colorbar(shrink=0.8,aspect=10)
    plt.title("pC1_x")
    
    
    plt.show()