本文介绍python实现最小二乘拟合/平方逼近问题的方法
最优平方逼近问题的定义为:
可以使用正规方程组求解:
解方程即得
(
c
0
,
c
1
,
.
.
.
,
c
n
)
(c_0, c_1, ..., c_n)
(c0,c1,...,cn),即为拟合式的系数。
当拟合式为多项式时:
正规方程组可以转化为:
用python语言描述正规方程组的求解过程如下:
def mox(n):
sum = 0
for i in x:
sum += pow(i, n)
return sum
def moy(n):
sum = 0
for i, j in zip(x, y):
sum += pow(i, n) * j
return sum
def p(i):
p = 0
for idx, k in enumerate(c):
p += k * pow(i, idx)
return p
def compute_error(c):
error = 0
for i, j in zip(x, y):
error += pow(p(i)-j, 2)
return math.sqrt(error)
def compute(n):
global c
A = np.zeros((n, n))
for i in range(n):
for j in range(n):
A[i, j] = A[j, i] = mox(i+j)
b = np.zeros(n)
for i in range(n):
b[i] = moy(i)
print('A = ', A)
print('b = ', b)
c = np.matmul(np.linalg.inv(A), b)
print('c = ', c)
error = compute_error(c)
print('error = ', error)
使用一道例题来检验一下:
输入例题:
x = [0.25*i for i in range(5)]
y = [1.0000, 1.2840, 1.6487, 2.1170, 2.7183]
compute(n=3)
输出结果:
A = [[5. 2.5 1.875 ]
[2.5 1.875 1.5625 ]
[1.875 1.5625 1.3828125]]
b = [8.768 5.4514 4.4015375]
c = [1.00513714 0.86418286 0.84365714]
error = 0.016556949339433698
拟合曲线和数据点的可视化结果为:
和手算的结果对比可发现他们是一致的: