基础回归问题及python实现

1.简单线性回归

简单线性问题要解决的问题就是已知数据$x=(x_1,x_2,\cdots .x_n)$,$y=(y_1,y_2,\cdots ,y_n)$具有线性关系（这一般从散点图上可以看出来），设$y=ax+b$，我们要找出最佳的a、b值。而所谓最佳，我们是取使以下“代价函数”（损失函数）最小的a、b值。$$代价函数w=\sum _{i=1}^n(a{x}_i+b-y_i{)}^2$$

1.1最小二乘法

最小二乘法是解决简单线性回归问题最常用的方法，也是实现起来最简单的方法。
注意到$w$是关于a、b的二元函数，有$$\frac{\partial w}{\partial a}=2\sum _{i=1}^n(a{x}_i+b-y_i)x_i=2(a\sum _{i=1}^n{x}_i^2+b\sum _{i=1}^n{x}_i-\sum _{i=1}^n{x}_i{y}_i)=0\cdots (1)$$$$\frac{\partial w}{\partial b}=2\sum _{i=1}^n(a{x}_i+b-y_i)=2(a\sum _{i=1}^n{x}_i+nb-\sum _{i=1}^n{y}_i)=0\cdots (2)$$由(2)式可以解出$$b=\frac{{\sum }_{i=1}^n{y}_i-a{\sum }_{i=1}^n{x}_i}{n}=\bar{y}-a\bar{x}$$将b代入（1）式并化简可解出$$a=\frac{{\sum }_{i=1}^n(x_i-\bar{x})(y_i-\bar{y})}{{\sum }_{i=1}^n(x_i-\bar{x}{)}^2}$$
有了以上结果，我们可以直接用程序实现a、b的计算，从而拟合出直线。
以下为一实例

import numpy as np
import matplotlib.pyplot as plt
#初始数据集
x = np.arange(20)
y = np.array([0.4, 0.8, 1.1, 2.1, 2.8, 2.7, 3.5, 4.6, 5.1, 4.5, 6.0, 5.5, 6.9, 6.8, 7.6, 8.0, 8.8, 8.5, 9.5, 9.3])
#计算x、y平均值
x_mean = np.mean(x)
y_mean = np.mean(y)
#计算a、b
m1 = np.dot(x-x_mean,y-y_mean)
m2 = np.dot(x-x_mean,x-x_mean)
a = m1/m2
b = y_mean-a*x_mean
#计算代价函数
w = np.dot(a*x+b-y,a*x+b-y)
print('a=',a,'\nb=',b)
print(w)
#作图显示
Y = a*x+b
plt.scatter(x,y,c='r')
plt.plot(x,Y,c='b')
plt.show()

结果及图像如下：
a= 0.4915037593984963
b= 0.5557142857142852
2.849496240601506

1.2梯度下降法

有关梯度下降法可以参考笔者另一文章：
梯度下降法

2.简单非线性回归

在上面的问题中，数据集y与x之间呈线性关系，但生活中还有很多情况并非如此，比如$y=a{x}^2+bx+c、y=a{e}^{bx}$等情况，那么我们也可以根据上述思想去作拟合。这里我们假设通过散点图假设的预测函数为$h(x)$，此时的代价函数为$$w=\sum _{i=1}^n(h(x_i)-y_i{)}^2$$现在我们要解决的问题便是如何选择$h(x)$中的参数
但随着预测函数的变化，我们会发现代价函数对于参数的偏导变得很复杂，利用最小二乘法的方法得到的方程组解起来并不容易甚至不可解，因此最小二乘法有一定局限性。对于此类函数应用梯度下降法可以解决。

3.多元线性回归

在上面的讨论中，因变量y与自变量x都是一维，在实际生活中，我们还可能遇到自变量有多个指标的情况，比如一名学生的总成绩与各个学科成绩均相关。此时y与x的关系式为：$$y={\theta }_1{x}_1+{\theta }_2{x}_2+\cdots +{\theta }_n{x}_n+{\theta }_0\cdots (\ast )$$
它表示自变量有n个特征。对于以上表达式，我们也可以应用梯度下降法进行求解，也可以直接使用正规方程解。下面对于正规方程解作一简单介绍。
设$$\theta =\left(\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right),X=\left(\begin{array}{ccccc}1& X_{11}& X_{21}& \cdots & X_{n1}\\ 1& X_{12}& X_{22}& \cdots & X_{n2}\\ ⋮& ⋮& ⋮& & ⋮\\ 1& X_{1m}& X_{2m}& \cdots & X_{nm}\end{array}\right),y=\left(\begin{array}{c}{y}_1\\ y_2\\ ⋮\\ y_m\end{array}\right)$$则$(\ast )$式可以写成$$y=X\theta$$
我们的结论是$$\theta =(X^{T}X{)}^{-1}{X}^{T}y$$此时代价函数取到最小

