直線当てはめ-出力誤差

出力誤差モデル

#-*- coding:utf-8 -*-
from scipy import stats
import numpy as np
from pylab import *
u=0 #期待値を0と置く

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]

if __name__=="__main__":
	print "分散の値を入力してください。"
	v=input() #分散を１と置く
	xlist=np.linspace(-1,1,100)
	print "パラメータaの値を入力してください。"
	a=input()
	print "パラメータbの値を入力してください。"
	b=input()

	wlist=a*xlist+b
	plot(xlist,wlist,)

	ylist=[]
	for i in xlist:
		ylist.append(a*i+b+stats.norm.rvs(u,v))

	a1,b1=g(xlist,ylist)
	ylist=matrix(ylist)
	xlist=matrix(xlist)
	zlist=a1[0]*xlist+b1[0]
	plot(xlist,ylist,'r-*')
	plot(xlist,zlist,'g*')
	show()