本文介绍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
拟合曲线和数据点的可视化结果为:
和手算的结果对比可发现他们是一致的: