ニュートン法

newton.py

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

import re
import urllib2
import numpy


def newton(X, y, initial_theta, num_iters=1500):
    m = X.shape[0]
    theta = numpy.copy(initial_theta)
    xp = numpy.copy(X)
    xp = numpy.insert(xp, 0, 1, axis=1)
    inverse_hesse = numpy.linalg.inv(xp.T.dot(xp))
    for i in xrange(num_iters):
        grad = (1.0/m) * numpy.sum((xp.dot(theta) - y)*xp, axis=0,
                                   keepdims=True)
        theta = theta - inverse_hesse.dot(grad.T)
    return theta


if __name__ == '__main__':
    datasrc = 'https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data'
    data = urllib2.urlopen(datasrc).read().replace('\\n', '').splitlines()
    data = numpy.array([[float(v) for v in re.split(" +", row.strip())]
                       for row in data])
    X = data[:, :-1]
    Y = data[:, -1:]
    initial_theta = numpy.zeros((data.shape[1], 1))
    num_iters = 1500

    # Normalize
    X = (X - numpy.mean(X, axis=0)) / numpy.std(X, axis=0)

    print newton(X, Y, initial_theta, num_iters)