机器学习实验：回归 | Leeml-2020-hw1

Overview

实验已作为 Kaggle 竞赛发布，见 Kaggle。涉及的主要内容：

• 通过完成一次时间序列预测任务，回顾回归相关知识
• 梯度下降的步骤，求解方法
• 学习率参数选择
• 特征缩放对学习效果的影响
• 最小二乘法直接求解回归问题

EDA

任务说明：训练集数据提供丰原气象站某年每月前 20 天每小时的 18 项指标数值。测试集给出余下天数的前 9 小时的所有指标数值，要求预测出第 10 小时的 PM2.5 数值。

load data

TASK_DIR = ROOT_DIR + 'ml2020spring-hw1/'
data_train = pd.read_csv(TASK_DIR + "train.csv",encoding = 'big5')
data_test  = pd.read_csv(TASK_DIR + "test.csv" ,encoding = 'big5',names = ['日期', '測項', '0', '1', '2', '3', '4', '5', '6', '7', '8'])

data_train.info()

<class 'pandas.core.frame.DataFrame'>
RangeIndex: 4320 entries, 0 to 4319
Data columns (total 27 columns):
 #   Column  Non-Null Count  Dtype 
---  ------  --------------  ----- 
 0   日期      4320 non-null   object
 1   測站      4320 non-null   object
 2   測項      4320 non-null   object
 3   0       4320 non-null   object
 4   1       4320 non-null   object
 5   2       4320 non-null   object
 6   3       4320 non-null   object
 7   4       4320 non-null   object
 8   5       4320 non-null   object
 9   6       4320 non-null   object
 10  7       4320 non-null   object
 11  8       4320 non-null   object
 12  9       4320 non-null   object
 13  10      4320 non-null   object
 14  11      4320 non-null   object
 15  12      4320 non-null   object
 16  13      4320 non-null   object
 17  14      4320 non-null   object
 18  15      4320 non-null   object
 19  16      4320 non-null   object
 20  17      4320 non-null   object
 21  18      4320 non-null   object
 22  19      4320 non-null   object
 23  20      4320 non-null   object
 24  21      4320 non-null   object
 25  22      4320 non-null   object
 26  23      4320 non-null   object
dtypes: object(27)
memory usage: 911.4+ KB

data_train.head(20)

	日期	測站	測項	0	1	2	3	4	5	6	...	14	15	16	17	18	19	20	21	22	23
