機械学習（3） パーセプトロン

ここからは、この本を参考に進めていきます。

オンライン機械学習 (機械学習プロフェッショナルシリーズ)

しばらく二値分類を扱います。線形を前提とします。すなわち、空間を超平面で2つの部分に分割し、ランダムに点を打ち、その座標とどちらの部分に属するかを得ます。そのようなデータを一つずつ処理して、超平面の近似式を求めます。

実際のアルゴリズムは、本を読むか、下のコードを読み解くかするとして、
ざっくり言うと、仮の超平面の方程式を作っておいて、分類が成功していれば超平面はそのまま、失敗していれば補正します。補正するのですが、方程式の係数の絶対値がどんどん大きくなるので、結果的に少し補正することになるのです。

さて、ある次元で考えて、超平面をランダムに生成し、x_i ∊ [-1, 1]の点をランダムに生成し、何点処理したら、どれくらい超平面の方程式が正しくなるかを調べます。具体的には、補正された超平面と正しい超平面との成す角度を得ます。これが小さいほど正確に超平面が得られたとみなします。

2次元、10次元、100次元で1000サンプル、10000サンプル、100000サンプルをそれぞれ1000回ずつ二乗和を出して平均しました。

次元	1000サンプル	10000サンプル	100000サンプル
2	0.0433	0.0147	0.0064
10	0.1718	0.0800	0.0369
100	0.4307	0.1888	0.0865

どうやら、サンプル数が100倍になると精度が5倍、次元が100倍になると精度が1/5倍になるっぽいです。
とにかくたくさんサンプルが要ります。

from itertools import combinations, islice, izip
from math import sqrt, log, cos, sin, pi
import random
import sys
import time

def average(v):
	n, s = reduce(lambda (x, y), e: (x + 1, y + e), v, (0, 0))
	return s / n

# inner product
def ip(xs, ys):
	return sum(x * y for x, y in izip(xs, ys))

def Box_Muller():
	r1 = random.random()
	r2 = random.random()
	return sqrt(-2.0 * log(r1)) * cos(2.0 * pi * r2)

def random_hyper_sphere(D):
	xs = [ Box_Muller() for _ in range(D) ]
	r = sqrt(sum(x * x for x in xs))
	return [ x / r for x in xs ]

def make_hyper_plane(D):
	return random_hyper_sphere(D) + [ 0.0 ]

def gen_points(D):
	while True:
		yield [ random.uniform(-1.0, 1.0) for _ in range(D) ]

def perceptron(D, N):
	W = make_hyper_plane(D)
	
	w = [ 0.0 ] * (D + 1)
	for x in islice(gen_points(D), N):
		y = 1 if ip(W, x) > 0 else -1
		y1 = 1 if ip(w[:-1], x) + w[-1] > 0 else -1
		if y * y1 <= 0:
			w = [ v + y * x1 for v, x1 in izip(w, x) ] + [ w[-1] + y ]
	
	mag = sqrt(sum(v * v for v in w[:-1]))
	w1 = [ v / mag for v in w ]
	return sum((v - V) ** 2 for v, V in izip(w1, W))

def mean_diff(D, N, M):
	return average(perceptron(D, N) for _ in xrange(M))

t0 = time.clock()
D = int(sys.argv[1])	# dimension
N = int(sys.argv[2])	# number of samples
M = int(sys.argv[3])	# number of trial
print mean_diff(D, N, M)
print time.clock() - t0