多項式曲線フィッティング

fitting.py

#!/usr/bin/env python

import numpy as np
import numpy.random
import numpy.linalg
import matplotlib.pyplot as pyplot

RANGE = 2 * np.pi

def noise_sin(x):
    return np.sin(x) + numpy.random.randn() / 10

def fold(x, w):
    return sum([w[i] * np.power(x, i) for i in range(len(w))])

def error(x1, x2):
    return np.sqrt(sum([np.power(e[0] - e[1], 2) for e in zip(x1, x2)]) / len(x1))

def solve(x, t, m):
    A = [[sum([np.power(x_n, i + j) for x_n in x]) for j in range(m)] for i in range(m)]
    T = [sum([np.power(x[n], i) * t[n] for n in range(len(x))]) for i in range(m)]
    return np.linalg.solve(A, T)

def random_sample_point(num):
    return [numpy.random.rand() * RANGE for _ in range(num)]

def main():
    pyplot.figure(figsize = (16, 12))
    line_range = np.linspace(0.0, RANGE, 1000)
    sin_line = np.sin(line_range)

    num = int(raw_input("sample n >> "))
    training_x = random_sample_point(num)
    training_y = [noise_sin(x) for x in training_x]
    test_x = random_sample_point(num)
    test_y = [noise_sin(x) for x in test_x]

    training_error = []
    test_error = []

    model = range(1, 16)
    for i in model:
        pyplot.subplot(4, 4, i)
        pyplot.title("M = %d" % (i - 1))
        pyplot.plot(line_range, sin_line)
        pyplot.plot(training_x, training_y, "*")
        w = solve(training_x, training_y, i)
        func = lambda x: fold(x, w)
        line = map(func, line_range)
        y = map(func, training_x)
        training_error.append(error(training_y, y))
        pyplot.plot(line_range, line)

        test_estimate_y = map(func, test_x)
        test_error.append(error(test_y, test_estimate_y))

    pyplot.subplot(4, 4, 16)
    pyplot.plot(model, training_error, label = "training")
    pyplot.plot(model, test_error, label= "test")
    pyplot.legend(loc = "best")
    pyplot.show()

if __name__ == "__main__":
    main()