0	2014/1/1	豐原	AMB_TEMP	14	14	14	13	12	12	12	...	22	22	21	19	17	16	15	15	15	15
1	2014/1/1	豐原	CH4	1.8	1.8	1.8	1.8	1.8	1.8	1.8	...	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8
2	2014/1/1	豐原	CO	0.51	0.41	0.39	0.37	0.35	0.3	0.37	...	0.37	0.37	0.47	0.69	0.56	0.45	0.38	0.35	0.36	0.32
3	2014/1/1	豐原	NMHC	0.2	0.15	0.13	0.12	0.11	0.06	0.1	...	0.1	0.13	0.14	0.23	0.18	0.12	0.1	0.09	0.1	0.08
4	2014/1/1	豐原	NO	0.9	0.6	0.5	1.7	1.8	1.5	1.9	...	2.5	2.2	2.5	2.3	2.1	1.9	1.5	1.6	1.8	1.5
5	2014/1/1	豐原	NO2	16	9.2	8.2	6.9	6.8	3.8	6.9	...	11	11	22	28	19	12	8.1	7	6.9	6
6	2014/1/1	豐原	NOx	17	9.8	8.7	8.6	8.5	5.3	8.8	...	14	13	25	30	21	13	9.7	8.6	8.7	7.5
7	2014/1/1	豐原	O3	16	30	27	23	24	28	24	...	65	64	51	34	33	34	37	38	38	36
8	2014/1/1	豐原	PM10	56	50	48	35	25	12	4	...	52	51	66	85	85	63	46	36	42	42
9	2014/1/1	豐原	PM2.5	26	39	36	35	31	28	25	...	36	45	42	49	45	44	41	30	24	13
10	2014/1/1	豐原	RAINFALL	NR	NR	NR	NR	NR	NR	NR	...	NR	NR	NR	NR	NR	NR	NR	NR	NR	NR
11	2014/1/1	豐原	RH	77	68	67	74	72	73	74	...	47	49	56	67	72	69	70	70	70	69
12	2014/1/1	豐原	SO2	1.8	2	1.7	1.6	1.9	1.4	1.5	...	3.9	4.4	9.9	5.1	3.4	2.3	2	1.9	1.9	1.9
13	2014/1/1	豐原	THC	2	2	2	1.9	1.9	1.8	1.9	...	1.9	1.9	1.9	2.1	2	1.9	1.9	1.9	1.9	1.9
14	2014/1/1	豐原	WD_HR	37	80	57	76	110	106	101	...	307	304	307	124	118	121	113	112	106	110
15	2014/1/1	豐原	WIND_DIREC	35	79	2.4	55	94	116	106	...	313	305	291	124	119	118	114	108	102	111
16	2014/1/1	豐原	WIND_SPEED	1.4	1.8	1	0.6	1.7	2.5	2.5	...	2.5	2.2	1.4	2.2	2.8	3	2.6	2.7	2.1	2.1
17	2014/1/1	豐原	WS_HR	0.5	0.9	0.6	0.3	0.6	1.9	2	...	2.1	2.1	1.9	1	2.5	2.5	2.8	2.6	2.4	2.3
18	2014/1/2	豐原	AMB_TEMP	16	15	15	14	14	15	16	...	24	24	23	21	20	19	18	18	18	18
19	2014/1/2	豐原	CH4	1.8	1.8	1.8	1.8	1.8	1.8	1.8	...	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8	1.8

20 rows × 27 columns

训练集由每月前 20 天每小时的 18 项观测指标组成，除 RAINFALL 观测指标外其余指标均为数值型。下载观测整个数据集发现，RAINFALL 指标值不为 NR 时，以数值型记录了该小时的降雨量，所以需要对 NR 值进行处理，将其改为 0。

data_test

	日期	測項	0	1	2	3	4	5	6	7	8
0	id_0	AMB_TEMP	21	21	20	20	19	19	19	18	17
1	id_0	CH4	1.7	1.7	1.7	1.7	1.7	1.7	1.7	1.7	1.8
2	id_0	CO	0.39	0.36	0.36	0.4	0.53	0.55	0.34	0.31	0.23
3	id_0	NMHC	0.16	0.24	0.22	0.27	0.27	0.26	0.27	0.29	0.1
4	id_0	NO	1.3	1.3	1.3	1.3	1.4	1.6	1.2	1.1	0.9
...	...	...	...	...	...	...	...	...	...	...	...
4315	id_239	THC	1.8	1.8	1.8	1.8	1.7	1.7	1.7	1.7	1.7
4316	id_239	WD_HR	80	92	95	95	96	97	96	96	84
4317	id_239	WIND_DIREC	76	99	93	97	93	94	98	97	65
4318	id_239	WIND_SPEED	2.2	3.2	2.5	3.6	5	4.2	5.7	4.9	3.6
4319	id_239	WS_HR	1.7	2.8	2.6	3.3	3.5	5	4.9	5.2	3.6

4320 rows × 11 columns

根据任务描述，训练集记录的是每月前 20 天的数据，那么每月余下的数据便打乱后包含于测试集，同时只保留前 9 小时的数据，目标是根据这 9 个小时的数据预测出第 10 个小时的 PM 2.5 值。同理，测试集的 NR 稍后也需处理为 0。

data_train = data_train.replace('NR',0.0)
data_test  = data_test.replace('NR',0.0)

我们的任务是：根据训练集数据，使用回归的方法，对模型进行训练，使得模型能够根据测试集中前 9 小时的数据预测出第 10 小时的 PM2.5 值。

