作者 乐乐lelele
建议先看一下该作者的讲解,比较详细
最小二乘法拟合最终通过求解N次方程的系数来实现的 Y=X*A A=INV(X)*Y
#coding=utf-8
#独立测试脚本
#最小二乘法曲线拟合 测试脚本
#作者:Kare
#参考资料:https://blog.csdn.net/StupidAutofan/article/details/78924601
import numpy as np
from matplotlib import pyplot as plt
from KCurveFit import KCurveFit
#rank = 6
rank = 6
#初始化值
# x=[2,4,5,6,6.8,7.5,9,12,13.3,15]
# y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5]
x = np.arange(-1,1,0.02)
y = ((x*x-1)**3+1)*(np.cos(x*2)+0.6*np.sin(x*1.3))
len_x = len(x)
X0 = np.zeros((rank + 1, len_x)) #构造X0
for i in np.arange(rank + 1):
for j in np.arange(len_x):
X0[i, j] = np.power(x[j], rank - i)
#需要用到的函数
# np.transpose()
# np.linalg.inv
# np.dot()
#求系数
coefficient = np.zeros(len_x)
X = np.transpose(X0) #X0本身已经是转置,参见公式类型
XT = np.transpose(X) #求转置
XTX = np.dot(XT, X)
rx = np.linalg.matrix_rank(XTX)
XInv = np.linalg.inv(XTX) #求逆矩阵 先构造方阵
#系数A = INV((XT*X0)) * XT * Y
Y = np.transpose(y)
v1 = XInv
v2 = np.dot(v1, XT)
v3 = np.dot(v2, Y)
A = v3
print('系数矩阵\n', A)
# x0 = np.linspace(0,20,100)
x0 = np.arange(-1,1,0.002)
y0 = np.zeros(np.shape(x0))
for i in np.arange(len(x0)):
for j in np.arange(rank + 1):
y0[i] = y0[i] + A[j] * x0[i]**(rank-j)
plt.figure(1)
plt.grid(True)
plt.plot(x,y,'r.')
plt.hold
plt.plot(x0,y0,'b-')
plt.show()
#coding=utf-8
#最小二乘法曲线拟合
#函数
#作者:Kare
#参考资料:https://blog.csdn.net/StupidAutofan/article/details/78924601
import numpy as np
from matplotlib import pyplot as plt
def KCurveFit(x, y, rank=6, isPlot = True):
'''
:param x: x数据
:param y: y数据
:param rank: 曲线阶数
:param isPlot: 是否进行画图
:return: 计算正常则返回数组形式的系数,错误则返回0
'''
len_x = len(x)
X0 = np.zeros((rank + 1, len_x)) # 构造X0,X0是N次方程组的矩阵表示
for i in np.arange(rank + 1):
for j in np.arange(len_x):
X0[i, j] = x[j] ** (rank-i)
# 需要用到的函数
# np.transpose() 求转置
# np.linalg.inv 求逆矩阵
# np.dot() 矩阵乘积
# np.linalg.matrix_rank 求矩阵的秩
# 求系数 其实就是一个N次方程组求解的过程
X = np.transpose(X0) # X0本身已经是转置,参见公式类型
XT = np.transpose(X) # 求转置
XTX = np.dot(XT, X)
XInv = np.linalg.inv(XTX) # 求逆矩阵 先构造方阵
# 系数A = INV((XT*X)) * XT * Y
Y = np.transpose(y)
A = XInv.dot(XT).dot(Y) #得到系数
#验证参数是否正确
if isPlot:
# x0 = np.linspace(0,20,100)
minx = np.min(x)
maxx = np.max(x)
x0 = np.arange(minx, maxx, 0.001)
y0 = np.zeros(np.shape(x0))
for i in np.arange(len(x0)):
for j in np.arange(rank + 1):
y0[i] = y0[i] + A[j] * x0[i] ** (rank - j)
plt.figure(1)
plt.grid(True)
plt.plot(x, y, 'r.')
plt.hold
plt.plot(x0, y0, 'b-')
plt.show()
return A
#coding=utf-8
#调用函数脚本
#最小二乘法曲线拟合 测试脚本
#作者:Kare
#参考资料:https://blog.csdn.net/StupidAutofan/article/details/78924601
import numpy as np
# from matplotlib import pyplot as plt
from KCurveFit import KCurveFit
#rank = 6
rank = 6
#初始化值
# x=[2,4,5,6,6.8,7.5,9,12,13.3,15]
# y=[-10,-6.9,-4.2,-2,0,2.1,3,5.2,6.4,4.5]
x = np.arange(-1,1,0.02)
y = ((x*x-1)**3+1)*(np.cos(x*2)+0.6*np.sin(x*1.3))
A = KCurveFit(x, y, rank)
print('计算结果,系数 \n', A)