织梦CMS - 轻松建站从此开始!

罗索

线性回归

jackyhwei 发布于 2020-04-02 10:30 点击:次 
一、回归问题的定义 回归是监督学习的一个重要问题,回归用于预测输入变量和输出变量之间的关系。回归模型是表示输入变量到输出变量之间映射的函数。回归问题的学习等价于函数拟合:
TAG: 线性回归  

一、回归问题的定义

回归是监督学习的一个重要问题,回归用于预测输入变量和输出变量之间的关系。回归模型是表示输入变量到输出变量之间映射的函数。回归问题的学习等价于函数拟合:使用一条函数曲线使其很好的拟合已知函数且很好的预测未知数据。

回归问题分为模型的学习和预测两个过程。基于给定的训练数据集构建一个模型,根据新的输入数据预测相应的输出。

回归问题按照输入变量的个数可以分为一元回归和多元回归;按照输入变量和输出变量之间关系的类型,可以分为线性回归和非线性回归。

 

一元回归:y = ax + b

多元回归:

二、回归问题的求解

2.1解析解

2.1.1 最小二乘法

最小二乘法又称最小平方法,它通过最小化误差的平方和寻找数据的最佳函数匹配。利用最小二乘法可以简便地求得未知的数据,并使得这些求得的数据与实际数据之间误差的平方和为最小。最小二乘法还可用于曲线拟合。

2.1.2利用极大似然估计解释最小二乘法

现在假设我们有m个样本,我们假设有:

误差项是IID,根据中心极限定理,由于误差项是好多好多相互独立的因素影响的综合影响,我们有理由假设其服从高斯分布,又由于可以自己适配theta0,是的误差项的高斯分布均值为0,所以我们有

所以我们有:

也即:

 表示在theta给定的时候,给我一个x,就给你一个y

那么我们可以写出似然函数:

由极大似然估计的定义,我们需要L(theta)最大,那么我们怎么才能是的这个值最大呢?两边取对数对这个表达式进行化简如下:

需要 l(theta)最大,也即最后一项的后半部分最小,也即:

所以,我们最后由极大似然估计推导得出,我们希望 J(theta) 最小,而这刚好就是最小二乘法做的工作。而回过头来我们发现,之所以最小二乘法有道理,是因为我们之前假设误差项服从高斯分布,假如我们假设它服从别的分布,那么最后的目标函数的形式也会相应变化。

好了,上边我们得到了有极大似然估计或者最小二乘法,我们的模型求解需要最小化目标函数J(theta),那么我们的theta到底怎么求解呢?有没有一个解析式可以表示theta?

2.1.3 theta的解析式的求解过程

我们需要最小化目标函数,关心 theta 取什么值的时候,目标函数取得最小值,而目标函数连续,那么 theta 一定为 目标函数的驻点,所以我们求导寻找驻点。

求导可得:

最终我们得到参数 theta 的解析式:

 关于向量、矩阵求导知识参见http://www.cnblogs.com/futurehau/p/6105236.html

上述最后一步有一些问题,假如 X'X不可逆呢?

我们知道 X'X 是一个办正定矩阵,所以若X'X不可逆或为了防止过拟合,我们增加lambda扰动,得到

从另一个角度来看,这相当与给我们的线性回归参数增加一个惩罚因子,这是必要的,我们数据是有干扰的,不正则的话有可能数据对于训练集拟合的特别好,但是对于新数据的预测误差很大。

2.1.4正则化

L2-norm: (Ridge回归)

L1-norm: (Lasso回归)

  J(theta) = J + lambda * sum(|theta|)

L1-norm 和 L2-norm都能防止过拟合,一般L2-norm的效果更好一些。L1-norm能够产生稀疏模型,能够帮助我们去除某些特征,因此可以用于特征选择。

L1-norm 和 L2-norm的直观理解:摘自http://lib.csdn.net/article/machinelearning/42049

 

  

今天又看到一个比较好的解释。可以把加入正则理解为加入约束条件,(类似于逆向拉格朗日)。那么,比如上边的图,L2约束就是一个圆,L1约束就是一个方形。那些关于w的圈圈都是等值线,代表了损失时多少,我们现在要求的就是在约束的条件下寻找最小的损失。所以其实就是找约束的图形和等值线的交点。

L1的缺点:如果有几个变量相关性比较大,那么它会随机的选择某一个。优化:Elastic Net

 2.2 梯度下降算法

我们在上边给出了最小二乘法求解线性回归的参数theta,实际python 的 numpy库就是使用的这种方法。

