最小二乘回归策略

本文总结了最小二乘法针对一元线性方程、多元线性方程与多元非线性方程的回归策略，内容笔记向。

一元线性回归

a = 0，单变量x，单个未知数 b，求解一元一次方程

模型： $y=ax+b$

目标函数：$y=c$

n个样本点：$(x_i,y_i),i=1,2,3,...,n$

偏差平方和：$e=\sum _{i=1}^n(c-y_i{)}^2$

求解过程：

​ 偏差平方和函数： $$e=\sum _{i=1}^n(c-y_i{)}^2$$ ​ 对c求导，令导数为0： $$\frac{de}{dc}=\sum _{i=1}^{n}2(c-y_i)=0$$ ​ 解关于c的一元一次方程可得： $$c=\frac{\begin{array}{}\text{(1)}& \sum _{i=1}^n{y}_i\end{array}}{n}=\bar{y}$$ ​ 由此可知，拟合一条与x轴平行的直线的最佳参数，就是样本点y值的均值

a $\ne$ 0，单变量x，两个未知数 a和b，解二元一次方程

模型： $y=ax+b$

目标函数：$y=ax+b$

n个样本点：$(x_i,y_i),i=1,2,3,...,n$

偏差平方和：$e=\sum _{i=1}^n((a{x}_i+b))-y_i{)}^2$

求解过程：

偏差平方和函数： $$e=\sum _{i=1}^n((a{x}_i+b))-y_i{)}^2$$ 对a、b的偏导数： $$\begin{array}{}\text{(2)}& & \frac{\partial e}{\partial a}=\sum _{i=1}^{n}2x_i(a{x}_i+b-y_i)\text{(3)}& \text{(4)}& & \frac{\partial e}{\partial b}=\sum _{i=1}^{n}2(a{x}_i+b-y_i)\text{(5)}& \end{array}$$

