引言：本文是《Python与机器学习应用--经济学视角》专栏的第一篇文章。正如专栏题目所示，笔者在介绍主流的机器学习算法的时候，会基于本人Stata/R语言以及经济学/计量经济学的背景进行。
这里还假设大家已经拥有一定的计量经济学基础，并且略懂Stata，R和Python。在写作的过程中，笔者会经常穿插着将主流的机器学习术语和计量经济学中的概念对比讲解；同时也会比较Stata和R语言，以及Python在实现同一个算法的优缺。
简单线性回归，也叫普通最小二乘法（Ordinary Least Square，OLS）不仅是计量经济学的基础模型，而是众多机器学习教材介绍的第一个算法。一句话，机器学习中的线性回归和计量经济学或统计学中的OLS，本质上是同一个东西，只是具体的术语和表达会稍有区别。
用Python介绍OLS，那当然要展示些特别的地方：1）手动计算损失函数值（SSR）；2）手动计算R方；3）基于损失函数，手动解出OLS的解释解。（注：基于最小化损失函数值求解，是机器学习算法通用的思想）
下面，我会基于Python的机器学习包sciki-learn，结合具体案例，逐步向大家介绍OLS的建模实现过程。以下案例用于披萨的价格预测——基于披萨的尺寸。我们先通过现有的数据，建立回归模型。并在新数据（披萨尺寸）的基础上，对披萨的价格作出预测。
1 准备工作。输入并查看数据
> import numpy as np #导入用于数组和矩阵运算的包numpy
> import matplotlib.pyplot as plt #导入用于画图的包matplotlib，并简写为plt
#输入数据
> X = np.array([[6],[8],[10],[14],[18],[19],[21],[23],[24],[25]]).reshape(-1, 1)
#reshape(-1, 1)表示转换为一列
> X #这里自变量用矩阵来存储。大小字母表示矩阵
> y = [7, 9, 13, 17.5, 18, 19.5, 21, 22, 24.3, 26]
> y #因变量用向量来储存。小写字母表示向量
#通过画图大概了解X和y的关系
> plt.figure() #定义一个空的图
> plt.title('Pizza diameter and price plot') #图的标题
> plt.xlabel('Diameter (inches)') #横坐标标题
> plt.ylabel('Prices (dollars)') #纵坐标标题
> plt.plot(X, y, 'k.') #k.表示画成散点图。如果不加则画成线型图
> plt.axis([0, 25, 0, 25]) #定义x轴和y轴的刻度长度
> plt.grid(True) #加上网格
> plt.show() #输出图像
披萨的大小和价格关系散点图
#2 基于现有数据，建立OLS模型
> from sklearn.linear_model import LinearRegression
> m1 = LinearRegression() #建立一个空的线型回归模型（实例）
> m1.fit(X, y) #用已有的训练数据集，拟合模型
#Out[61]: LinearRegression()
> m1.intercept_ #查看模型的截距项
#Out[62]: 2.6717877094972096
> m1.coef_ #查看模型的系数
#Out[63]: array([0.89632216])
#3 使用测试数据进行预测
> test_diameter1 = np.array([[22]])
> price_hat1 = m1.predict(test_diameter1)[0]
> price_hat1
#Out[66]: 22.390875232774672 #基于我们训练出的OLS模型，22寸的披萨应该卖22.39美金
#4 计算损失函数值，SSR
#在OLS中，损失函数就是我们的残差平方和SSR
> m1.predict(X) #全部的y的预测值
> ssr = np.mean( (m1.predict(X) - y)**2 )
> ssr
#Out[20]: 1.2323189013035378
#5 基于最小化损失函数，解出OLS的解释解
> import numpy as np
> X = np.array([[6],[8],[10],[14],[18],[19],[21],[23],[24],[25]]).reshape(-1, 1)
> x_bar = X.mean()
> x_bar
#Out[24]: 16.8
#计算自变量x的方差
> var_x = ((X-x_bar)**2).sum() / (X.shape[0]-1)
> var_x
#Out[27]: 47.733333333333334
#计算x和y的协方差
> y = [7, 9, 13, 17.5, 18, 19.5, 21, 22, 24.3, 26] #列表型
> y = np.array([7, 9, 13, 17.5, 18, 19.5, 21, 22, 24.3, 26]) #列表型转换为数组型
#np数组计算均值方差等会更方便
> y_bar = y.mean()
> y_bar
#Out[32]: 17.73
> cov_xy = np.multiply( (X-x_bar).transpose(), y-y_bar).sum() / (X.shape[0]-1)
> cov_xy
#Out[34]: 42.78444444444444
> coef = cov_xy/var_x
> coef
#Out[36]: 0.8963221601489757
> intercept = y_bar - x_bar*coef
> intercept
#Out[38]: 2.671787709497208
#和前面拟合模型的结果一致
#m1.intercept_ = 2.6717877094972096; m1.coef = 0.89632216
#6 评价模型 -- 手动计算R方
> X_train = np.array([[6],[10],[18],[21],[24]]).reshape(-1, 1) #抽取一半样本
> y_train = np.array([7, 13, 18, 21, 24.3]) #抽取一半样本
> X_test = np.array([[8],[14],[19],[23],[25]]).reshape(-1, 1) #剩一半样本
> y_test = np.array([7, 13, 18, 21, 24.3]) #剩一半样本
> y_test_bar = y_test.mean()
> m2 = LinearRegression()
> m2.fit(X_train, y_train)
> r2 = m2.score(X_test, y_test) #用函数计算R方
> r2
#Out[69]: 0.900523046779419
#注意，这里计算R2的思想和在拟合模型中计算会有所区别。这里的y和ybar用的是测试数据，yhat
#则基于训练数据拟合的模型所计算出来
> sst = np.mean((y_test - y_test_bar)**2)
> sst #Out[79]: 37.142399999999995
> y_test_hat = m2.predict(X_test) #基于训练数据拟合的模型，计算测试数据的yhat
> ssr = np.mean((y_test - y_test_hat)**2)
> ssr #Out[81]: 3.6948127873001106
> 1 - ssr / sst #Out[83]: 0.900523046779419。和前面算出来的r2一样
小结：和Stata/R对比起来，Python的算法实现显得更底层，更零碎，效率更低。但好处是让用户更清楚的了解OLS算法的细节。
对于线性回归模型，Python也有从计量经济学角度开发的包，比如statsmodels。statsmodels的封装程度更高，使用会让人感觉更像R语言（坦白说Python里面很多统计或计量包，就是参考或模仿R语言创建的）。
Ref
- sciki-learn 官方帮助手册,
contents - scikit-learn 0.23.2 documentationscikit-learn.org
-----全文结束---------------------------------------------------------------------------------
