2.3 模型拟合

一、最优化问题

在了解模型优化拟合前，我们需要了解优化的概念。

举一个简单的例子，我们一定都见过如下的函数：

$$f(x) = (x-1)^2$$

求得此函数的最小值，就是一个简单的优化问题。通过计算，我们可以得出当x=1时，此函数取得最小值。

接下来，尝试通过编程来让计算机完成这个工作：

import numpy as np
import matplotlib.pyplot as plt 
from scipy.optimize import minimize

func = lambda x: (x-1)**2 # 我们定义这个损失函数

# 用minimize去优化求解， 这里x0是参数搜索的初始值
# fun这里是minimize的一个形式参数, 我们给其赋值为上面定义的func函数。函数本身也是可以作为一个值传入到另外一个函数
res = minimize(fun=func, x0=2)
print(res.x)

结果显示，计算机给出的x = 0.999999999，这个结果非常精确但并不完全准确。这是因为在计算机的逻辑中，它是通过非常小的步长进行循环然后比较每次函数的大小，最后找出那个最小函数取值时的x值。计算机在这里的工作逻辑有点类似于我们使用枚举法，最后得到一个数值解。

在一定的约束条件下，找到目标函数的最大值或最小值的问题即为最优化问题（Optimization Problem），在计算机中，使用minimize()函数完成。其广泛应用于工程、经济学、运筹学、计算机科学等领域，用于解决资源分配、成本最小化、利润最大化等实际问题。

二、损失函数/成本函数/目标函数

在最优化问题中，很重要的就是找到最优化的目标是什么，在机器学习领域，一般会用损失函数（Loss Function）、成本函数（Cost Function）或目标函数（Objective Function）来描述即将要优化的对象。

• 损失函数指的是我用我的模型去描述数据的时候，损失了多少的准确信息。采用模型得到的预测结果来代表真实结果将造成多少的数据损失。

• 成本函数指的是用我的模型去描述真实数据时，这个想要以模型代替真实数据的代价（cost）有多高。

• 目标函数则指的是模型要最优化的目标。如果目标函数取得最值，那么这个模型就是效果最好的。

一般来说，损失函数、成本函数、目标函数在机器学习领域都有使用且内涵没有差异，但可能根据具体上下文有些许细微差别，在本书中，认为三者只是同一内容的不同表达，并且在后文中仅使用损失函数来指代这一内容。

损失函数有多种类型，最朴素的思想，就是我们可以将所有的误差（error）相加，但是因为不同的error有符号问题的存在，所以一般通过一系列数学变换使得所有的误差均为正直。

比如常用的损失函数有，均方误差（MSE），均绝对误差（MAE）等。

在求解最优化问题过程中，我们一般最小化损失函数以达到最好的模型拟合效果。

三、求解线性回归的三种方法

假设真模型为：

$$y = ax + b$$

其中a = 0.5, b = 3。我们选取高斯噪声标准差为5，通过如下代码模拟数据并绘制二维散点图：

import numpy as np
import matplotlib.pyplot as plt 
from scipy.optimize import minimize

a = 0.5 # coefficient
b = 3   # intercept
noiseVar = 1 # noise的标准差
x = np.arange(1, 100, 0.1)
# generate data
y =  a * x + b + noiseVar* np.random.randn(x.size)

print(x.size)

plt.plot(x, y, 'o', color='r')

我们会得到如下的散点图：

2.3.1 模拟数据结果

问题来了，当我只有这张散点图和数据点的资料时，我该如何求解线性回归？

方法一：手写损失函数求解

在代码编写中给出损失函数的形式，然后在后续的编程中通过计算机数值解来求得线性回归。具体操作步骤如下：

先定义此线性回归的损失函数：

def lossfun(params):
    # params是个(2,)的数组, 第一个是斜率，第二个是截距
    a = params[0]
    b = params[1]

    y_pred = a * x + b

    return ((y-y_pred)**2).sum() # 计算squared error作为损失函数

随后我们来最优化这个损失函数：

from scipy.optimize import minimize
res = minimize(fun=lossfun, x0=(1, 1))
print('Linear coefficient is', res.x[0])
print('Intercept is', res.x[1])

结果显示a = 0.5009765933689194, b = 2.9547084373733297。这个数值解与我们的真模型的参数已经非常相近了。

方法二：利用线性回归的解析解求解。

解析解可以理解为我们利用求解一元二次方程的公式来得到自变量结果。在线性回归中我们同样可以利用公式来直接求解最优值。

在这个线性回归中，可以用如下最小二次法求解公式计算解析解：

$$\begin{aligned} &\mathbf{X}^\top \left( \mathbf{y} - \mathbf{X} \hat{\mathbf{w}}^{\mathrm{OLS}} \right) = 0 \\ &\mathbf{X}^\top \mathbf{y} - \mathbf{X}^\top \mathbf{X} \hat{\mathbf{w}}^{\mathrm{OLS}} = 0 \\ &\hat{\mathbf{w}}^{\mathrm{OLS}} = \left( \mathbf{X}^\top \mathbf{X} \right)^{-1} \mathbf{X}^\top \mathbf{y} \end{aligned}$$

具体代码操作步骤如下：

x2 = np.vstack((x, np.ones(x.size))) # x2的第二排都是1
# 因为numpy都是横向量，我们划成列向量
x2 = x2.T # x2 is 99 x 2
y_hat = y[:, np.newaxis].copy() # y is 99 x 1

res = np.linalg.inv(x2.T @ x2) @ x2.T @ y_hat 

print('Coefficient is ',res[0][0])
print('Intercept is ',res[1][0])

结果显示a = 0.5009766233665898, b = 2.954706934764938。可以看到，这个结果同样也是比较精确的。

方法三：利用算法包直接求解。

有些算法包已经提供了对于特定模型的专用优化函数，我们可以直接调用这些函数来得到最优解。

具体操作步骤如下：

x2 = np.vstack((x, np.ones(x.size))) # x2的第二排都是1

# 利用lstsq函数求解
res = np.linalg.lstsq(x2.T, y, rcond=None)[0]

print('Coefficient is ',res[0])
print('Intercept is ',res[1])

结果显示a = 0.5009766233665904, b = 2.9547069347648884。同样得到了较为不错的结果。