DDL挑战——时间压力下的交作业之道

背景

　　最近我在上一门名为“地球物理观测解析学”的课程。该课程主要教授学生如何使用Python对地球物理观测结果进行分析。在探讨非线性最小二乘法时，老师按照惯例布置了一项Python作业。然而，这次的作业有点特殊——它涉及的数据集并非是地球物理学的观测数据，而是关于该课程学生提交作业数量（y）与距离提交截止的时间之间（x）的关系（图1）。

　　乍一看，只需要使用一个指数函数模型y(x) = A \exp{(\beta x)}就能轻松解决掉这次作业。然而，在不断试错后，我发现在这组数据集背后蕴含着一些更为深奥的秘密。

最初设想的几种模型

　　根据图1，很容易就能发现这次涉及到的观测量是随着时间单调递增的。因此，如表1所示，我引用了几个在计算机科学和物理学中经常用到的模型。不过，为了彰显出单调递增函数模型的适配性，我还特别增加了一个2次函数模型作为参考。

表1 最初设想的几种模型的参数和表达式

　　最终，经过迭代10次的回归分析，我得出了如表3所示的回归结果。

表3 最初设想的几种模型的回归结果

　　将这些回归结果与原始观测数据绘制成图表，便可得到图2。从图2可以看出，由于2次函数并非在定义域上一直呈现单调递增的性质，因此最终得到的回归结果与实际情况相差较大。而其他的3个模型由于都具有“先缓后急”、定义域上单调递增的性质，因此最终的回归结果基本与实际观测结果相吻合。

　　事实上，基本上所有同学在完成了上述回归分析后，便将结果总结成了PDF文档将作业提交了上去。然而，如果稍微思考一下，我们便能联想到一个极为重要的事实——通常情况下，绝大多数学生不太可能在凌晨、午休或者正常上课时间提交作业。对于这种情况的思考，让我觉得有必要继续深入展开探索。

考虑了1日周期的指数函数模型

　　如表1所示，由于指数函数模型和观测数据之间的残差范数||\bm{d}||最小，因此我决定在指数函数模型的基础上进一步展开调查。考虑到上一章末尾中提到的情况，我决定为指数函数模型只添加一个周期为1日的修正项。为了减少计算量，该修正项不考虑平日和周末的区别，也不考虑作业提交截止时间到来前为学生们施加的“紧张感”。

　　正如上文所说，由于观测量是随着时间单调递增的，因此周期修正项也应当具备这种性质。在考察了多种模型后，我决定在指数函数模型中添加如下的周期修正项。

\Psi (t) \equiv\gamma \int_0^t \sin^2(\pi t' + \varphi) dt'

　　将周期修正项\Psi(t)进行绘图，我们便得到了如图3所示的函数图像。从图3可以看出，该周期修正项长得有点像是阶梯函数，确实具备单调递增、随1天周期发生变化等性质。

　　综上所述，考虑了1日周期的指数函数模型共有5个参数（\alpha, \beta, \gamma, \varphi, \delta），其表达式可写为

n(t) = \alpha \exp (\beta t) + \gamma \int_0^t \sin^2(\pi t' + \varphi) dt' + \delta

　　在该模型的基础上进行迭代10次的回归分析，我们得到了如下回归结果。

n(t) = 47.42\exp (1.66 t) + 4.13 \int_0^t \sin^2(\pi t' -0.88) dt' + 22.26

　　该回归曲线的残差仅有16.77，与一般的指数函数模型相比减少了15%左右。根据其初始相位\varphi = -0.88，我们还可以推测出在每天24时的\left| \frac{-0.88}{2\pi} \right| = 3.35小时前，即每天20时39分左右提交作业的人数是最多的。

　　如图4所示，将回归曲线和观测结果绘制成图像，通过与图2中所示的4种模型进行比较，我们可以发现考虑了1日周期的指数函数模型更加贴合观测数据，并且在一定程度上展现出了周期变化。

总结

　　本文通过对“地球物理观测解析学”课程中作业提交量与距离提交截止的时间进行非线性回归分析，探讨了不同数学模型的适用性。结果表明，指数函数模型与实际情况最为相符。在该模型的基础上考虑周期为1日的修正项后，可以发现修正后的模型更加贴合观测数据，并且在一定程度上展现出了周期变化。

附：原理

[buy]以下内容适合对高等数学和回归分析等主题有一定了解的读者阅读。[/buy]

　　根据数据集的特点，我们可以将提交作业数n视为距离提交截止时间t的函数n(t)。由于该数据集共有66组数据，因此观测方程可写为

d_j = n_j^{\mathrm{obs}} - n_j \:(j = 1, 2, \dots, 66)

　　将\alpha_i\:(i=1, 2, \dots, n)视为n(t)的参数，利用泰勒展开，可以得到如下的近似式。

n(\alpha_1 + \Delta \alpha_1, \alpha_2 + \Delta \alpha_2, \dots, \alpha_n + \Delta \alpha_n)

  \cong n(\alpha_1, \alpha_2, \dots, \alpha_n)

  + \sum^n_{i=1} \frac{\partial n}{\partial \alpha_i}\Delta \alpha_i

　　因此，关于正规方程式\bm{d} = \bm{y} - A\bm{x}，则有

\bm{d} = \begin{pmatrix}

    d_1 \\ \vdots \\ d_{66} \\\end{pmatrix}
\\

  \bm{y} =

  \begin{pmatrix}

    n_1^{\mathrm{obs}} - n_1 \\ \vdots \\ n_{66}^{\mathrm{obs}} - n_{66} \\\end{pmatrix}
\\

  A =

  \begin{pmatrix}

    \frac{\partial n}{\partial \alpha_1}|_{t=t_1} &

    \frac{\partial n}{\partial \alpha_2}|_{t=t_1} & \dots &

    \frac{\partial n}{\partial \alpha_n}|_{t=t_1} \\

    \frac{\partial n}{\partial \alpha_1}|_{t=t_2} &

    \frac{\partial n}{\partial \alpha_2}|_{t=t_2} & \dots &

    \frac{\partial n}{\partial \alpha_n}|_{t=t_2} \\

    \vdots & \vdots & \vdots & \vdots \\

    \frac{\partial n}{\partial \alpha_1}|_{t=t_{66}} &

    \frac{\partial n}{\partial \alpha_2}|_{t=t_{66}} & \dots &

    \frac{\partial n}{\partial \alpha_n}|_{t=t_{66}} \\\end{pmatrix}
\\

  \bm{x} = \begin{pmatrix}

    \Delta \alpha_1 \\ \Delta \alpha_2 \\ \vdots \\ \Delta \alpha_n \\\end{pmatrix}

　　通过将残差向量\bm{d}的欧几里德范数（L2范数）||\bm{d}||^2最小化，可以求出参数的变化量\bm{x}为

\bm{x} = (A^T A)^{-1} A^T \bm{y}

　　通过\alpha_i \rightarrow \alpha_i + \Delta \alpha_i对参数进行更新，重复上述步骤10次，便可大致求出精度较好的回归曲线系数。