同时令e对a和b的偏导等于0： $$\{\begin{array}{l}\begin{array}{}\text{(6)}& & \sum _{i=1}^n{x}_i(a{x}_i+b-y_i)=0\text{(7)}& \text{(8)}& & \sum _{i=1}^n(a{x}_i+b-y_i)=0\end{array}\end{array}$$ 解关于a、b的二元一次方程组得： $$\{\begin{array}{l}\begin{array}{}\text{(9)}& & a=\frac{n\sum x_i{y}_i-\sum x_i\sum y_i}{n\sum x_i^2-(\sum x_i{)}^2}\text{(10)}& \text{(11)}& \text{(12)}& & b=\frac{\sum y_i-a\sum x_i}{n}\end{array}\end{array}$$

多元线性回归

多项式函数：n个变量x，(n+1)个未知数a，解多元线性方程组

模型：$y=a_0+a_1{x}_1+a_2{x}_2+...+a_n{x}_n$

目标函数： $y=a_0+a_1{x}_1+a_2{x}_2+...+a_n{x}_n$

m组样本点：$(x_{1i},x_{2i},\dots ,x_{ni},y_i),i=1,2,...,m$

偏差平方和：$e=\sum _{i=1}^m(a_0+a_1{x}_{1i}+a_2{x}_{2i}+...+a_n{x}_{ni}-y_i{)}^2$

求解过程：

偏差平方和函数： $$\begin{array}{}\text{(13)}& e=\sum _{i=1}^m(a_0+a_1{x}_{1i}+a_2{x}_{2i}+...+a_n{x}_{ni}-y_i{)}^2\end{array}$$ 令e对 $a_0$ ~ $a_n$ 的各个偏导为0，稍微变化形式可得到关于未知数的正规方程组（n+1个方程，n+1个未知数）： $$\{\begin{array}{l}\begin{array}{}\text{(14)}& \sum _{i=1}^m(a_0+a_1{x}_{1i}+a_2& x_{2i}+...+a_n{x}_{ni})=\sum _{i=1}^m{y}_i\text{(15)}& \text{(16)}& \sum _{i=1}^m(a_0+a_1{x}_{1i}+a_2& x_{2i}+...+a_n{x}_{ni})x_{1i}=\sum _{i=1}^m{x}_{1i}{y}_i\text{(17)}& \text{(18)}& & ⋮\text{(19)}& \text{(20)}& \sum _{i=1}^m(a_0+a_1{x}_{1i}+a_2& x_{2i}+...+a_n{x}_{ni})x_{ni}=\sum _{i=1}^m{x}_{ni}{y}_i\end{array}\end{array}$$ 用矩阵的形式表示为： $$\begin{array}{r}\left[\begin{array}{ccccc}n& \sum x_{1i}& \sum x_{2i}& \dots & \sum x_{ni}\\ \\ \sum x_{1i}& \sum x_{1i}^2& \sum x_{1i}{x}_{2i}& \dots & \sum x_{1i}{x}_{ni}\\ \\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \\ \sum x_{ni}& \sum x_{1i}{x}_{ni}& \sum x_{2i}{x}_{ni}& \sum \dots & \sum x_{ni}^2\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_n\end{array}\right]=\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ \\ x_{11}& x_{12}& x_{13}& \dots & x_{1m}\\ \\ ⋮& ⋮& ⋮& \ddots & ⋮\\ \\ x_{n1}& x_{n2}& x_{n3}& \dots & x_{nm}\end{array}\right]\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\\ y_m\end{array}\right]\end{array}$$ 令样本矩阵以及参数矩阵如下： $$\begin{array}{}\text{(21)}& X=\left[\begin{array}{cccccc}& 1& x_{11}& x_{21}& \dots & x_{n1}\\ \\ & 1& x_{12}& x_{22}& \dots & x_{n2}\\ \\ & ⋮& ⋮& ⋮& \ddots & ⋮\\ \\ & 1& x_{1m}& x_{2m}& \dots & x_{nm}\end{array}\right]，Y=\left[\begin{array}{c}{y}_1\\ y_2\\ y_3\\ ⋮\\ y_m\end{array}\right]，A=\left[\begin{array}{c}{a}_0\\ a_1\\ a_2\\ ⋮\\ a_n\end{array}\right]\end{array}$$ 则上面的矩阵方程可以表示为： $$\begin{array}{}\text{(22)}& (X^{T}X)A=XY\end{array}$$ 由于 $X^{T}X$ 为满秩矩阵 ，即非奇异、可逆，因此将上式等号两侧同时左乘 $(X^{T}X{)}^{-1}$，可得最优参数矩阵： $$\begin{array}{}\text{(23)}& A=(X^{T}X{)}^{-1}XY\end{array}$$

多元非线性回归

一元非线性回归可通过换元法转化为一元线性回归

四参数正弦函数：四个未知数，解多元非线性方程组

解法：视为整体，在初始值处一阶泰勒展开，结合牛顿-拉弗森迭代法每次构造一个线性方程组，用高斯消元法求解，再代回去迭代直至符合误差预期

模型：$y=Asin(wt+\theta )+c$

目标函数：$y=Asin(wt+\theta )+c$

样本点：$(t_i,y_i),i=1,2,3...n$

偏差平方和：$e=\sum _{i=1}^n[(Asin(w{t}_i+\theta )+c)-y_i{]}^2$

求解过程：

偏差平方和函数： $$\begin{array}{}\text{(24)}& & e=\sum _{i=1}^n[(Asin(w{t}_i+\theta )+c)-y_i{]}^2\text{(25)}& \end{array}$$ 对各个待求参数求偏导： $$\begin{array}{}\text{(26)}& & \frac{\partial e}{\partial A}=\sum _{i=1}^{n}2[Asin(w{t}_i+\theta )+c-y_i]sin(w{t}_i+\theta )\text{(27)}& \text{(28)}& & \frac{\partial e}{\partial w}=\sum _{i=1}^{n}2[Asin(w{t}_i+\theta )+c-y_i]A{t}_{i}cos(w{t}_i+\theta )\text{(29)}& \text{(30)}& & \frac{\partial e}{\partial \theta }=\sum _{i=1}^{n}2[Asin(w{t}_i+\theta )+c-y_i]Acos(w{t}_i+\theta )\text{(31)}& \text{(32)}& & \frac{\partial e}{\partial c}=\sum _{i=1}^{n}2[Asin(w{t}_i+\theta )+c-y_i]\end{array}$$ 令各个偏导数为0，等价于求解非线性方程组（注意和上式有差异，约掉了常系数）： $$\{\begin{array}{l}\begin{array}{}\text{(33)}& & f_1(X)=\sum _{i=1}^n[Asin(w{t}_i+\theta )+c-y_i]sin(w{t}_i+\theta )=0\text{(34)}& \text{(35)}& & f_2(X)=\sum _{i=1}^n[Asin(w{t}_i+\theta )+c-y_i]t_{i}cos(w{t}_i+\theta )=0\text{(36)}& \text{(37)}& & f_3(X)=\sum _{i=1}^n[Asin(w{t}_i+\theta )+c-y_i]cos(w{t}_i+\theta )=0\text{(38)}& \text{(39)}& & f_4(X)=\sum _{i=1}^n[Asin(w{t}_i+\theta )+c-y_i]=0\end{array}\end{array}$$ 这里采用牛顿-拉弗森迭代法的思想，线性化迭代求解该方程组的4个未知量 $A、w、\theta ,、c$

用矩阵形式表示参数集以及偏导数集，并假设初值为 $X_0$ $$X=\left[\begin{array}{c}A\\ w\\ \theta \\ c\end{array}\right]，F(X)=\left[\begin{array}{c}{f}_1(X)\\ f_2(X)\\ f_3(X)\\ f_4(X)\end{array}\right]$$ 构建 $F(X)$ 在 $X_0$ 处的雅克比矩阵： $$\begin{array}{}\text{(40)}& J_F(X)=\left[\begin{array}{cccc}\frac{\partial f_1(X_0)}{\partial A}& \frac{\partial f_1(X_0)}{\partial w}& \frac{\partial f_1(X_0)}{\partial \theta }& \frac{\partial f_1(X_0)}{\partial c}\\ \\ \frac{\partial f_2(X_0)}{\partial A}& \frac{\partial f_2(X_0)}{\partial w}& \frac{\partial f_2(X_0)}{\partial \theta }& \frac{\partial f_2(X_0)}{\partial c}\\ \\ \frac{\partial f_3(X_0)}{\partial A}& \frac{\partial f_3(X_0)}{\partial w}& \frac{\partial f_3(X_0)}{\partial \theta }& \frac{\partial f_3(X_0)}{\partial c}\\ \\ \frac{\partial f_4(X_0)}{\partial A}& \frac{\partial f_4(X_0)}{\partial w}& \frac{\partial f_4(X_0)}{\partial \theta }& \frac{\partial f_4(X_0)}{\partial c}\end{array}\right]\end{array}$$

其中f1对各个参数的偏导： $$\begin{array}{}\text{(41)}& & \frac{\partial f_1}{\partial A}=\sum _{i=1}^{n}si{n}^2(w{t}_i+\theta )\text{(42)}& \text{(43)}& & \frac{\partial f_1}{\partial w}=\sum _{i=1}^n[A{t}_{i}sin2(w{t}_i+\theta )]+\sum _{i=1}^n[(c-y_i)t_{i}cos(w{t}_i+\theta )]\text{(44)}& \text{(45)}& & \frac{\partial f_1}{\partial \theta }=\sum _{i=1}^n[Asin2(w{t}_i+\theta )]+\sum _{i=1}^n[(c-y_i)cos(w{t}_i+\theta )]\text{(46)}& \text{(47)}& & \frac{\partial f_1}{\partial c}=\sum _{i=1}^{n}sin(w{t}_i+\theta )\end{array}$$

f2对各个参数的偏导： $$\begin{array}{}\text{(48)}& & \frac{\partial f_2}{\partial A}=\frac{1}{2}\sum _{i=1}^n{t}_{i}sin2(w{t}_i+\theta )\text{(49)}& \text{(50)}& & \frac{\partial f_2}{\partial w}=\sum _{i=1}^n[A{t}_i^{2}cos2(w{t}_i+\theta )]-\sum _{i=1}^n[(c-y_i)t_i^{2}sin(w{t}_i+\theta )]\text{(51)}& \text{(52)}& & \frac{\partial f_2}{\partial \theta }=\sum _{i=1}^n[A{t}_{i}cos2(w{t}_i+\theta )]-\sum _{i=1}^n[(c-y_i)t_{i}sin(w{t}_i+\theta )]\text{(53)}& \text{(54)}& & \frac{\partial f_2}{\partial c}=\sum _{i=1}^n{t}_{i}cos(w{t}_i+\theta )\end{array}$$

f3对各个参数的偏导： $$\begin{array}{}\text{(55)}& & \frac{\partial f_3}{\partial A}=\frac{1}{2}\sum _{i=1}^{n}sin2(w{t}_i+\theta )\text{(56)}& \text{(57)}& & \frac{\partial f_3}{\partial w}=\sum _{i=1}^n[A{t}_{i}cos2(w{t}_i+\theta )]-\sum _{i=1}^n[(c-y_i)t_{i}sin(w{t}_i+\theta )]\text{(58)}& \text{(59)}& & \frac{\partial f_3}{\partial \theta }=\sum _{i=1}^n[Acos2(w{t}_i+\theta )]-\sum _{i=1}^n[(c-y_i)sin(w{t}_i+\theta )]\text{(60)}& \text{(61)}& & \frac{\partial f_3}{\partial c}=\sum _{i=1}^{n}cos(w{t}_i+\theta )\end{array}$$

f4对各个参数的偏导： $$\begin{array}{}\text{(62)}& & \frac{\partial f_4}{\partial A}=\sum _{i=1}^{n}sin(w{t}_i+\theta )\text{(63)}& \text{(64)}& & \frac{\partial f_4}{\partial w}=\sum _{i=1}^n[A{t}_{i}cos(w{t}_i+\theta )]\text{(65)}& \text{(66)}& & \frac{\partial f_4}{\partial \theta }=\sum _{i=1}^n[Acos(w{t}_i+\theta )]\text{(67)}& \text{(68)}& & \frac{\partial f_4}{\partial c}=n\end{array}$$

$F(X)$ 在初值 $X_0$ 处的一阶泰勒展开： $$\begin{array}{}\text{(69)}& F(X)\approx F(X_0)+J_F(X_0)(X-X_0)\end{array}$$ 可知：令 等式右侧为0 求解出的 $X$ 会比 $X_0$ 更加趋近于 $F(X)=0$ 的真实解，只要初值在这个真实解的收敛域内，经过不断的迭代求得的近似解会越来越趋近于真实解，我们令 $\mathrm{\Delta }X=X-X_0$ 可得： $$\begin{array}{}\text{(70)}& & F(X_0)+J_F(X_0)\mathrm{\Delta }X=0\end{array}$$ 转化一下形式： $$\begin{array}{}\text{(71)}& J_F(X_0)\mathrm{\Delta }X=-F(X_0)\end{array}$$ 这样我们便将非线性方程的求解问题转化成了线性方程的求解问题，上述式子就等价于： $$\begin{array}{}\text{(72)}& \text{(73)}& \left[\begin{array}{cccc}\frac{\partial f_1(X_0)}{\partial A}& \frac{\partial f_1(X_0)}{\partial w}& \frac{\partial f_1(X_0)}{\partial \theta }& \frac{\partial f_1(X_0)}{\partial c}\\ \frac{\partial f_2(X_0)}{\partial A}& \frac{\partial f_2(X_0)}{\partial w}& \frac{\partial f_2(X_0)}{\partial \theta }& \frac{\partial f_2(X_0)}{\partial c}\\ \frac{\partial f_3(X_0)}{\partial A}& \frac{\partial f_3(X_0)}{\partial w}& \frac{\partial f_3(X_0)}{\partial \theta }& \frac{\partial f_3(X_0)}{\partial c}\\ \frac{\partial f_4(X_0)}{\partial A}& \frac{\partial f_4(X_0)}{\partial w}& \frac{\partial f_4(X_0)}{\partial \theta }& \frac{\partial f_4(X_0)}{\partial c}\end{array}\right]\left[\begin{array}{c}\mathrm{\Delta }A\\ \mathrm{\Delta }w\\ \mathrm{\Delta }\theta \\ \mathrm{\Delta }c\end{array}\right]=\left[\begin{array}{c}-f_1(X_0)\\ -f_2(X_0)\\ -f_3(X_0)\\ -f_4(X_0)\end{array}\right]\text{(74)}& \end{array}$$ 用高斯消元法或奇异值分解法解这个线性方程，可以求出 $\mathrm{\Delta }A，\mathrm{\Delta }w,\mathrm{\Delta }\theta ,\mathrm{\Delta }c$，进而本次迭代的结果就是 $(X_0+\mathrm{\Delta }X)$，然后我们将$(X_0+\mathrm{\Delta }X)$ 作为新的 $X_0$ ，代入上面的线性方程，不断迭代求解，$\mathrm{\Delta }X$ 就会越来越小，可以获得越来越趋近于真实解的 $X$。

Python API -- Leastsq功能测试

需要预先安装scipy库： sudo pip install scipy

这个程序是我在学最小二乘的时候写的一个测试程序，测试的目标是scipy库里leastsq这个以最小二乘法为基本原理的回归函数，目标函数可以是线性的也可以是非线性的。为了两个都测试，我选择了常量函数y=C作为线性模型，四参数正弦函数y=Asin(wt+Θ)+C作为非线性模型。其余各个函数的意义如下：

leastsq：API，输入参数依次为：误差函数、拟合初值、样本点集、迭代次数上限

func_sin：求正弦函数值。输入自变量x的值与四个参数组成的列表p，返回由p组成的正弦函数在x处的值。

error_sin：求拟合值与样本值的偏差量。输入拟合出来的参数列表p0与样本因变量y，返回二者的偏差量。

fit_sin：调用Leastsq函数进行拟合，获取拟合结果并求积分，返回拟合结果或积分值。

func_c/error_c/fit_c：模型换成y=C，与上面几个对照

read：纯粹是为了个人debug方便写的一个读取语句和捕获输入异常的功能函数，建议根据个人需求改写main函数。

import numpy as np
from scipy.optimize import leastsq


def func_sin(x, p):
    a, f, theta, c = p
    return a * np.sin(2 * np.pi * f * x + theta) + c


def error_sin(p0, x, y):
    return y - func_sin(x, p0)


def fit_sin(x, y, p, t, flag, max_fev=500):
    x_sample = np.array(x)
    y_sample = np.array(y)
    x_max = max(x_sample)

    p0 = np.array(p)
    para = leastsq(error_sin, p0, args=(x_sample, y_sample), maxfev=max_fev)

    a1, a2, a3, a4 = para[0]
    a1 = round(float(a1), 3)
    a2 = round(float(a2), 3)
    a3 = round(float(a3), 3)
    a4 = round(float(a4), 3)

    # result_api = round(integrate(expr, (x_t, x_max, x_max+t)), 3)  # 积分api,速度奇慢
    tmp = a1 / (2*np.pi*a2)
    temp_x = x_max + t
    lower = - tmp - tmp*np.cos(2*np.pi*a2*x_max + a3) + a4*x_max
    upper = - tmp - tmp*np.cos(2*np.pi*a2*temp_x + a3) + a4*temp_x
    result = round((upper-lower), 3)
    # print("\nresult_api = {}, result_straight = {}".format(result_api, result))

    if flag:
        return [a1, a2, a3, a4], result
    else:
        return result


def func_c(x, c):
    return c * np.ones(x.shape)


def error_c(p, x, y):
    return y - func_c(x, p)


def fit_c(x, y, p, t, flag, max_fev=500):
    x_sample = np.array(x)
    y_sample =np.array(y)
    para = leastsq(error_c, p, args=(x_sample, y_sample), maxfev=max_fev)
    a = round(float(para[0]), 3)
    result = round(t*a, 3)
    if flag:
        return a, result
    else:
        print("拟合结果: y_c =", a)
        return result


# (独立测试使用) 按规模读取键盘数字输入, 捕获异常
def read(sentence, quantity):
    assert quantity > 0
    sentence_error = "格式错误, 请按照要求重新输入: "
    inf = []
    try:
        inf = list(map(float, input(sentence).split()))
        while len(inf) != quantity:
            inf = list(map(float, input(sentence_error).split()))
    except ValueError:
        while True:
            try:
                inf = list(map(float, input(sentence_error).split()))
                while len(inf) != quantity:
                    inf = list(map(float, input(sentence_error).split()))
                break
            except ValueError:
                continue

    if quantity == 1:
        return inf[0]
    else:
        return inf


# 作为独立程序测试
if __name__ == "__main__":

    # 样本点要求: 跨度小 或者 跨度大密度大, 这里的跨度指的是样本横坐标跨越的周期个数

    from sympy import *
    from matplotlib import pyplot as plt
    import time

    sentence_1 = "请用空格分割, 按顺序输入希望的测试点 数目 & 横坐标最小值 & 横坐标最大值 & 横坐标分度值(如:15 -3.14 0 0.1): "
    sentence_2 = "输入w测试正弦模型, e测试直线模型, 其它键取消测试: "
    sentence_3 = "请用空格分割, 请按顺序输入正弦模型y = Asin(2πf*t+Θ)+c的四个参数A,f,Θ,c的真实值: "
    sentence_4 = "请用空格分割, 输入参数A,f,Θ,c的拟合初值: "
    sentence_5 = "请输入直线模型y = c中参数c的真实值: "
    sentence_6 = "请输入参数c的拟合初值: "
    sentence_7 = "请输入预测时长(单位:s): "
    sentence_8 = "分度值不可为0, 请重新输入: "
    sentence_9 = "请输入拟合次数上限: "

    x_samples, y_samples, x_fit, y_fit = [], [], [], []

    info = read(sentence_1, 4)
    if info[3] == 0:
        info = read(sentence_8, 4)
    amount, x_min, x_max, x_division = int(info[0]), info[1], info[2], info[3]

    key = input(sentence_2)

    if key == 'w':
        a, f, theta, c = read(sentence_3, 4)  # 真实值
        p0 = read(sentence_4, 4)  # 拟合初值
        x_samples = np.linspace(x_min, x_max, amount)  # 样本x
        y_origin = func_sin(x_samples, [a, f, theta, c])  # 样本y(无噪声)
        y_samples = y_origin + 0.8 * np.random.randn(amount)  # 样本y(带噪声)
        t = read(sentence_7, 1)  # 预测时长
        max_fev = int(read(sentence_9, 1))  # 拟合次数上限

        time_1 = time.time()
        paras, result = fit_sin(x_samples, y_samples, p0, t, true, max_fev)  # 拟合
        time_2 = time.time()

        x_t = symbols("x")
        expr = paras[0] * sin(round(2 * np.pi * paras[1], 3) * x_t + paras[2]) + paras[3]

        print("\n模型: y = Asin(2πf*t+Θ)+c")
        print("样本点:\nx={}\ny(无噪声)={}\ny(含噪声)={}".format(x_samples, y_origin, y_samples))
        print("函数原型参数: ", [a, f, theta, c])
        print("拟合初值: ", p0)
        print("拟合结果: ", paras)
        print("拟合表达式: y =", expr)
        print("拟合耗时: {} ms".format(round((time_2-time_1)*1000, 3)))
        print("积分区间: [{}, {}]".format(x_max, x_max+t))
        print("积分值: ", result)

        if x_max > 0:
            x_fit = np.arange(x_min, 2*(x_max+t), x_division)
            x_bound = np.array([x_min, 2*(x_max+t)])
        else:
            x_fit = np.arange(x_min, abs(-x_max+2*t), x_division)
            x_bound = np.array([x_min, abs(-x_max+2*t)])

        y_fit = func_sin(x_fit, paras)
        y_bound = np.array([c, c])

        plt.figure(figsize=(8, 6))
        plt.scatter(x_samples, y_samples, color='red', label='Sample Points', linewidth=2)  # 带噪声的样本点图
        plt.plot(x_samples, y_samples, color='gray', label="Sample Lines", linewidth=2)  # 带噪声的样本曲线图
        plt.plot(x_fit, y_fit, color='orange', label="Fitting Line", linewidth=2)  # 拟合的曲线图
        plt.plot(x_bound, y_bound, color='gray', linewidth=1, linestyle='--')  # 对称轴
        plt.xlabel("Time")
        plt.ylabel("Velocity")
        plt.plot(x_max, func_sin(x_max, paras), "ob", color='blue', label="current")
        plt.plot(x_max+t, func_sin(x_max+t, paras), "ob", color='purple', label='Prediction')

        plt.plot([x_max, x_max], [func_sin(x_max, paras), c],
                 color='blue', linestyle='--')  # 作最新点到对称轴的垂直虚线
        plt.plot([x_max+t, x_max+t], [func_sin(x_max+t, paras), c],
                 color='purple', linestyle='--')  # 作预测点到对称轴的垂直虚线
        x_series = np.linspace(x_max, x_max+t, amount)  # 积分区域
        y_1 = c
        y_2 = func_sin(x_series, paras)
        plt.fill_between(x_series, y_1, y_2, facecolor='white', hatch='/', interpolate=true)  # 填充阴影

        plt.legend()
        plt.show()
        exit(0)

    elif key == 'e':
        c = read(sentence_5, 1)
        p0 = read(sentence_6, 1)
        x_samples = np.linspace(x_min, x_max, amount)  # 样本x
        y_origin = func_c(x_samples, c)  # 样本y(无噪声)
        y_samples = y_origin + 0.8 * np.random.randn(amount)  # 样本y(带噪声)
        t = read(sentence_7, 1)  # 预测时长
        max_fev = int(read(sentence_9, 1))  # 拟合次数上限

        time_1 = time.time()
        para, result = fit_c(x_samples, y_samples, p0, t, true, max_fev)  # 拟合
        time_2 = time.time()

        print("\n模型: y = c")
        print("样本点:\nx={}\ny={}".format(x_samples, y_samples))
        print("函数原型参数: ", c)
        print("拟合初值: ", p0)
        print("拟合结果: ", para)
        print("拟合表达式: y =", para)
        print("拟合耗时: {} ms".format((time_2-time_1)*1000))
        print("积分区间: [{}, {}]".format(x_max, x_max+t))
        print("积分值: ", result)

        if x_max > 0:
            x_fit = np.arange(x_min, 2*(x_max+t), x_division)
            x_bound = np.array([x_min, 2.5*(x_max+t)])
        else:
            x_fit = np.arange(x_min, abs(-x_max+2*t), x_division)
            x_bound = np.array([x_min, 1.25 * abs(-x_max+2*t)])

        y_fit = func_c(x_fit, para)
        y_bound = np.array([c, c])

        fig = plt.figure(figsize=(8, 6))
        ax = fig.add_subplot(111)
        ax.spines['right'].set_color('none')  # 去掉右方的y轴
        ax.spines['top'].set_color('none')  # 去掉上方的x轴
        ax.spines['left'].set_position(('data', x_min))  # 将y轴移动至x最小值处
        ax.spines['bottom'].set_position(('data', min(y_samples)))  # 将x轴移动至y样本最小值处
        ax.set_xticks(np.arange(min(x_samples), x_max+2*t, 1))  # 设置x刻度
        ax.scatter(x_samples, y_samples, color='red', label='Original Points', linewidth=1)  # 有噪声的样本点图
        plt.plot(x_bound, y_bound, color='green', label='Real Value', linewidth=1)  # 实际
        ax.plot(x_fit, y_fit, color='orange', label="Fitting Line", linewidth=2)  # 有噪声样本拟合的曲线图
        plt.plot(x_max, func_c(np.array([x_max]), para), "ob", color='blue', label="current")  # 最新点
        plt.plot(x_max+t, func_c(np.array([x_max+t]), para), "ob", color='purple', label='Prediction')  # 预测点
        ax.plot([x_max, x_max], [func_c(np.array([x_max]), para), min(y_samples)],
                color='blue', linestyle='--')  # 作最新点到x轴的垂直虚线
        ax.plot([x_max+t, x_max+t], [func_c(np.array([x_max+t]), para), min(y_samples)],
                color='purple', linestyle='--')  # 作预测点到x轴的垂直虚线
        x_series = np.linspace(x_max, x_max+t, amount)  # 积分区域
        y_1 = min(y_samples) * np.ones(x_series.shape)
        y_2 = func_c(x_series, para)
        plt.fill_between(x_series, y_1, y_2, where=y_1<y_2, facecolor='white', hatch='/', interpolate=true)  # 填充阴影
        plt.xlabel("Time/s")
        plt.ylabel("Velocity")
        plt.legend()
        plt.show()
        exit(0)

    else:
        print("已取消测试")
        exit(0)

运行结果如下：

输入参数依次为：[15, -3.14, 0, 0.1], 'w', [8, 0.34, 0.52, 0], [7, 0.2, 0.2, 0], 0.7, 1000， 结果如下：

由于个人需要求了个积分，因此输入的参数中的"0.7"表示积分区间的长度，积分区间为 [样本横坐标最大值， 样本横坐标最大值+0.7]，由阴影部分表示。其余点/线的意义如下：

• 红色点：由输入参数的正弦函数生成的带高斯噪声的样本点
• 蓝色点：横坐标最大的样本点
• 紫色点：积分点
• 灰色折线：样本点依次连线
• 黄色曲线：拟合结果

输出的数据结果如下，其中带噪声的样本点是由无噪声样本点的基础上加上符合标准正态分布的噪声产生的:

程序可改进的地方：因为选择了y=C这样一个与自变量不相关的线性模型，所以最佳拟合值就是样本均值，不需要再调用api了！