线性回归与曲线拟合

效果

下面是随便编造的一组数据在1~7阶上的拟合效果。

数学原理

最小二乘

假设拟合有m个样本的n维数据，设拟合线性方程如下

则实际的方程如下，

h(x)是误差项，服从

移项， 则有：

对它作极大似然估计（这是参数估计的一般步骤），似然方程如下：

因为我们只关心theta, 因此把sigma看成常数，此时当以下函数取最小值时，L有最大值：

这就是最小二乘法的由来，即用实际值减去拟合的值 的平方和表示误差，我们希望误差最小。即把问题转化成了：在theta取何值时，误差最小, 。我们把J用矩阵表示：

注意，其中X是把每个向量按行摆放，共m行，n列。

求解

为了求在theta为何值时，J取得最小值，我们需要通过求导来得到它的驻点，取得驻点时的theta就是我们要求的值。以下是求导过程：(矩阵求导自行百度)

但是，往往XTX不一定可逆，因此需要增加一个扰动，使它可逆。扰动被增加在主对角线上，有防止过拟合和特征选择的功能。以下是一种方法：（可以证明这个方法一定会使矩阵可逆， 略。除了L2正则，还有L1正则，略。）

其中必须使Lambda > 0 , 一般取一个很小的值。

好了，到这里，就把参数求出来的。但是你会发现这仅仅适用于线性的情况。但其实很容易扩展，例如

代码

以下代码是开头示例的源代码。

import numpy as np
import matplotlib.pyplot as plt

#随便造的数据
x = np.array([1, 2, 3, 4, 5, 6, 7, 8]);
y = np.array([1, 3, 7, 13, 26, 68, 150, 300]);


#计算参数函数-times为阶数
def cal(x, y, times):
    t = np.ones((times + 1, 8));
    for i in range(times):
        t[i + 1] = np.power(x, i + 1);
    
    x = t.T;
    #带入公式
    theta = np.linalg.inv((x.T.dot(x) + 0.001)).dot(x.T).dot(y)
    return theta;

#根据参数计算y值的函数-times为阶数
def getY(theta, x, times):
    t = np.ones((1, len(x)));
    for i in range(times + 1):
        t += np.power(x, i) * theta[i];
    
    return t;


#先把点画出来
plt.scatter(x, y)

a = np.linspace(1, 8, 100);

#每别拟合1~7阶
for i in range(7):
    theta = cal(x, y, i + 1);
    #print(theta);
    b = getY(theta, a, i + 1);
    plt.plot(a, b.reshape(100));

labels = np.arange(1, 8, 1);
plt.legend(labels);
plt.show()