Python实现最小二乘法曲线拟合

感谢 最小二乘法曲线拟合以及Matlab实现 

作者 乐乐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)