Python
PRML
機械学習

PRML2-4章の実装

More than 1 year has passed since last update.

前書き

春頃にPRMLの2-4章を読んで、いくつか実装を行いました。
別のブログに記事を書いていたものをこちらに移行します。
使用言語はPython(3系)で、実装したのは以下です。

  • 2章
    • カーネル密度推定
    • k-最近傍法(密度推定)
  • 3章
    • ベイズ線形回帰
    • エビデンス関数の最大化
  • 4章
    • 線形識別関数(最小二乗、フィッシャー、パーセプトロン)
    • ロジスティック回帰

カーネル密度推定

下記コードは動きますが、多次元の場合にクラスが動くかは確認していません。
nykergotoさんの「カーネル密度推定 - nykergoto’s blog」を参考にしました。

#!/usr/bin/env python3

import sys
import numpy as np
from matplotlib import pyplot as plt

class KernelDensityEstimation:
    '''kernel density estimation'''

    def __init__(self, sample, h=0.5):
        self.sample = sample
        self.h = h
        try:
            self.D = self.sample.shape[1]
        except IndexError as e:
            self.D = 1

    def estimate(self, x):
        p = np.zeros(x.size)
        for i in range(x.size):
            vec = self.sample - x[i]
            p[i] = np.mean((np.exp(-np.linalg.norm(vec, axis=1) / (2 * self.h ** 2))) / ((2 * np.pi * self.h ** 2) ** (self.D / 2)))
        return p

def main():
    sample = np.random.randn(1000,1) #np.random.rand(データ数, 次元)
    kernel = KernelDensityEstimation(sample)
    x = np.linspace(-5, 5, 100)
    plt.hist(sample,normed=True, facecolor='green',alpha=0.3)
    plt.plot(x,kernel.estimate(x),"b--")
    plt.show()

if __name__ == "__main__":
    main()

kernel_density_estimation-300x225.png

k-最近傍法(密度推定)

#!/usr/bin/env python3

import numpy as np
from scipy import special
from matplotlib import pyplot as plt

class Kneighbors:
    '''k-neighbors'''

    def __init__(self, sample, k=30):
        self.sample = sample
        self.k = k
        try:
            self.D = self.sample.shape[1]
        except IndexError as e:
            self.D = 1
        self.n = sample.size / self.D

    def estimate(self, x):
        p = np.zeros(x.size)
        for i in range(x.size):
            vec = self.sample - x[i]
            length = np.linalg.norm(vec, axis=1)
            sorted_length = np.sort(length)
            p[i] = self.k * special.gamma(self.D/2 + 1) / (self.n * (np.pi ** (self.D/2)) * (sorted_length[self.k - 1] ** (self.D)))
        return p


def main():
    sample = np.random.randn(1000, 1)
    kn = Kneighbors(sample)
    x = np.linspace(-5, 5, 100)
    plt.hist(sample,normed=True, facecolor='green',alpha=0.3)
    plt.plot(x, kn.estimate(x),"b--")
    plt.show() 


if __name__ == "__main__":
    main()

k_neighbors-300x225.png

ベイズ線形回帰

ほとんど技術評論社さんの「機械学習 はじめよう:第14回 ベイズ線形回帰を実装してみよう」の写しになってます。
一応クラスに変えた程度。

#!/usr/bin/env python3

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns

def phi(x): 
    s = 0.1
    return np.append(1, np.exp(-(x - np.arange(0, 1 + s, s)) ** 2 / (2 * s * s)))

class BayesianLinearRegression:

    def __init__(self, X, t, alpha=0.1, beta=9.0):
        self.alpha = alpha
        self.beta = beta
        self.train(X, t)

    def train(self, X, t):
        PHI = np.array([phi(x) for x in X])
        Sigma_N = np.linalg.inv(self.alpha * np.identity(PHI.shape[1]) + self.beta * np.dot(PHI.T, PHI))
        self.mu_N = self.beta * np.dot(Sigma_N, np.dot(PHI.T, t))

    def estimate(self, xlist):
        return [np.dot(self.mu_N, phi(x)) for x in xlist]

def main():
    X = np.array([0.02, 0.12, 0.19, 0.27, 0.42, 0.51, 0.64, 0.84, 0.88, 0.99])
    t = np.array([0.05, 0.87, 0.94, 0.92, 0.54, -0.11, -0.78, -0.79, -0.89, -0.04])
    blr = BayesianLinearRegression(X, t)
    xlist = np.arange(0, 1, 0.01) 
    plt.plot(xlist, blr.estimate(xlist), 'b') 
    plt.plot(X, t, 'o') 
    plt.show()