由于涉及的项较多，为简化任务，本次实验的模型均采用线性模型，那么考虑模型复杂程度可以分为以下两种模型：

1. 仅考虑 PM2.5 指标，根据前 9 小时的 PM2.5 指标预测第 10 小时的值；
2. 考虑所有指标对 PM2.5 的影响，根据前 9 小时所有指标预测第 10 小时的值。

测试集中，前 9 小时的数据是连续的 9 小时，所以可以从训练集中提取连续的 9 小时作为训练样本，之后的 1 小时的 PM2.5 观测值作为标签，获得训练集。

data_train.iloc[:,3:] = data_train.iloc[:,3:].astype(float)

X_train = []
y_train = []
x_train_pm25 = []
for date in data_train.日期.unique():
    data_daily = data_train[data_train['日期']==date]
#     display(data_daily)
    properties = data_daily.測項
#     display(properties)
    for i in range(0+3,16+3):
        X_train.append(data_daily.iloc[:,i:(i+9)].values)
        x_train_pm25.append(data_daily.iloc[9,i:(i+9)].values)
        y_train.append(data_daily[data_daily['測項']=='PM2.5'].iloc[:,(i+8)].values[0])

就这样，我们获取到了基础的训练数据集，可以开始后续的模型训练。

Regression - Model1 : Simple Linear Model

Definition

上面的分析已经提到，我们的模型为线性模型。以最简单的模型 1 为例，最后的模型应为如下形式：

$$f(\boldsymbol{x}) = w_1 \cdot x_1 + w_2 \cdot x_2 + \cdots + w_9 \cdot x_9 + b$$

其中 $y$ 为第 10 小时 PM2.5 预测值，$x_1,x_2,\cdots,x_9$ 为前 9 小时的 PM2.5 观测值。

为了后续讨论方便，尽可能将模型的元素使用向量或矩阵的形式表示，原模型可写为：

$$f(\boldsymbol{x}) = \boldsymbol{w}^T\boldsymbol{x}+b$$

为便于后续讨论，将模型进一步简化表示，将偏置 $b$ 收入向量形式 $\boldsymbol{\hat{w}}=(\boldsymbol{w},b)$。全体数据集表示为一个矩阵，行列数量对应样本数、属性值数+1，其中每行对应一个样本，最后一个元素恒置为 1：

$$\boldsymbol{X} = \left(\begin{matrix} x_{11} & x_{12} & \cdots & x_{1d} & 1 \\ x_{21} & x_{22} & \cdots & x_{2d} & 1 \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ x_{m1} & x_{m2} & \cdots & x_{md} & 1 \end{matrix}\right) = \left(\begin{matrix} \boldsymbol{x}_1^T & 1 \\ \boldsymbol{x}_2^T & 1 \\ \vdots & \vdots \\ \boldsymbol{x}_m^T & 1 \\ \end{matrix}\right)$$

其中 $d=9$。

我们的目标是，学习得到 $f(\boldsymbol{x}_i)$ 使得 $f(\boldsymbol{x}_i)\simeq y_i$

Loss Function

模型在不断优化的过程中，需要对模型效果进行评估。在回归任务中，均方误差（MSE）是最常见的性能度量。

采用均方误差作为损失函数，我们原来学习函数的问题转化为求取所有样本点到我们拟合的模型曲线上的欧式距离之和，并迭代优化模型使之最小。

为了方便表示，将样本点标签（label，即我们预测的第 10 小时 PM2.5 值）也写成向量形式：$\boldsymbol{y} = (y_1;y_2;\cdots;y_m)$，损失函数可以表示为：

$$L = \dfrac{1}{m}(\boldsymbol{y} - \boldsymbol{X\hat{w}})^T(\boldsymbol{y} - \boldsymbol{X\hat{w}})$$

Gradient Descent

随机下降过程是利用样本学习参数的过程，将参数向梯度下降方向（损失函数值在某一点减少的方向）移动。

梯度下降可以分为以下步骤：

1. 随机选取一个初始点 $\boldsymbol{\hat{w}}^0$ ；
2. 计算在该处参数 $\boldsymbol{\hat{w}}$ 对损失函数 $L$ 的（偏）微分 $\dfrac{\partial L}{\partial \boldsymbol{\hat{w}}}|_{\boldsymbol{\hat{w}}=\boldsymbol{\hat{w}}^i}$；
3. 更新 $\boldsymbol{\hat{w}}^i$：$\boldsymbol{\hat{w}}^{i+1} \leftarrow \boldsymbol{\hat{w}}^i - \eta\dfrac{\partial L}{\partial \boldsymbol{\hat{w}}}|_{\boldsymbol{\hat{w}}=\boldsymbol{\hat{w}}^i}$。其中 $\eta$ 为学习率，微分值为负则增加参数，若微分值为正则减少参数。重复步骤2-3，直至微分值为 $0$，得到 $\boldsymbol{\hat{w}}$ 的最优值 $\boldsymbol{\hat{w}}^\star$ 。

梯度下降较为关键的步骤就是梯度（损失函数偏微分）的计算和参数更新，计算梯度的方法与损失函数的选择相关，前面我们选择了 MSE 作为损失函数。$L$ 对 $\boldsymbol{\hat{w}}$ 求偏导得：

$$\dfrac{\partial L}{\partial \boldsymbol{\hat{w}}} = \dfrac{2}{m} \boldsymbol{X}^T(\boldsymbol{X}\boldsymbol{\hat{w}}-\boldsymbol{y})$$

def grad_descent(X,Y,epoch=20,eta=0.0001,w_old=None):
    """
    Args:
        X: 输入数据集, Numpy 数组 (Sample_count, Attr_count)
        Y: 预测标签, Numpy 数组 (Sample_count, 1)
        epoch: 迭代轮数
        eta: 学习率
    """
    X = np.column_stack((X,np.ones(X.shape[0])))
    Y = Y.reshape((Y.shape[0],1))
    w = w_old if w_old is not None else np.random.random(size=(X.shape[1],1)) # 随机初始化参数 w（Step 1）
    loss = []
    for iter in range(epoch):
        print("\r{}/{}".format(iter+1,epoch),end="")
        gradient = 2 * X.T.dot(X.dot(w)-Y) / X.shape[0]                       # 计算 L 对 w 的偏微分，即梯度。（Step 2）
        w = w - eta * gradient                                                # 更新参数 w （Step 3）
        loss.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)/ X.shape[0])   # 计算损失 L
    return w, loss

w,loss = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float))

20/20

通过迭代我们学习得到了模型参数 $\boldsymbol{\hat{w}}$ 和每一轮函数的损失值，我们可以通过可视化的方式直观展现函数损失

def plot_loss(loss,title="Loss Value per Iter"):
    fig = plt.figure()
    ax = fig.add_subplot()
    ax.set_title(title)
    ax.set_xlabel("iteration")
    ax.set_ylabel("Loss/MSE")
    ax.plot(range(len(loss)),np.array(loss).reshape(-1))
plot_loss(loss)
print(loss[-1])

[26.05574238]

可以看到，随着训练轮数的增加，损失值在不断下降，说明我们的模型在不断学习，向数据逼近。

参数调整：训练轮数 epoch

如果我们进一步增加训练轮数，可以使得损失函数进一步收敛：

w,loss = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float),epoch=80,w_old=w)
print()
plot_loss(loss)
print(loss[-1])

80/80
[4.08627535]

函数收敛的速度不断降低，所以如果无限度增大训练轮数确实可以使得损失函数较小，但是有可能带来过拟合现象。所以调整训练轮数在适当值即可。

参数调整：学习率 eta

学习率能够控制参数更新速度，过大的学习率可能导致损失函数错过全局最优，甚至无法收敛；而较小的学习率可能导致损失函数下降较慢

