第四章 最小二乗法-「これなら分かる最適化数学」

#二次元データで正規方程式を解く
#-*- coding:utf-8 -*-
from numpy import *
from pylab import *
#例題4.2
#x=[4,15,30,100,200]
#y=[-17,-4,-7,50,70]

#例題4.3
#x=[2.8,2.9,3.0,3.1,3.2,3.2,3.2,3.3,3.4]
#y=[30,26,33,31,33,35,37,36,33]

#正規方程式
def g(A,B):
	x_11=np.dot(A,A)
	x_12=sum(A)
	x_21=sum(A)
	x_22=len(A)
	y_11=np.dot(A,B)
	y_21=sum(B)
	A=matrix([[x_11,x_12],[x_21,x_22]])
	b=matrix([[y_11],[y_21]])
	C=dot(A.I,b)
	return C[0],C[1]

def main():
	a,b=g(x,y)
	plot(x,y,'r*')
	xlist=matrix(linspace(0,200,1000))
	ylist=a*xlist+b
	plot(xlist,ylist,'b*')
	show()

if __name__=="__main__":
	main()

つづく。