curvefitting.py

#!/usr/bin/env python
# -*- coding: utf-8 -*-

import math
import csv
import numpy as np
import matplotlib.pyplot as plt

# M次多項式近似
M = 0
def y(x, w):
    ret = 0
    for i in range(len(w)):
        ret += w[i] * (x ** i)
    return ret

def sum_matrix(mat):
    # 行列の各要素の和を返却
    
    ret = 0
    # 行成分を取り出す
    for line in mat:
        # 行から要素を取り出す
        for element in line:
           ret += element
    return ret

def make_mat_A(xn, tn):
    A = [[0 for i in range(M + 1)] for row in range(M + 1)]
    for i in range(M + 1):
        for j in range(M + 1):
            sum_xn = 0
            N = len(xn)
            for n in range(N):
                sum_xn += xn[n] ** (i+j)
            A[i][j] = sum_xn
            
    return A

def make_vec_b(xn, tn):
    N = len(xn)
    b = []
    for i in range(M + 1):
        temp = 0
        for n in range(N):
            temp += ((xn[n] ** i) * tn[n])
        b.append(temp)

    return b
    
def read_data(file_name):
    # CSVデータの読み込み
    data = csv.reader(file(file_name, 'r'), delimiter=' ')
    data = [map(float, x) for x in data]
    
    # CSVデータの整理
    xn = []
    tn = []
    for row0, row1 in data:
        xn.append(row0)
        tn.append(row1)

    return (xn, tn)

def main():
    xn, tn = read_data('curvefitting.txt')
    xlist = np.arange(0.0, 1.01, 0.01)


    # m = 0 .. 10次多項式曲線
    for m in range(10):
        global M
        M = m
        A = make_mat_A(xn, tn)
        # 行列Aの逆行列
        inv_A = np.matrix(A).I

        b = make_vec_b(xn, tn)
        b = np.array(b)

        # 内積 A^(-1)・b
        w = np.dot(inv_A, b).tolist()[0]
        
        model = [y(x, w) for x in xlist]
        plt.plot(xlist, model)


    # サンプルデータを赤色点で打点
    plt.plot(xn, tn, 'bo')
    plt.plot(xlist, np.sin(2.0 * math.pi * xlist), 'g-')
    plt.show()

if __name__ == '__main__':
    main()