# PRML 5.1〜5.3：2層 Neural Network をベタに実装した結果

#### やる事

(x, y) 平面上のデータをロジスティクス回帰で二項分類します。

PRMLのFigure 5.4に相当する処理になります。

これをPythonでベタに実装すると次のコードになります。

# -*- coding: utf-8 -*-

import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
from pandas import Series, DataFrame
from numpy.random import multivariate_normal, rand

N1 = 10
Mu1 = [40,10]
N2 = 10
Mu2 = [10,-10]
N3 = 15
Mu3 = [5,8]

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

def prepare_dataset():
cov1 = np.array([[30,0],[0,80]])
cov2 = np.array([[80,0],[0,10]])
cov3 = np.array([[60,0],[0,20]])

df1 = DataFrame(multivariate_normal(Mu1,cov1,N1),columns=['x','y'])
df1['type'] = 1
df2 = DataFrame(multivariate_normal(Mu2,cov2,N2),columns=['x','y'])
df2['type'] = 1
df3 = DataFrame(multivariate_normal(Mu3,cov3,N3),columns=['x','y'])
df3['type'] = 0

df = pd.concat([df1,df2,df3],ignore_index=True)
df = df.reindex(np.random.permutation(df.index)).reset_index()
return df[['x','y','type']]

def run_simulation(tset, M, subplot1, subplot2):
tset1 = tset[tset['type']==1]
tset2 = tset[tset['type']==0]
ymin, ymax = tset.y.min()-5, tset.y.max()+10
xmin, xmax = tset.x.min()-5, tset.x.max()+10

bias = 0.5 * (tset.x.mean() + tset.y.mean())
w = {}

def f(x,y):
a = [0.0]*M
z = [0.0]*M
for j in range(1,M):
a[j] = w[(1,(j,0))]*bias + w[(1,(j,1))]*x + w[(1,(j,2))]*y
z[j] = np.tanh(a[j])
a_k = w[(2,(1,0))]*1.0
for j in range(1,M):
a_k += w[(2,(1,j))]*z[j]
return a_k, a

def initval():
return rand()-0.5

w[(2,(1,0))]=initval()
for j in range(1,M):
w[(1,(j,0))] = initval()
w[(1,(j,1))] = initval()
w[(1,(j,2))] = initval()
w[(2,(1,j))] = initval()

eta = 0.05
err_pre = 10000
for iter_num in range(1000):
for index, point in tset.iterrows():
x, y, type = point.x, point.y, point.type
a_k, a = f(x,y)
d_k = sigmoid(a_k) - type
w[(2,(1,0))] -= eta * d_k * 1
d = [0] * M
for j in range(1,M):
d[j] = d_k*w[(2,(1,j))]*(1.0-np.tanh(a[j])*np.tanh(a[j]))
w[(1,(j,0))] -= eta * d[j] * bias
w[(1,(j,1))] -= eta * d[j] * x
w[(1,(j,2))] -= eta * d[j] * y
w[(2,(1,j))] -= eta * d_k * np.tanh(a[j])
err = 0
for index, point in tset.iterrows():
x, y, type = point.x, point.y, point.type
a_k, a = f(x,y)
p = sigmoid(a_k)
err += -(type*np.log(p)+(1-type)*np.log(1.0-p))
if np.abs(err_pre - err) < err_pre * 0.00001:
print "break"
break
err_pre = err

print err_pre
RES = 80
X = np.linspace(xmin,xmax,RES)
Y = np.linspace(ymin,ymax,RES)
field1 = DataFrame(np.zeros(shape=(len(X),len(Y))))
field2 = DataFrame(np.zeros(shape=(len(X),len(Y))))
for x, xval in enumerate(X):
for y, yval in enumerate(Y):
a_k, a = f(xval,yval)
field1.ix[y,x] = np.sign(sigmoid(a_k)-0.5)
field2.ix[y,x] = sigmoid(a_k)

subplot1.set_title('M=%d, Err=%1.2f' % (M-1, err_pre))
subplot1.set_xticks([])
subplot1.set_yticks([])
subplot1.set_ylim([ymin-1, ymax+1])
subplot1.set_xlim([xmin-1, xmax+1])
subplot1.scatter(tset1.x, tset1.y, marker='o')
subplot1.scatter(tset2.x, tset2.y, marker='x')
subplot1.imshow(field1.values, extent=(xmin,xmax,ymax,ymin),
alpha=0.4)

subplot2.set_title('M=%d, Err=%1.2f' % (M-1, err_pre))
subplot2.set_xticks([])
subplot2.set_yticks([])
subplot2.set_ylim([ymin-1, ymax+1])
subplot2.set_xlim([xmin-1, xmax+1])
subplot2.scatter(tset1.x, tset1.y, marker='o')
subplot2.scatter(tset2.x, tset2.y, marker='x')
subplot2.imshow(field2.values, extent=(xmin,xmax,ymax,ymin),
alpha=0.4, vmin=0, vmax=1, cmap=plt.cm.gray_r)

# Main
if __name__ == '__main__':
fig = plt.figure()
tset = prepare_dataset()
for i in range(4):