4.利用sklearn库解决线性回归问题

python中的机器学习库sklearn对多种算法进行了封装，方便我们直接使用。

4.1简单线性回归的例子

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression as LR
#产生测试样本集：
x = np.random.uniform(-10,10,100)
X = x.reshape(-1,1)
y = 3*x+2+np.random.normal(0,1,100)
#进行线性回归：
lr = LR()
lr.fit(X,y)
y_pr = lr.predict(X)
print(lr.coef_,lr.intercept_)#打印系数、截距
plt.scatter(x,y,c='r')#作散点图
plt.plot(np.sort(x),y_pr[np.argsort(x)],c='b')#作预测直线图
plt.show()

运行结果：
[2.998317] 1.9778929070034075

4.2多元线性回归的例子

以波士顿房价的例子为例：

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
from sklearn.datasets import load_boston
#加载样本：
boston = load_boston()
#定义样本数据集：
x = boston['data']
y = boston['target']
#进行回归分析：
lr = LinearRegression()
lr.fit(x,y)
print(lr.coef_,lr.intercept_)#输出各特征系数与截距
#将预测值与真实值作图比较：
y_pr = lr.predict(x)
plt.scatter(y,y_pr)
X = np.arange(10,50,0.01)
Y = X
plt.plot(X,Y,c='r')
plt.title('y_predict-y_true')
plt.show()

结果：
[-1.08011358e-01 4.64204584e-02 2.05586264e-02 2.68673382e+00
-1.77666112e+01 3.80986521e+00 6.92224640e-04 -1.47556685e+00
3.06049479e-01 -1.23345939e-02 -9.52747232e-01 9.31168327e-03
-5.24758378e-01] 36.459488385089855

可以看出除去y=50左右的一些样本外，其余样本基本分布在y=x上下。

5.多项式回归

在一元非线性回归中，有一类比较特殊的关系，即因变量y与自变量x满足：$$y=\sum _{i=0}^n{a}_i{x}^i$$
解决此类多项式回归问题，我们可以用“增维”的思想：即将自变量由一维变换到多维。我们把x的每个幂次看作是一个特征，这样一个n次多项式便产生了n个特征，接下来将y于这n个特征作多元线性回归分析即可，我们通过下面的例子来进一步学习。

import numpy as np
import matplotlib.pyplot as plt
from sklearn.linear_model import LinearRegression
#生成数据集：
x = np.random.uniform(-2,2,100)
X = x.reshape(-1,1)
y = 2*x**3-3.5*x**2+1.7*x+12+np.random.normal(0,1,100)
#增加特征数量：
X1 = X**2
X2 = X**3
X3 = np.hstack([X2,X1,X])
#多元线性回归分析：
lr = LinearRegression()
lr.fit(X3,y)
print(lr.coef_,lr.intercept_)
#作图比较：
y_pr = lr.predict(X3)
plt.scatter(x,y,c='b')
plt.plot(np.sort(x),y_pr[np.argsort(x)],c='r')
plt.show()

运行结果：
[ 1.89709582 -3.57319401 1.84037251] 12.102429122160789

可以看出运行出的系数与我们产生样本集的系数基本相同，产生偏差是因为我们加上了误差导致的，预测曲线基本符合散点图也能看出这一点。
以上程序是我们通过增维方法来写的，便于我们理解问题的解决过程，在应用中我们可以用另一种方法来解决多项式回归问题。

import numpy as np
import matplotlib.pyplot as plt

#生成数据集：
x = np.random.uniform(-2,2,100)
y = 2*x**3-3.5*x**2+1.7*x+12+np.random.normal(0,1,100)
#回归分析：
f = np.polyfit(x,y,3)#前两个参数为两个样本集，第三个参数为欲用几次多项式进行回归分析
print(f)#输出系数
#作图比较：
p = np.poly1d(f)
y_pr = p(x)
plt.scatter(x,y,c='b')
plt.plot(np.sort(x),y_pr[np.argsort(x)],c='r')
plt.show()

结果：
[ 1.90769024 -3.38104398 1.88981261 11.91198708]

可以看出与上面一种方法基本相同，均比较成功地进行了回归分析。