サンプルデータ生成
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()