当然了,这需要我们的参数的维度不大,当维度大的时候,使用解析解就不适用了,这里讨论梯度下降算法。

2.2.1梯度下降法步骤:

  初始化theta

  沿着负梯度方向迭代,更新后的theta使得J(theta)更小。

  其中α表示学习率

  一个优化技巧:不同的特征采用不同的学习率 Adagrad

  

梯度下系那个示意图如下:

每次迭代找到一个更好的,最后得到一个局部最优解,不一定是全局最优,但是是堪用的。

2.2.2 具体实现

梯度方向:

2.2.2.1 批量梯度下降算法:

由于在线性回归中,目标函数收敛而且为凸函数,是有一个极值点,所以局部最小值就是全局最小值。

2.2.2.2随机梯度下降算法:

拿到一个样本就下降一次。实际中随机梯度下降算法其实是更适用的。出于一下亮点考虑:

1.由于噪声的存在,不能保证每个变量都让目标函数下降,但总的趋势是下降的。但是另一方面,由于噪声的存在,随机梯度下降算法往往有利于跳出局部最小值。

2.流式数据的处理

2.2.2.3 mini-batch

拿到若干个样本的平均梯度之后在更新梯度方向。

如上图所示,一个批量梯度下降,一个随机梯度下降,最终效果其实是相近的。

 2.2.2.4 上升一个高度把三种梯度下降算法联系起来

期望损失:理论上模型关于自变量因变量的平均意义下的损失,学习的目标就是选择期望损失最小的模型。

经验风险:模型关于训练样本集的平均损失。因为我们不可能得到所有的样本来计算期望损失,所以我们使用经验风险来代替期望损失。

那么怎么来处理选择这些样本呢?

BGD:我拥有的所有者n个样本的平均损失

SGD:单个样本处理

mini-batch:多个样本处理

 

三、实际线性回归时候的数据使用

此处分析不仅仅局限于线性回归。

实际中可以把数据分为训练数据和测试数据,然后根据不同模型在测试数据上的表现来选择模型。

另外一些情况,比如上边加上正则化之后,我们不能由训练数据得到lambda,那么我们需要把训练数据进一步划分为训练数据和验证数据。在训练数据上学习theta和lambda,然后在验证数据上选择lambda,然后再在测试数据上验证选择不同模型。

实际中采用交叉验证充分利用数据,例如五折交叉验证。

 

 四、几个系数定义说明

对于m个样本:    

某模型的估计值为: 

定义:

总平方和 TSS(Total Sum of Squares) :  即样本伪方差的m倍 Var(Y) = TSS / m

残差平方和 RSS(Residual Sum of Squares):   RSS也记作误差平方和SSE (Sum of Squares for Error)

可解释平方和ESS(Explained Sum of Squares) :  ESS又称为回归平方和SSR(Sum of Squares for Regression)

决定系数:

 

TSS >= RSS + ESS, 在无偏估计的时候取等号。

R^2越大,拟合效果越好。

 

需要额外说明的是,这里所谓的线性回归,主要是针对的参数theta,并不是针对x,我们可以对原来的数据进行处理,比如平方得到x^2的数据,然后把这个看作一个影响因素,这样最终得到的y关于x的图形就不是线性的,但当然这也是线性回归。

另外,还有局部加权的线性回归的概念,这部分内容在SVM中进一步解释。

 

五. Linear Regression for binary classfication

考虑到线性分类问题求解不方便,所以可不可以通过线性回归来求解线性分类问题呢? 

 两者的差别主要在于损失函数。

平方损失是0/1损失的一个上限。

