logistic_regression.py

#!/usr/bin/python
# -*- coding: utf-8 -*-
"""
ロジスティック回帰の実験
参考 
「人工知能に関する断層録」
http://aidiary.hatenablog.com/entry/20100430/1272590402
"""

import numpy as np
import matplotlib.pyplot as plt

# 重みの初期化
def init_weights(i, j):
    return np.random.uniform(-1.0, 1.0, (i, j))

# バイアス項の追加
def add_bias(d):
    return np.hstack((np.ones((d.shape[0], 1)), d))

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

def phi(x):
    return x

def function(x, w):
    return np.dot(x, w.T)

def train(x, t, w):
    P = phi(x)
    y = sigmoid(function(P, w))

    R = np.diag(np.array(y*(1-y)).ravel())
    H = np.dot(np.dot(P.T, R), P)
    dE = np.dot(P.T, y-t)
    
    w -= np.dot(np.linalg.inv(H), dE).T
    
    return y, w

# set datas
N = 100
c1, c2 = [], []
m1, m2, m3 = [-1, 2], [1, -1], [8, -6]
cov = [[1,0.8], [0.8, 1]]
c1.extend(np.random.multivariate_normal(m1, cov, N/2))
c2.extend(np.random.multivariate_normal(m2, cov, N/2-1))
c2.extend(np.random.multivariate_normal(m3, cov, 1))
datas = np.vstack((c1, c2))
# set labels
labels = np.vstack((np.zeros((N/2, 1)), np.ones((N/2, 1))))

# add bias
datas = add_bias(datas)

# weights initialize
ip_num, op_num = datas.shape[1], labels.shape[1]
weights = init_weights(op_num,ip_num)

# training
for epoch in range(20):
    print("---epoch {}---".format(epoch))
    predict, weights = train(datas, labels, weights)

# plot area   
xmax = 10
ymax = 10
# line
x = np.arange(-xmax, xmax, 0.01)
y = -(weights[0,0]+weights[0,1]*x)/weights[0,2]
plt.plot(x, y, "g-")
# datas and labels
for (d, l) in zip(datas[:,1:], labels):    
    plt.scatter(d[0], d[1], c="r" if l > 0 else "b")
plt.xlim(-xmax, xmax)
plt.ylim(-ymax, ymax)
plt.show()