线性回归分析：找出实验数据近似的线性函数关系

【摘要：课程通过一个具体案例，借助于EXECL分析工具，讲述线性回归分析必要的步骤和方法，并编写Python程序求出实验数据的拟合函数。】

线性回归分析

在数据分析中，线性回归分析是一种预测性的建模技术，它可以从给出的一组数据发现变量间的因果关系，建立变量间近似的线性函数模型。回归分析按照涉及变量的多少，分为一元回归和多元回归分析。
许多工程问题，常常需要根据两个变量的几组实验数据，来找出这两个变量的近似线性函数关系，用于预测和控制两个变量间的相互变化。

某厂1~4月份用水量用水量分析

例2 下表是某厂1~4月份用水量（单位：百吨）的一组数据：

根据上面的一组数据建立月份和用水量之间的近似线性函数，若设月份为x，用水量为y，近似的线性函数为y = f(x)，也就是，要找出一个能使上述数据大体适合的函数关系y = f(x)。
下面给出具体的回归分析步骤。
首先，要确定f(x)的函数类型，看f(x)是否近似于线性函数，当f(x)近似于线性函数时，才适合做线性回归分析。
我们可以使用EXECL来确定f(x)函数类型，新建EXECL文档，输入上面的数据，月份占一列，用水量占一列。

然后选择所有数据，让EXECL插入散点图。

从用水量数据的散点图可以看出，这些点的连线大致接近于一条直线，可以认为f(x)是线性函数。
设：f(x) = ax + b
其中a和b是待定系数，分析的主要工作就是依据给出的数据确定a和b的值，从而确定f(x)函数。

上图的红色趋势线就是我们要找的f(x)函数图形。从图中可以看出，f(x)没有经过图中所有的点，因为这些点本来就不在同一条直线上，函数f(x)反映了用水量随着月份变化的趋势。
如何确定a和b呢？可以按照下面的思路来确定a和b。
让f(x)在x1，x2，x3，x4处的函数值与表中数据y1，y2，y3，y4相差都很小，就是要使偏差：
y(i) – f(x(i)) (i=1,2,3,4)
都很小，如下图所示：

实际上，我们要求的是f(x)在x1，x2，x3，x4处的函数值与表中数据y1，y2，y3，y4的总偏差最小，即：
| (y1+y2+y3+y4) -(f(x1)+f(x2)+f(x3)+f(x4)) | = β
若β取得最小值，则说明f(x)就是我们要求的红色趋势线。为了计算方便，对每个偏差求平方后再相加，这样就可以去掉绝对值了：
(y1^2+y2^2+y3^2+y4^2) - (f(x1)^2+f(x2)^2+f(x3)^2+f(x4)^2) = β
求总偏差公式可进一步描述为:

其中，y(i)为表中所列的经验数据（用水量），x(i)为月份。
这种根据偏差的平方和作为最小的条件来选择系数a、b的方法叫做最小二乘法，是线性回归分析经常采用的方法。
下面的问题是如何确定a和b的系数，可以使M取得最小值。将经验数据和月份代入求总偏差公式：
(4.5-(a+b))^2 + (4-(2a+b))^2 + (3-(3a+b))^2 + (2.5-(4a+b))^2 = M
在上面的公式中，我们希望使所有偏差的平方和最小，如何求最小值M呢？可以通过微积分的方法得到，把偏差的平方和看作函数，它有a和b两个变量，求这个函数的最小值。
该函数是二元二次函数，分别求变量a和b的偏导函数，令偏导数为0，M取得最小值。求导后得到下面的方程组：
60.0*a + 20.0*b - 63.0 = 0
20.0*a + 8.0*b - 28.0 = 0

程序编写

from sympy import diff
from sympy import symbols
# 定义计算偏导的函数
def func(a,b):
return (4.5-(a+b)) ** 2 + \
(4-(2*a+b)) ** 2 + \
(3-(3*a+b)) ** 2 + \
(2.5-(4*a+b)) ** 2
a = symbols("a")
b = symbols("b")
# 计算变量a的偏导函数
print(diff(func(a,b),a))
# 计算变量b的偏导函数
print(diff(func(a,b),b))

程序需要导入sympy库，SymPy库是一个计算机代数系统，它支持符号计算、高精度计算、模式匹配、绘图、解方程、微积分、组合数学、离散数学、几何学、概率与统计、物理学等方面的功能。
求解上面的方程组，即可求得a和b的值，确定f(x)函数。求解方程组可以使用消元法，也可以利用Python的NumPy库来求解方程组。
NumPy是Python中科学计算的基础软件包，它提供了众多数学运算工具，这些数学运算工具包括：线性代数中的矩阵和向量运算、傅里叶变换、多维数组运算、数据统计运算以及丰富的数学函数库。

# 导入numpy库
import numpy as np
# 程序入口
if __name__ == '__main__':

# 建立2X2矩阵
ta = np.array([[60.0,20.0],[20.0,8.0]])
# 建立2维向量
tb = np.array([63.0,28.0])
# 解线性方程组
x = np.linalg.solve(ta,tb)
# 输出x和y的值
print("a = %.2f b = %.2f" % (x[0],x[1]))
运行上面的Python程序，求得:
a = -0.70 , b = 5.25
f(x)函数为：
f(x) = -0.7x + 5.25

验证f(x)函数，可以得到下面的数据表：