if __name__ == "__main__":
    main()

bayesian_linear_regression-300x206.png

エビデンス関数の最大化

PRML3章からエビデンス関数の最大化の実装
Gordian_knot さんの 「PRML第3章 エビデンス近似 Python実装」が参考になりました。
実装するだけならPRMLの式をそのままです。

#!/usr/bin/env python3

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns

def phi(x): 
    s = 0.1
    return np.append(1, np.exp(-(x - np.arange(0, 1 + s, s)) ** 2 / (2 * s * s)))

class BayesianLinearRegression:
    '''Bayesian Linear Regression'''
    def __init__(self, X, t, alpha=0.1, beta=9.0):
        self.alpha = alpha
        self.beta = beta
        self.PHI = np.array([phi(x) for x in X])
        self._train(t)

    def _train(self, t):
        self.S_N = np.linalg.inv(self.alpha * np.identity(self.PHI.shape[1]) + self.beta * np.dot(self.PHI.T, self.PHI))
        self.m_N = self.beta * np.dot(self.S_N, np.dot(self.PHI.T, t))

    def estimate(self, xlist):
        return [np.dot(self.m_N, phi(x)) for x in xlist]

class EvidenceApproximation(BayesianLinearRegression):
    '''Evidence Approximation'''
    def __init__(self, X, t, alpha=0.1, beta=9.0, iter_max=100):
        self.alpha = alpha
        self.beta = beta
        self.PHI = np.array([phi(x) for x in X])
        self.N = X.size
        self.iter_max = iter_max
        self._approximate(t)

    def _approximate(self, t):
        for i in range(self.iter_max):
            self._train(t)
            eigenvalues = np.linalg.eigvals(self.beta * self.PHI.T.dot(self.PHI))
            self.gamma = np.sum(eigenvalues / (self.alpha + eigenvalues))
            self.alpha = self.gamma / self.m_N.T.dot(self.m_N)
            self.beta = (self.N - self.gamma) / np.sum((t - self.PHI.dot(self.m_N)) ** 2)

    def get_params(self):
        return self.alpha, self.beta

def main():
    X = np.array([0.02, 0.12, 0.19, 0.27, 0.42, 0.51, 0.64, 0.84, 0.88, 0.99])
    t = np.array([0.05, 0.87, 0.94, 0.92, 0.54, -0.11, -0.78, -0.79, -0.89, -0.04])
    blr = BayesianLinearRegression(X, t)
    ea = EvidenceApproximation(X, t)
    a, b = ea.get_params()
    xlist = np.arange(0, 1, 0.01)
    plt.plot(xlist, blr.estimate(xlist), 'b', label='alpha=0.1, beta=9.0') 
    plt.plot(xlist, ea.estimate(xlist), 'g', label=f'alpha={a:.3f}, beta={b:.3f}') 
    plt.plot(X, t, 'o')
    plt.legend()
    plt.show()

if __name__ == "__main__":
    main()

evidence_approximation-300x206.png

青はベイズ前節のものと同じグラフで、緑が今回のものです。
データへのフィッティングがゆるくなっています。

線形識別関数とロジスティック回帰

最小二乗、フィッシャーの識別関数、パーセプトロンアルゴリズムの実装については以下を参考にしました。

#!/usr/bin/env python3

import numpy as np
from matplotlib import pyplot as plt
import seaborn as sns

def least_squares(X, T):
    N = X.shape[0]
    ones = np.ones((N, 1))
    X_ = np.hstack((ones, X))
    W_ = np.linalg.pinv(X_).dot(T)
    W_t = W_.T
    a = - (W_t[0, 1] - W_t[1, 1]) / (W_t[0, 2] - W_t[1, 2])
    b = - (W_t[0, 0] - W_t[1, 0]) / (W_t[0, 2] - W_t[1, 2])
    return lambda x: a * x + b

def fisher(x1, x2):
    m1 = np.mean(x1, axis=0)
    m2 = np.mean(x2, axis=0)
    m = np.mean(np.vstack((x1, x2)), axis=0)
    S_w = (x1 - m1).T.dot(x1 - m1) + (x2 - m2).T.dot(x1 - m2)
    w = np.linalg.inv(S_w).dot(m2 - m1)
    a = - (w[0] / w[1])
    b = - a * m[0] + m[1]
    return lambda x: a * x + b

