最小2乗法 - 数値計算とかの備忘録（仮）

ラグランジュ補間、スプライン補間と続いて今度は最小2乗法をやってみた。データからそれっぽい直線や曲線を求めるやつ。次のサイトが参考になる。

・http://szksrv.isc.chubu.ac.jp/lms/lms1.html
・http://www.geocities.jp/supermisosan/saisyounizyouhou2.html

pylabについては次のサイトが参考になる。

・http://www.ike-dyn.ritsumei.ac.jp/~uchida/scipy-lecture-notes/intro/matplotlib/matplotlib.html

コードはこんな感じ。今回はpylabを使ったのでちょっとだけ長め。

import pylab

def lsm(x, y, m):    #最小2乗法
	matrix = []
	ans = []
	for i in range(m+1):
		p = []
		t = 0.0
		for j in range(i, i + m + 1):
			s = 0.0
			for k in x:
				s += k ** j
			p.append(s)
		for j in range(len(x)):
			t += y[j] * (x[j] ** i)
		p.append(t)    #各行の右端の列には各行の式の解が入っている
		matrix.append(p)

	hakidashi(matrix, m+1)

	for i in matrix:
		ans.append(i[-1])

	return ans    #ansには次数の低い順に係数が格納されている

def hakidashi(matrix, dim):    #掃き出し法
	for i in range(dim):
		num = matrix[i][i]
		for j in range(dim + 1):
			matrix[i][j] = matrix[i][j] / num
		for j in range(dim):
			if i == j:
				pass
			else:
				a = matrix[j][i]
				for k in range(i, dim + 1):
					matrix[j][k] = matrix[j][k]- a * matrix[i][k]

def show(numlist, m):    #求めた関数を表示
	a = m
	ch = 'y = '
	while a >= 0:
		if a > 0:
			ch += '%f*x^%d' % (numlist[a], a)
			if numlist[a - 1] > 0:
				ch += '+'
		else:
			ch += '%f' % numlist[a]
		a -= 1
	print ch
	return ch

def make_graph(datx, daty, start, end, numlist, m, formula):    #グラフ描画
	a = m
	y = 0
	x = pylab.arange(start, end , 0.01)    #0.01刻みでyを求めることにする
	while a >= 0:    #求めた係数を使って関数を作る
		y += numlist[a] * (x ** a)
		a -= 1
	pylab.plot(x, y, label = formula)
	pylab.plot(datx, daty, 'ro', label = 'data')
	#pylab.legend(loc = 'upper right')    #凡例をつけることができる
	pylab.show()

def main():
	x = [0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0]
	y = [2.4824,1.9975,1.6662,1.3775,1.0933,0.7304,0.4344,0.2981,-0.0017,-0.0026]
	m = 2    #近似する関数の次数を設定
	ans = lsm(x, y, m)
	ch = show(ans, m)
	make_graph(x, y, x[0], x[-1] , ans, m, ch)

if __name__ == '__main__':
	main()