由于我们的参数矩阵是随机生成，为便于对比不同学习率的效果，我们控制每一次的初始参数都相同

w = np.random.random(size=(len(x_train_pm25[0])+1,1))
w1,loss1 = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float),epoch=50,eta=0.0001,w_old=w)
w2,loss2 = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float),epoch=50,eta=0.001,w_old=w)
w3,loss3 = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float),epoch=50,eta=0.00001,w_old=w)

50/50

plot_loss(loss1,title="Loss per Iter: $\eta=10^{-4}$")
print(loss1[-1])
plot_loss(loss2,title="Loss per Iter: $\eta=10^{-3}$")
print(loss2[-1])
plot_loss(loss3,title="Loss per Iter: $\eta=10^{-6}$")
print(loss3[-1])

[11.36702501]
[8.90902162e+109]
[37.52494183]

可以看到，当 $\eta=0.001$ 时，损失函数最终飙升，无法收敛，此时学习率过大；当 $\eta=0.000001$ 时，损失函数最初的收敛速度明显不如 $\eta=0.0001$ 时，且经过同样较大的迭代轮数后，损失函数值的收敛效果不如后者，学习率较小。

所以，需要经过多次测试，比较模型在不同学习率下的表现，方能获得较好的学习效果。

模型预测

通过参数调整测试，我们已经确定了较好的参数，这样就可以获取最终的模型，对目标数据进行预测。

$$f(\boldsymbol{X}) = \boldsymbol{\hat{w}X}$$

w,loss = grad_descent(np.array(x_train_pm25,dtype=float),np.array(y_train,dtype=float),eta=0.00013,epoch=300)
print("Final Loss: {}".format(loss[-1]))

300/300Final Loss: [0.639463]

X_test_pm25 = data_test[data_test['測項']=='PM2.5'].iloc[:,2:].values.astype(float)
X = np.column_stack((X_test_pm25,np.ones(X_test_pm25.shape[0])))
Y = X.dot(w)

result = pd.merge(data_test[data_test['測項']=='PM2.5'].iloc[:,0].reset_index(drop=True),pd.DataFrame(Y),left_index=True,right_index=True,how='outer')
result.columns = ['id','value']
result.to_csv(OUT_DIR + "submission.csv",index=False)

结果

最终的 Public Score 是 6.28859（击败 Simple Baseline: 6.55912），Private Score 是 8.53027（击败 Simple Baseline: 8.73773）

Regression - Model2 : Take more feature into consideration

模型构建

模型方面其实没有变化，但是这一步中，我们将考虑数据中给出的所有特征，并不局限于 PM 2.5 数据。

def grad_descent(X,Y,epoch=20,eta=0.0001,w_old=None):
    """
    Args:
        X: 输入数据集, Numpy 矩阵 (Sample_count, Attr_count)
        Y: 预测标签, Numpy 矩阵 (Sample_count, 1)
        epoch: 迭代轮数
        eta: 学习率
    """
    X = np.column_stack((X,np.ones(X.shape[0])))
    Y = Y.reshape((Y.shape[0],1))
    w = w_old if w_old is not None else np.random.random(size=(X.shape[1],1)) # 随机初始化参数 w（Step 1）
    loss = []
    for iter in range(epoch):
        print("\r{}/{}".format(iter+1,epoch),end="")
        gradient = 2 * X.T.dot(X.dot(w)-Y) / X.shape[0]                       # 计算 L 对 w 的偏微分，即梯度。（Step 2）
        w = w - eta * gradient                                                # 更新参数 w （Step 3）
        loss.append((Y-X.dot(w)).T.dot(Y-X.dot(w)).reshape(-1)/ X.shape[0])   # 计算损失 L
    
    return w, loss

