最小二乗法

# -*- coding: UTF-8 -*-
from math import *
import numpy as np

############################################
# 線形連立方程式を解く（逆行列を求める）
#      w : 方程式の左辺及び右辺
#      n : 方程式の数
#      m : 方程式の右辺の列の数
#      eps : 逆行列の存在を判定する規準
#      return : =0 : 正常
#               =1 : 逆行列が存在しない
#      coded by Y.Suganuma
############################################

def gauss(w, n, m, eps) :

	nm  = n + m
	ind = 0

	for i1 in range(0, n) :

		y1 = 0.0
		m1 = i1 + 1
		m2 = 0
					# ピボット要素の選択
		for i2 in range(i1, n) :
			y2 = abs(w[i2][i1])
			if y1 < y2 :
				y1 = y2
				m2 = i2
					# 逆行列が存在しない
		if y1 < eps :
			ind = 1
			break
					# 逆行列が存在する
		else :   # 行の入れ替え
			for i2 in range(i1, nm) :
				y1        = w[i1][i2]
				w[i1][i2] = w[m2][i2]
				w[m2][i2] = y1
						# 掃き出し操作
			y1 = 1.0 / w[i1][i1]

			for i2 in range(m1, nm) :
				w[i1][i2] *= y1

			for i2 in range(0, n) :
				if i2 != i1 :
					for i3 in range(m1, nm) :
						w[i2][i3] -= (w[i2][i1] * w[i1][i3])

	return ind

#########################################
# 最小二乗法
#      m : 多項式の次数
#      n : データの数
#      x,y : データ
#      return : 多項式の係数（高次から）
#               エラーの場合はNULLを返す
#      coded by Y.Suganuma
#########################################

def least(m, n, x, y) :

	m += 1
	z  = np.empty(m, np.float)
	w  = np.empty((m, m+1), np.float)
	A  = np.empty((n, m), np.float)

	for i1 in range(0, n) :
		A[i1][m-2] = x[i1]
		A[i1][m-1] = 1.0
		x1         = A[i1][m-2]
		x2         = x1
		for i2 in range(m-3, -1, -1) :
			x2        *= x1
			A[i1][i2]  = x2

	for i1 in range(0, m) :
		for i2 in range(0, m) :
			w[i1][i2] = 0.0
			for i3 in range(0, n) :
				w[i1][i2] += A[i3][i1] * A[i3][i2]

	for i1 in range(0, m) :
		w[i1][m] = 0.0
		for i2 in range(0, n) :
			w[i1][m] += A[i2][i1] * y[i2]

	sw = gauss(w, m, 1, 1.0e-10)

	if sw == 0 :
		for i1 in range(0, m) :
			z[i1] = w[i1][m]
	else :
		z = np.empty(0, np.float)

	return z

----------------------------------

# -*- coding: UTF-8 -*-
import numpy as np
import sys
from math import *
from function import least

############################
# 最小二乗法（多項式近似）
#      coded by Y.Suganuma
############################

line = sys.stdin.readline()
s    = line.split()
m    = int(s[0])   # 多項式の次数
n    = int(s[1])   # データの数

x    = np.empty(n, np.float)
y    = np.empty(n, np.float)

for i1 in range(0, n) :   # データ
	line  = sys.stdin.readline()
	s     = line.split()
	x[i1] = float(s[0])
	y[i1] = float(s[1])

z = least(m, n, x, y)

if z.size  > 0 :
	print("結果")
	for i1 in range(0, m+1) :
		print("   " + str(m-i1) + " 次の係数 " + str(z[i1]))
else :
	print("***error  逆行列を求めることができませんでした")

------------データ例１（コメント部分を除いて下さい）--------
1 10   // 多項式の次数とデータの数
0  2.450   // 以下，(x, y)
1  2.615
2  3.276
3  3.294
4  3.778
5  4.009
6  3.920
7  4.267
8  4.805
9  5.656

-------------------------データ例２-------------------------
2 6
10  28.2
15  47.0
20  44.4
26  32.8
32  20.8
40   0.8