最小二乗法フィッティング

サンプルデータ生成

from numpy.random import *

fout = open("out.txt",'w')

a = 20
b = -15

for x in range(100):
    y = a*x + b + normal(0,20)
    fout.write("{0},{1}\n".format(x,y))

fout.close()

正規方程式を解く

import numpy

fin = open('out.txt','r')

a11 = 0
a12 = 0
a21 = 0
a22 = 0
b1 = 0
b2 = 0

for line in fin:
    line = line.rstrip()
    x = float(line.split(',')[0])
    y = float(line.split(',')[1])
    a11 = a11+x**2
    a12 = a12+x
    a21 = a21+x
    a22 = a22 + 1
    b1 = b1 + x*y
    b2 = b2 + y

A = numpy.matrix([[a11,a12],[a21,a22]])
b = numpy.matrix([b1,b2]).transpose() 
x = numpy.linalg.solve(A,b)

print x

fin.close()

逆行列を公式で展開すると、連立方程式の求解が楽になる

import numpy

fin = open('out.txt','r')

a11 = 0
a12 = 0
a21 = 0
a22 = 0
b1 = 0
b2 = 0

sxx = 0
sxy = 0
sy = 0
sx = 0
N = 0

for line in fin:
    line = line.rstrip()
    x = float(line.split(',')[0])
    y = float(line.split(',')[1])
    a11 = a11+x**2
    a12 = a12+x
    a21 = a21+x
    a22 = a22 + 1
    b1 = b1 + x*y
    b2 = b2 + y
    N = N +1
    sxx = sxx + x**2
    sy = sy + y
    sxy = sxy + x*y
    sx = sx + x

A = numpy.matrix([[a11,a12],[a21,a22]])
b = numpy.matrix([b1,b2]).transpose() 
x = numpy.linalg.solve(A,b)

print x

a = (-sxy*N+sx*sy)/(sxx*N-sx**2)
b = (sx*sxy-sy*sxx)/(sxx*N-sx**2)

print a,b

fin.close()