python手写多项式拟合、曲线拟合

深藏阁楼爱情的钟 2022-05-26 06:19 273阅读 0赞

上篇博客写完之后，终于发现自己线性回归入门！

然后洗澡的时候就在想一个问题，线性回归会了，写线性拟合是完全没问题的，但是np库的多项式拟合到底是怎么做出来的呢？突然灵光一闪多项式拟合？多变量的线性回归？好像发现了什么？重新理清一下思路。什么是多项式拟合？对的，这个问题以前没有好好思考过，以前的直观感觉是，不就是给一些x，y的点，让一个多项式去拟合，使得这个多项式的曲线看起来大致符合那些点。等等，好像忽略了什么，什么叫大致符合？对的，什么叫大致符合？答案已在心中了！就是让这些点距离这个拟合曲线的和最小！慢着，这其实就是线性回归的损失函数的定义嘛！最小化，损失函数！多项式多个参数怎么确定？多变量线性回归！每个一个次方项当成独立的一个变量！比如a\*x^2+b\*x+c这个标准的二次项，x^2当成x1，x当成x2，1当成x3！好了，整个思路有了，尝试一下实现吧！

import numpy as np
    
    def my_fit(x, y, power):
        size = len(x)
        c = np.ones((1, size))
        x_arr = np.array(x)
        x_a=x_arr.copy()
        x_a.resize((1,size))
        x_mat = np.append(x_a, c, axis=0)
        y_arr = np.array(y)
        y_arr.resize((1, size))
        y_mat = np.mat(y_arr)
        for i in range(2, power+1):
            temp_x = x_arr**i
            temp_x.resize((1, size))
            x_mat = np.append(temp_x, x_mat, axis=0)
        x_mat=np.mat(x_mat)
        w = (x_mat * x_mat.T).I * x_mat * y_mat.T
        w0=w.T
        w0.resize(w0.size)
        return w0
    
    def f(x):
        return 2*x ** 2 +x+3
    
    x = np.linspace(-5, 5)
    y = f(x) + 0.5*np.random.randn(len(x))  # 加入噪音
    y_fit = np.polyfit(x, y, 2)  # 二次多项式拟合
    print(y_fit)
    w=my_fit(x,y,2)
    print(w)

激动人心的时候到了，看下输出：

[1.99205043 1.00302882 3.08335708]
    [1.99205043 1.00302882 3.08335708]

竟然和np库的结果完全一样！惊了！虽说实现代码不算优美，但是基本思路没啥区别！懒得调参用梯度下降了，直接解出来。

很好，还有一个疑问，scipy的曲线拟合怎么做的呢？一样嘛！当成多变量线性回归就好！

code：

import numpy as np
    from scipy.optimize import curve_fit
    
    
    def f_fit(x, a, b, c):
        return a*np.sin(x)+b*x+c
    
    
    def f_test(x):
        return 2*np.sin(x)+3*x+1
    
    
    def my_curve_fit1(fit_fun, x, y, p_num):  #p_num参数的个数
        size = len(x)
        if p_num <= 0 or p_num > size:
            print('no parameter to fit')
            return
        test_list = [0]*p_num
        test_list[0] = 1
        x_arr = np.array(fit_fun(x, *test_list))
        x_arr.resize((1, size))
        for i in range(1, p_num):
            test_list[i-1] = 0
            test_list[i] = 1
            temp_x = np.array(fit_fun(x, *test_list))
            temp_x.resize((1, size))
            x_arr = np.append(x_arr, temp_x, axis=0)
        x_mat = np.mat(x_arr)
        y_arr = np.array(y)
        y_arr.resize((1, size))
        y_mat = np.mat(y_arr)
        w = (x_mat * x_mat.T).I * x_mat * y_mat.T
        w0 = w.T
        w0.resize(w0.size)
        return w0
    
    #为了强行接近库函数的接口,无需参数的个数，但实现有点捞
    #很容易知道参数个数多于样本x的个数无法拟合
    def my_curve_fit2(fit_fun, x, y):
        test_list = [0]
        size = len(x)
        p_num = 1  #参数个数默认为1
        for i in range(0, size):
            try:
                fit_fun(x, *test_list)
                break
            except:
                p_num += 1
                test_list.append(0)
        if p_num > size:
            print('can not fit')
            return
        test_list[0] = 1
        x_arr = np.array(fit_fun(x, *test_list))
        x_arr.resize((1, size))
        for i in range(1, p_num):
            test_list[i - 1] = 0
            test_list[i] = 1
            temp_x = np.array(fit_fun(x, *test_list))
            temp_x.resize((1, size))
            x_arr = np.append(x_arr, temp_x, axis=0)
        x_mat = np.mat(x_arr)
        y_arr = np.array(y)
        y_arr.resize((1, size))
        y_mat = np.mat(y_arr)
        w = (x_mat * x_mat.T).I * x_mat * y_mat.T
        w0 = w.T
        w0.resize(w0.size)
        return w0
    
    
    
    
    
    x = np.linspace(-2*np.pi, 2*np.pi)
    y = f_test(x)+0.3*np.random.randn(len(x))   #加入噪音
    p_fit, prov = curve_fit(f_fit, x, y)  #曲线拟合
    my_fit1 = my_curve_fit1(f_fit, x, y, 3)
    my_fit2 = my_curve_fit2(f_fit, x, y)
    print('sicpy库的曲线拟合')
    print('a,b,c', p_fit)
    print('手写曲线拟合1')
    print('a,b,c', my_fit1)
    print('手写曲线拟合2')
    print('a,b,c', my_fit2)

输出：

sicpy库的曲线拟合
    a,b,c [1.99131057 3.00473361 1.05173153]
    手写曲线拟合1
    a,b,c [1.99131057 3.00473361 1.05173153]
    手写曲线拟合2
    a,b,c [1.99131057 3.00473361 1.05173153]

结果竟然也就是一样！厉害了，竟然一样。好了，实现的小窍门用到了列表参数传递，想单独提出来的那个项，就令该项的参数为1，其他参数为0，然后就把每个项都单独提出来。一用多变量线性回归，搞掂！然后为了强行接近scipy库的曲线拟合的接口，不知道有几个参数拟合，怎么办呢？好像没有什么好的方法，只好用try来试了，反正传的参数个数不对就会报错的，就一个一个参数往上加咯。

以上。