所以,使用Linear Regression来求解Linear Regression也是有道理的。

 

 使用 sklearn 库来进行训练数据测试数据的划分学习 线性回归库的调用:

  1. #!/usr/bin/python 
  2. # -*- coding:utf-8 -*- 
  3.  
  4. import numpy as np 
  5. import matplotlib.pyplot as plt 
  6. import pandas as pd 
  7. from sklearn.model_selection import train_test_split 
  8. from sklearn.linear_model import Lasso, Ridge 
  9. from sklearn.model_selection import GridSearchCV 
  10.  
  11.  
  12. if __name__ == "__main__": 
  13.     # pandas读入 
  14.     data = pd.read_csv('8.Advertising.csv')    # TV、Radio、Newspaper、Sales 
  15.     x = data[['TV', 'Radio', 'Newspaper']] 
  16.     # x = data[['TV', 'Radio']] 
  17.     y = data['Sales'] 
  18.     print x 
  19.     print y 
  20.  
  21.     x_train, x_test, y_train, y_test = train_test_split(x, y, random_state=1
  22.     # print x_train, y_train 
  23.     model = Lasso() 
  24.     # model = Ridge() 
  25.  
  26.     alpha_can = np.logspace(-3, 2, 10)  #10^(-3) ~ 10^(2) 等比10个数 
  27.     lasso_model = GridSearchCV(model, param_grid={'alpha': alpha_can}, cv=5) #5折交叉验证 
  28.     lasso_model.fit(x, y) 
  29.     print '验证参数: ', lasso_model.best_params_ 
  30.  
  31.     y_hat = lasso_model.predict(np.array(x_test)) 
  32.     mse = np.average((y_hat - np.array(y_test)) ** 2)  # Mean Squared Error 
  33.     rmse = np.sqrt(mse)  # Root Mean Squared Error 
  34.     print mse, rmse 
  35.  
  36.     t = np.arange(len(x_test)) 
  37.     plt.plot(t, y_test, 'r-', linewidth=2label='Test'
  38.     plt.plot(t, y_hat, 'g-', linewidth=2label='Predict'
  39.     plt.legend(loc='upper right'
  40.     plt.grid() 
  41.     plt.show() 
  42.  

Advertising 完整版代码:

  1. #!/usr/bin/python 
  2. # -*- coding:utf-8 -*- 
  3.  
  4. import csv 
  5. import numpy as np 
  6. import matplotlib.pyplot as plt 
  7. import pandas as pd 
  8. from sklearn.model_selection import train_test_split 
  9. from sklearn.linear_model import LinearRegression 
  10.  
  11.  
  12. if __name__ == "__main__": 
  13.     path = '8.Advertising.csv' 
  14.     # # 手写读取数据 - 请自行分析,在8.2.Iris代码中给出类似的例子 
  15.     # f = file(path) 
  16.     # x = [] 
  17.     # y = [] 
  18.     # for i, d in enumerate(f): 
  19.     #     if i == 0: #第一行是标题栏 
  20.     #         continue 
  21.     #     dd = d.strip() #去除首位空格 
  22.     #     if not d: 
  23.     #         continue 
  24.     #     d = map(float, d.split(',')) #每个数据都变为float 
  25.     #     x.append(d[1:-1]) 
  26.     #     y.append(d[-1]) 
  27.     # print x 
  28.     # print y 
  29.     # x = np.array(x) #显示的更好看 
  30.     # y = np.array(y) 
  31.     # print x 
  32.     # print y 
  33.  
  34.     # # Python自带库 
  35.     # f = file(path, 'rb') 
  36.     # print f 
  37.     # d = csv.reader(f) 
  38.     # for line in d: 
  39.     #     print line 
  40.     # f.close() 
  41.  
  42.     # # numpy读入 
  43.     # p = np.loadtxt(path, delimiter=','skiprows=1
  44.     # print p 
  45.     # print p.shape 
  46.     # print p[1,2] 
  47.     # print type(p[1,2]) 
  48.  
  49.     # pandas读入 
  50.     data = pd.read_csv(path)    # TV、Radio、Newspaper、Sales 
  51.     # x = data[['TV', 'Radio', 'Newspaper']] 
  52.     x = data[['TV', 'Radio']] 
  53.     y = data['Sales'] 
  54.     # print x 
  55.     # print y 
  56.  
  57.     # 绘制1 
  58.     plt.plot(data['TV'], y, 'ro', label='TV'
  59.     plt.plot(data['Radio'], y, 'g^', label='Radio'
  60.     plt.plot(data['Newspaper'], y, 'mv', label='Newspaer'
  61.     plt.legend(loc='lower right'
  62.     plt.grid() 
  63.     plt.show() 
  64.  
  65.     # 绘制2 
  66.     plt.figure(figsize=(9,12)) #设置图的大小 宽9inch 高12inch 
  67.     plt.subplot(311) 
  68.     plt.plot(data['TV'], y, 'ro') 
  69.     plt.title('TV') 
  70.     plt.grid() 
  71.     plt.subplot(312) 
  72.     plt.plot(data['Radio'], y, 'g^') 
  73.     plt.title('Radio') 
  74.     plt.grid() 
  75.     plt.subplot(313) 
  76.     plt.plot(data['Newspaper'], y, 'b*') 
  77.     plt.title('Newspaper') 
  78.     plt.grid() 
  79.     plt.tight_layout() # 紧凑显示图片,居中显示 
  80.     plt.show() 
  81.  
  82.     x_train, x_test, y_train, y_test
  83.  = train_test_split(x, y, train_size = 0.75, random_state=1) #random_state 种子 
  84.     # print x_train, y_train 
  85.     linreg = LinearRegression() 
  86.     model = linreg.fit(x_train, y_train) 
  87.     print model 
  88.     print linreg.coef_ #系数 
  89.     print linreg.intercept_ #截距 
  90.  
  91.     y_hat = linreg.predict(np.array(x_test)) 
  92.     mse = np.average((y_hat - np.array(y_test)) ** 2)  # Mean Squared Error 
  93.     rmse = np.sqrt(mse)  # Root Mean Squared Error 
  94.     print mse, rmse 
  95.  
  96.     t = np.arange(len(x_test)) 
  97.     plt.plot(t, y_test, 'r-', linewidth=2label='Test'
  98.     plt.plot(t, y_hat, 'g-', linewidth=2label='Predict'
  99.     plt.legend(loc='upper right'
  100.     plt.grid() 
  101.     plt.show() 

Advertising 正则化 交叉验证相关代码:

线性回归多项式拟合:

  1. #!/usr/bin/python 
  2. # -*- coding:utf-8 -*- 
  3.  
  4. import numpy as np 
  5. from sklearn.linear_model import LinearRegression, RidgeCV 
  6. from sklearn.preprocessing import PolynomialFeatures 
  7. import matplotlib.pyplot as plt 
  8. from sklearn.pipeline import Pipeline 
  9. import matplotlib as mpl 
  10.  
  11.  
  12. if __name__ == "__main__": 
  13.     np.random.seed(0) # 指定种子 
  14.     N = 9 
  15.     x = np.linspace(0, 6, N) + np.random.randn(N) 
  16.     x = np.sort(x) 
  17.     y = x**2 - 4*x - 3 + np.random.randn(N) 
  18.     # print x 
  19.     # print y 
  20.     x.shape = -1, 1 
  21.     y.shape = -1, 1 
  22.     # print x 
  23.     # print y 
  24.  
  25.     model_1 = Pipeline([ 
  26.         ('poly', PolynomialFeatures()), 
  27.         ('linear', LinearRegression(fit_intercept=False))]) 
  28.     model_2 = Pipeline([ 
  29.         ('poly', PolynomialFeatures()), 
  30.         ('linear', RidgeCV(alphas=np.logspace(-3, 2, 100), fit_intercept=False))]) 
  31.     models = model_1, model_2 
  32.     mpl.rcParams['font.sans-serif'] = [u'simHei'] 
  33.     mpl.rcParams['axes.unicode_minus'] = False 
  34.     np.set_printoptions(suppress=True
  35.  
  36.     plt.figure(figsize=(9, 11), facecolor='w'
  37.     d_pool = np.arange(1, N, 1)  # 阶 
  38.     m = d_pool.size 
  39.     clrs = []  # 颜色 
  40.     for c in np.linspace(16711680, 255, m): 
  41.         clrs.append('#%06x' % c) 
  42.     line_width = np.linspace(5, 2, m) 
  43.     titles = u'线性回归', u'Ridge回归' 
  44.     for t in range(2): 
  45.         model = models[t] 
  46.         plt.subplot(2, 1, t+1) 
  47.         plt.plot(x, y, 'ro', ms=10zorder=N
  48.         for i, d in enumerate(d_pool): 
  49.             model.set_params(poly__degree=d) 
  50.             model.fit(x, y) 
  51.             lin = model.get_params('linear')['linear'] 
  52.             if t == 0: 
  53.                 print u'%d阶,系数为:' % d, lin.coef_.ravel() 
  54.             else: 
  55.                 print u'%d阶,alpha=%.6f,系数为:' % (d, lin.alpha_), lin.coef_.ravel() 
  56.             x_hat = np.linspace(x.min(), x.max(), num=100
  57.             x_hat.shape = -1, 1 
  58.             y_hat = model.predict(x_hat) 
  59.             s = model.score(x, y) 
  60.             print s, ' ' 
  61.             zorder = N - 1 if (d == 2) else 0 
  62.             plt.plot(x_hat, y_hat, color=clrs[i]
  63. lw=line_width[i], label=(u'%d阶,score=%.3f' % (d, s)), zorderzorder=zorder) 
  64.         plt.legend(loc='upper left'
  65.         plt.grid(True) 
  66.         plt.title(titles[t], fontsize=16
  67.         plt.xlabel('X', fontsize=14
  68.         plt.ylabel('Y', fontsize=14
  69.     plt.tight_layout(1, rect=(0, 0, 1, 0.95)) 
  70.     plt.suptitle(u'多项式曲线拟合', fontsize=18
  71.     plt.show() 

 不使用库:

BGD 与 SGD:

  1. # -*- coding: cp936 -*- 
  2. import numpy as np 
  3. import matplotlib.pyplot as plt 
  4.  
  5.  
  6. def linear_regression_BGD(x, y, alpha, lamda): 
  7.     m = np.alen(x) 
  8.     ones = np.ones(m) 
  9.     x = np.column_stack((ones, x)) 
  10.     n = np.alen(x[0]) 
  11.     theta = np.ones(n) 
  12.     x_traverse = np.transpose(x) 
  13.  
  14.     for i in range(1000): 
  15.         hypothesis = np.dot(x, theta) 
  16.         loss = hypothesis - y 
  17.         cost = np.sum(loss ** 2) 
  18.         print i, cost 
  19.         gradient = np.dot(x_traverse, loss) 
  20.         thetatheta = theta - alpha * gradient 
  21.     return theta 
  22.  
  23. def liear_regression_SGD(x, y, alpha, lamda): 
  24.     m = np.alen(x) 
  25.     ones = np.ones(m) 
  26.     x = np.column_stack((ones, x)) 
  27.     n = np.alen(x[0]) 
  28.     theta = np.ones(n) 
  29.     for j in range(1, m): 
  30.         hypothesis = np.dot(x[j], theta) 
  31.         loss = hypothesis - y[j] 
  32.         gradient = np.dot(loss, x[j]) 
  33.         thetatheta = theta - alpha * gradient 
  34.     return theta 
  35.  
  36.  
  37.  
  38. if __name__ == '__main__': 
  39.     N = 10 
  40.     # x = np.linspace(0, 10, N) + np.random.randn(N) 
  41.     # y =  3 * x + 5 + np.random.randn(N) 
  42.  
  43.     x = np.linspace(0, 10, N) + np.random.randn(N) 
  44.     y = 4 * x * x + 3 * x + 5 + np.random.randn(N) 
  45.     xx_square = x * x 
  46.     x_predict = np.column_stack((x, x_square)) 
  47.  
  48.     # theta = linear_regression_BGD(x_predict, y, 0.00001,0.1) # 批量梯度下降 
  49.     theta = liear_regression_SGD(x_predict, y, 0.0001, 0.1) # 随机梯度下降 
  50.     plt.plot(x, y, 'ro') 
  51.  
  52.     x = np.linspace(x.min(), x.max(), 10 * N) # 构建测试数据 
  53.     ones = np.ones(10 * N) 
  54.     # x_predict = np.column_stack((ones, x))x 
  55.     x_test = np.column_stack((ones, x, x * x)) 
  56.     y = np.dot(x_test, theta) 
  57.  
  58.     plt.plot(x, y, 'b-') 
  59.     plt.show() 

Regression代码:

  1. import numpy as np 
  2. import matplotlib.pyplot as plt 
  3. def regression(data, alpha, lamda ): 
  4.     n = len(data[0]) - 1 
  5.     theta = np.zeros(n) 
  6.     times = 1 
  7.     for i in range(times): 
  8.         for d in data: 
  9.             x = d[:-1] 
  10.             y = d[-1] 
  11.             h_theta = np.dot(theta, x) - y 
  12.             thetatheta = theta - alpha * h_theta * x + lamda * theta 
  13.         #print i,theta 
  14.     return theta 
  15. def preduceData(): 
  16.    x = [] 
  17.    y = [] 
  18.    for i in range(0, 50): 
  19.        x.append(i) 
  20.        y.append(3 * x[i] + np.random.random_sample()*3) 
  21.    data = np.array([x, y]).T 
  22.    theta = regression(data, 0.001, 0.1) 
  23.    return x, y, theta 
  24. def myplot(x, y, theta): 
  25.     plt.figure() 
  26.     plt.plot(x, y, 'go') 
  27.     plt.plot(x, theta * x, 'r') 
  28.     plt.show() 
  29. x, y, theta = preduceData() 
  30. myplot(x, y, theta) 

 

(futurehau)
本站文章除注明转载外,均为本站原创或编译欢迎任何形式的转载,但请务必注明出处,尊重他人劳动,同学习共成长。转载请注明:文章转载自:罗索实验室 [http://www.rosoo.net/a/202004/17641.html]
本文出处:博客园 作者:futurehau 原文
顶一下
(0)
0%
踩一下
(0)
0%
------分隔线----------------------------
相关文章
发表评论
请自觉遵守互联网相关的政策法规,严禁发布色情、暴力、反动的言论。
评价:
表情:
用户名: 验证码:点击我更换图片