# PRML2-4章の実装

More than 1 year has passed since last update.

# 前書き

• 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()
```

# 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()
```

# ベイズ線形回帰

ほとんど技術評論社さんの「機械学習 はじめよう：第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()
```

# エビデンス関数の最大化

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

```#!/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()
```

データへのフィッティングがゆるくなっています。

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

```#!/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 &lt; 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 &lt; 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()
```

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

です。ロジスティック回帰が一番綺麗に分離できているように見えます。

# 後書き

また、多くの方の記事を参考にさせていただきました。