def perceptron(X, T, bias=1.0, iter_max=100):
    N = X.shape[0]
    ones = np.ones((N, 1)) * bias
    X_ = np.hstack((ones, X))
    W = np.array([1.0, 0.0, 0.0])
    for i in range(iter_max):
        temp = W.dot(X_.T) * T
        failures = [temp < 0.0]
        W = W + T[failures].dot(X_[failures])
    a = - W[1] / W[2]
    b = - (W[0]/ W[2]) * bias
    return lambda x: a * x + b

def sigmoid(a):
    return 1 / (1 + np.exp(-a))

class LogisticRegression:
    '''Logistic Regression'''

    def __init__(self, X, T):
        self.N = T.size
        ones = np.ones((self.N, 1))
        self.X = np.hstack((ones, X))
        self.T = T
        self.M = 3
        self.w = np.zeros(self.M)

    def _train(self, iter_max):
        for i in range(iter_max):
            Y = sigmoid(self.w.dot(self.X.T))
            R = np.diag(Y * (1.0 - Y))
            H = self.X.T.dot(R).dot(self.X)
            w_new = self.w - np.linalg.inv(H).dot(self.X.T).dot(Y - self.T)
            if np.linalg.norm(self.w):
                diff = np.linalg.norm(w_new - self.w) / np.linalg.norm(self.w)
                if diff < 0.1:
                    return self.w
            self.w = w_new
        return self.w

    def estimate(self, iter_max=100):
        W = self._train(iter_max)
        a = - W[1] / W[2]
        b = - W[0] / W[2]
        return lambda x: a * x + b


def funcs_subplot(x_num, y_num, index, xlist, func):
    ylist = [func(x) for x in xlist]
    plt.subplot(x_num, y_num, index)
    plt.scatter(x_red[:, 0], x_red[:, 1], color='r', marker='x')
    plt.scatter(x_blue[:, 0], x_blue[:, 1], color='b', marker='x')
    plt.plot(xlist, ylist, 'g-')
    plt.xlim(-4.0, 9.0)
    plt.ylim(-9.0, 4.0)
    return index

def main():
    xlist = np.linspace(-4, 8, 1000)
    index = 0

    # least_squares
    plt.figure(figsize=(6, 6))
    ls_f = least_squares(X, T)
    index = funcs_subplot(2, 2, index + 1, xlist, ls_f)

    # fisher
    fi_f = fisher(x_red, x_blue)
    index = funcs_subplot(2, 2, index + 1, xlist, fi_f)

    # perceptron
    pe_f = perceptron(X, T2)
    index = funcs_subplot(2, 2, index + 1, xlist, pe_f)

    # logistic regression
    lr = LogisticRegression(X, T3)
    lr_f = lr.estimate()
    index = funcs_subplot(2, 2, index + 1, xlist, lr_f)    

    plt.show()

if __name__ == "__main__":
    # generate data
    mu_red = [-1, 1]
    mu_blue = [1, -1]
    cov = [[1.2, 1.0], [1.0, 1.2]]
    mu_b_out= [7,-6]
    cov2 = [[1,0], [0,1]]
    x_red = np.random.multivariate_normal(mu_red, cov, 50)
    x_blue = np.vstack((np.random.multivariate_normal(mu_blue, cov, 45),
                            np.random.multivariate_normal(mu_b_out, cov2, 5)))
    X = np.vstack((x_red, x_blue))
    T_red = np.matlib.repmat(np.array([1, 0]), 50, 1)
    T_blue = np.matlib.repmat(np.array([0, 1]), 50, 1)
    T = np.concatenate((T_red, T_blue), axis=0)
    T2 = np.array([1] * 50 + [-1] * 50)
    T3 = np.array([1.0] * 50 + [0.0] * 50)

    main()

logistic_regression.png

  • 左上:最小二乗
  • 右上:フィッシャー
  • 左下:パーセプトロン
  • 右下:ロジスティック回帰

です。ロジスティック回帰が一番綺麗に分離できているように見えます。
最小二乗やフィッシャーでは、境界面から遠いデータの影響を受けて分類がうまくいっていません。

後書き

後から見返すと、微妙なところもあるコードですね。
また、多くの方の記事を参考にさせていただきました。