> For the complete documentation index, see [llms.txt](https://ruyuanzhang.gitbook.io/compmodcogpsy/llms.txt). Markdown versions of documentation pages are available by appending `.md` to page URLs; this page is available as [Markdown](https://ruyuanzhang.gitbook.io/compmodcogpsy/di-er-zhang-ji-suan-mo-xing-ji-chu/2.3-mo-xing-ni-he.md).

# 2.3 模型拟合

Edited by 陆新泉

如果我们已经有了一个模型，无论是线性还是非线性，下一步是如何拟合该模型，在了解模型拟合(model fitting)之前，我们先聊聊优化的问题。

## 优化问题

心理学的学生可能对于优化(optimization)这件事并不熟悉。优化这个问题是理论机器学习里面的重要问题，很多高校的人工智能课程也开发了相应课程专门讲优化。这里我们只是浅显的介绍一下优化的概念。

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

$$
f(x) = (x-1)^2\qquad (2-6)
$$

一个简单的问题，当$$x$$等于多少的时候，$$f(x)$$的值最小？这是一个小学生都能做的问题。但是如果要计算机来解这个问题，要麻烦很多，因为很多时候计算机不能很轻易的完成解析解，必须要通过数值优化的办法。

{% hint style="info" %}
解析解：解析解是通过解析方法得到的精确解。解析方法是指能够通过数学推导和运用已知的数学定理和公式，直接求得问题的解的方法。

数值解：数值解是通过数值计算方法近似求解数学问题的方法得到的解。数值计算方法将问题转化为数值计算的形式，并使用计算机进行近似计算。
{% endhint %}

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

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

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

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

我们仔细探究一下`res=minimize(fun=func, x0=2)` 这行代码。`minimize`这个函数是为了一个固定函数的极小值，翻译成大白话就是: 请帮我找一个自变量$$x$$的值，使得`func` 这个函数取值最小，并且从2开始搜索这个值。

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

在一定的约束条件下，找到目标函数的最大值或最小值的问题，即为最优化问题（Optimization Problem）。各个软件或者编程语言(例如matlab/python)都有自己的优化函数功能。本书中，我们全部使用python `scipy` 包里面的`minimize`函数完成。

## 损失函数

在计算建模中，很重要的就是找到最优化的目标函数是什么，在机器学习领域，一般会用损失函数(Loss Function），有时候也叫Cost Function或者Objective Function来描述即将要优化的对象。

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

例如，在上面求解的问题中，$$f(x)=(x-1)^2$$就是损失函数，我们要优化的就是这个函数。

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

比如常用的损失函数有, 均方误差(mean squared error)和均绝对误差(mean absolute error)等。

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

## 线性回归问题的损失函数

我们回到线性回归问题，我们可以定义最简单的线性模型:

$$
y = ax + b\qquad (2-7)
$$

我们通常认为线性回归问题的误差来自高斯噪声，那么损失函数就是最小二乘的形式。

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

假设真模型为：

$$
y = ax + b\qquad (2-7)
$$

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

```python
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')
plt.show()
```

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

<figure><img src="https://1379976374-files.gitbook.io/~/files/v0/b/gitbook-x-prod.appspot.com/o/spaces%2Fu8x1pCBjIDBIizdIV9Wv%2Fuploads%2FvY47KJBpaEg9iQb1XDdM%2Fimage.png?alt=media&amp;token=92ae5b51-6d98-4b02-b162-1a8d7e8eae3b" alt=""><figcaption><p>图2-5 模拟数据结果</p></figcaption></figure>

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

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

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

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

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

    y_pred = a * x + b # 计算预测的因变量值

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

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

```python
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}\qquad (2-8)
$$

具体代码操作步骤如下：

```python
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。可以看到，这个结果同样也是比较精确的。

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

我们还可以利用Python里面一些已经有的解决最小二乘法的函数来求解参数。例如，在python里面的numpy包中就有`numpy.linalg.lstsq` 这个函数用于求解最小二乘法。

具体操作步骤如下：

```python
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。同样得到了较为不错的结果。