def plot_loss(loss,title="Loss Value per Iter"):
    fig = plt.figure()
    ax = fig.add_subplot()
    ax.set_title(title)
    ax.set_xlabel("iteration")
    ax.set_ylabel("Loss/MSE")
    ax.plot(range(len(loss)),np.array(loss).reshape(-1))

w,loss = grad_descent(np.array(X_train,dtype=float).reshape(len(X_train),-1),np.array(y_train,dtype=float),epoch=300,eta=0.0000014)

300/300

plot_loss(loss)

在这部分，我们略去 Model1 中已经涉及的参数调节过程，直接选取我多次实验后选取的最优参数，随后直接进行预测。

模型预测

X_test = data_test.iloc[:,2:].values.astype(float)
id = data_test.日期.unique()
X_test = X_test.reshape(len(id),-1)

X = np.column_stack((X_test,np.ones(X_test.shape[0])))
Y = X.dot(w)
result = pd.merge(pd.DataFrame(id),pd.DataFrame(Y),left_index=True,right_index=True,how='outer')
result.columns = ['id','value']
result.to_csv(OUT_DIR + "submission.csv",index=False)

结果

确定 $\eta = 1.4 \times 10^{-6}$ 是学习率的最优参数后，在 Kaggle 上进行了多次提交，发现设定迭代次数在 50000 次以内时，结果得分在 10 ~ 30；而设定迭代次数到 100000 次时， 取得了 Public Score 为 6.406，Private Score 为 8.477，表现优于 Model1 ，但是显然出现了 参数收敛较慢 的问题。

考虑更多指标，本质上是增加了模型的复杂度，因为 Model2 的表示能力是大于且包含 Model1 的表示能力的（因为 Model1 的 9 个参数均是 Model2 的参数），所以 Model2 的理论性能是一定要高于 Model1 的。

如果提升学习率，会使得损失函数无法收敛，所以在其他方面对模型还有优化空间。

Regression - Model3 : Normalize the Feature

在 Model2 中考虑了更多指标，提升了模型表现，但是带来了另一个问题：参数收敛较慢，造成该现象的原因可能是，我们引入了新的参数，而新的参数与 Model1 的特征的量纲完全不同，即我们只考虑 PM2.5 指标值，它的范围可能是 0 - 200，而如果考虑 WIND_DIRECTION，它的范围则是 0 ~ 360。对于不同范围的特征，参数的敏感程度是不同的，就可能带来不同参数收敛速度不同，造成整体参数收敛速度较慢。而如果将所有参数都标准化，就可以在一定程度上提升参数的收敛速度。

Normalizaiton

我们将采用 Min-Max Normalization，需要获取各参数的最大值、最小值。

X_large = np.array(X_train[0])
for i in range(len(X_train)):
    X_large = np.column_stack((X_large,X_train[i]))
X_max = np.max(X_large,axis=1)
X_min = np.min(X_large,axis=1)

随后，对数据中的每一个元素都进行标准化操作。Min-Max Normalization 的方法是：

$$x^\star = \dfrac{x-\min{\boldsymbol({x})}}{\max{\boldsymbol({x})}-\min{\boldsymbol({x})}}$$

X_norm = []
for i in range(len(X_train)):
    X_norm.append(((X_train[i].T-X_min)/(X_max-X_min)).T)

模型构建

w,loss = grad_descent(np.array(X_norm,dtype=float).reshape(len(X_norm),-1),np.array(y_train,dtype=float),epoch=1000,eta=0.035)

1000/1000

plot_loss(loss)
print(loss[-1])

[20.63629716]

可以看到，对数据进行标准化操作后，参数收敛的速度没有比上一轮快，但是损失值远远比不标准化特征还要低，这可能是因为特征值被放缩导致欧式距离减小导致的。但是相比而言，损失值比不放缩数据更早收敛。

模型预测

由于我们对训练集输入进行了标准化，那么我们对测试集输入也需要进行相同的操作，而这里标准化的最大值、最小值并不是取决于测试集，而是取决于训练集，这样才能使得两个数据集具有对等的放缩映射，模型参数才能生效。

X_test = data_test.iloc[:,2:].values.astype(float)
id = data_test.日期.unique()
attr_count = int(X_test.shape[0]/len(id))
X_test_list = []
for i in range(int(X_test.shape[0]/attr_count)):
    X_i = X_test[i*attr_count:(i+1)*attr_count]
    X_test_list.append(((X_i.T-X_min)/(X_max-X_min)).T)

X_test = np.array(X_test_list).reshape(int(X_test.shape[0]/attr_count),-1)
X = np.column_stack((X_test,np.ones(X_test.shape[0])))
Y = X.dot(w)
result = pd.merge(pd.DataFrame(id),pd.DataFrame(Y),left_index=True,right_index=True,how='outer')
result.columns = ['id','value']
result.to_csv(OUT_DIR + "submission.csv",index=False)

结果

迭代同样次数情况下，取得 Public Score 6.286，Private Score 8.221 成绩。

Regression - Model4 : Ordinary Least Squares

如果数据量较大，能够实现输入矩阵求逆，就可以不使用梯度下降的方法逐次迭代寻找最优，而是利用最小二乘法直接求取最优值。

模型构建

前面我们已经得到

$$\dfrac{\partial L}{\partial \boldsymbol{\hat{w}}} = \dfrac{2}{m} \boldsymbol{X}^T(\boldsymbol{X}\boldsymbol{\hat{w}}-\boldsymbol{y})$$

我们可以令 $\dfrac{\partial L}{\partial \boldsymbol{\hat{w}}}=0$ ，即可得到最优解

$$\boldsymbol{\hat{w}}^\star = (\boldsymbol{X}^T\boldsymbol{X})^{-1}\boldsymbol{X}^T\boldsymbol{y}$$

def ols(X,Y):
    X = np.column_stack((X,np.ones(X.shape[0])))
    Y = Y.reshape((Y.shape[0],1))
    w = np.linalg.inv(X.T.dot(X)).dot(X.T.dot(Y))
    return w

w = ols(np.array(X_norm,dtype=float).reshape(len(X_norm),-1),np.array(y_train,dtype=float))

模型预测

Y = X.dot(w)
result = pd.merge(pd.DataFrame(id),pd.DataFrame(Y),left_index=True,right_index=True,how='outer')
result.columns = ['id','value']
result.to_csv(OUT_DIR + "submission.csv",index=False)

结果

采用最小二乘法直接求取闭式解的结果是：Public Score 6.290，Private Score 8.221，表现均比 Model2、3 好。

Summary

Results

本次竞赛就丰原站气象数据进行了回归预测，对模型进行了优化，最终表现总结如下表：

Model	iteration	parameters & Description	PubScore	PrivScore
GD Model 1	$10^2$	$\eta = 1.5 \times 10^{-4}$ ，仅考虑 PM2.5 指标	6.289	8.530
GD Model 2	$10^5$	$\eta = 1.4 \times 10^{-6}$ ，考虑所有提供的指标	6.406	8.422
GD Model 3	$10^5$	$\eta = 3.5 \times 10^{-2}$ ，对所有指标标准化	6.286	8.221
OLS Model	1	直接求解	6.290	8.220
Simple Baseline	Unknown	Unknown	6.559	8.738
Strong Baseline	Unknown	Unknown	5.560	7.142

尝试的所有模型均突破了 Simple Baseline，但是距离 Strong Baseline 仍然还有一定距离。

Highlights

• 相比于大部分他人做法，我采样获得的样本较多，获取了大多数连续的 8 小时作为样本，总计 3 840 样本。
• 考虑了学习率、迭代轮数、特征选择、特征放缩等因素优化模型。

Lowlights

• 可以使用交叉验证进行模型选择
• 没有进行本地测试准确率评估性能，仅参考 Loss
• 考虑的样本、特征过多可能带来了噪声，可以考虑进行降维或进一步特征选择
• 可以将学习率优化为自适应