Python机器学习与数据分析系列(2)-线性回归

最新推荐文章于 2024-08-13 01:54:14 发布

希铭

最新推荐文章于 2024-08-13 01:54:14 发布

阅读量3.3k

点赞数 3

分类专栏：机器学习数据分析线性回归文章标签：机器学习数据分析算法

本文链接：https://blog.csdn.net/hawkuu/article/details/77480525

版权

机器学习同时被 3 个专栏收录

3 篇文章 0 订阅

订阅专栏

数据分析

1 篇文章 0 订阅

订阅专栏

线性回归

1 篇文章 0 订阅

订阅专栏

1. 引言

这一章的主题是线性回归。我们按照这一顺序来展开教程：

什么是线性回归
如何用sklearn来实现线性回归的计算
线性回归的算法原理是什么
如何用python创建线性回归的算法
总结

2. 什么是线性回归

线性回归回归可以说是机器学习的入门算法，类似于学代码先学会如何打印“Hello World”，因此我们从实现线性回归入手。

先来说说线性回归的上延概念：机器学习。机器学习分为两种：一种是回归（regression）学习，一种是分类（classification）学习。为什么这里叫做“回归学习”和“分类学习”呢？因为机器学习算法的重点在于对数据集的学习，而不同的机器学习算法，只是实现这种学习的一种载体，因此机器学习与以往的统计分析有所同而有所不同。

（接着回来讲回归和分类。）

回归就是对一个数据集的趋势进行拟合，从而产生一个模型。
例如，我们手里有一段这样的数据集：

蓝色点代表数据集，它是一个二维点集，每个点代表一个样本（sample）。红色曲线代表线性回归拟合曲线。线性回归分析的结果就是创造这样一条拟合曲线，它反映了这个样本点集的整体”趋势”。

可以想象，当我们拥有了红色“拟合曲线”以后，就可以用于“预测”。例如，横坐标X轴（x-asix）给定一个数，那么就可以根据拟合曲线，得到纵坐标Y轴（y-axis）的值、

如果我们将X轴上的数当成自变量，Y轴上的数当成因变量，那么就相当于建立了从自变量-因变量的数学关系，上面的“预测”就可以翻译为：给定你一个样本的自变量的值，根据拟合曲线，就可以给出一个样本的因变量的值，从而实现对这个样本的“预测”。

分类就是对一个数据集中不同种类样本的分类，也产生一个模型。
例如，我们手里有这样一段数据：

黑色的点代表一类样本，红色的点代表另外一类样本，机器学习分类模型学习得到黑色、红色两类样本的分类模型。然后我们就可以将一个新的样本点（蓝色点），对它的归类进行“预测”。

由此可见，机器学习也是一门用于“预测”建模的学问，实现其“预测”能力的重点在于“学习”，而“学习”是基于数据集的。在后面的章节中，你会从代码层面了解这种学习机制。

然后我们也从表面上了解了什么是“线性回归”，这一章主要围绕它展开。

3. 用sklearn简单地实现的线性回归

3.1 获取数据

首先我们将需要的库安装一下：

$ sudo pip install numpy
$ sudo pip install scipy
$ sudo pip install sklearn
$ pip install matplotlib
$ pip install pandas
$ pip install quandl
# quandl 是个获取数据库的库，后面我们会用它来获取数据集

在粗略地了解了线性回归的原理之后，我们用sklearn的线性回归（linear_regression)函数解决一个股票价格预测的例子：

进入python命令行：

>> import quandl
>> df = quandl.get('WIKI/GOOGL')
>> print(df.head())
# df.head()是打印dataframe的头部信息
# df.tail()可以打印dataframe的尾部信息

输出得到：

              Open    High     Low   Close    Volume  Ex-Dividend  \
Date                                                                
2004-08-19  100.00  104.06   95.96  100.34  44659000            0   
2004-08-20  101.01  109.08  100.50  108.31  22834300            0   
2004-08-23  110.75  113.48  109.05  109.40  18256100            0   
2004-08-24  111.24  111.60  103.57  104.87  15247300            0   
2004-08-25  104.96  108.00  103.88  106.00   9188600            0   

            Split Ratio  Adj. Open  Adj. High  Adj. Low  Adj. Close  \
Date                                                                  
2004-08-19            1     50.000      52.03    47.980      50.170   
2004-08-20            1     50.505      54.54    50.250      54.155   
2004-08-23            1     55.375      56.74    54.525      54.700   
2004-08-24            1     55.620      55.80    51.785      52.435   
2004-08-25            1     52.480      54.00    51.940      53.000   

            Adj. Volume  
Date                     
2004-08-19     44659000  
2004-08-20     22834300  
2004-08-23     18256100  
2004-08-24     15247300  
2004-08-25      9188600

在这里，我们用quandl.get()获取了WIKI数据库中GOOGL这支股票的数据，这似乎是一个表格：row是不同天数，column是每天的股票指数，如开盘价、收盘价、最高价、成交量等。

让我们来看一下获取的数据格式是什么：

>> type(df)
<class 'pandas.core.frame.DataFrame'>

说明这是一个pandas的dataframe（数表）

另外，想要获取其它股票的信息，可以去Quandl的官网查找不同的’数据库/股票’（类似于’WIKI/GOOGL’)的获取方法。

3.2 清洗数据

我们先来看一看这段数据：

            Open    High     Low   Close    Volume  Ex-Dividend  \
Date                                                                
2004-08-19  100.00  104.06   95.96  100.34  44659000            0   
2004-08-20  101.01  109.08  100.50  108.31  22834300            0   

            Split Ratio  Adj. Open  Adj. High  Adj. Low  Adj. Close  \
Date                                                                  
2004-08-19            1     50.000      52.03    47.980      50.170   
2004-08-20            1     50.505      54.54    50.250      54.155    

            Adj. Volume  
Date                     
2004-08-19     44659000  
2004-08-20     22834300

Adj. Open/High/Low/Close/Volume 代表调整后的数据，至于原理是什么，其实我也不大清楚……但不管怎样，我们可以将它们看成是由许多样本组成的数据集（横坐标Date中的一条代表一个样本，纵坐标的每一个条目代表这个样本的属性）。

我们将Open/High/…/Adj.Open/Adj.High/…/Adj.Volume称为一个样本的属性、或是特征（features）。

那么这段数据集（dataset）有多少个样本（sample），每个样本有多少个特征（feature）呢？我们可以用：

>> print(len(df))
3272

>> print(len(df.columns))
12

3272，这个数显然会增长下去，代表从数据库开始记录到最近发生的一天。另外，len()这个函数非常有用，务必熟练掌握之，它和type()一样，是最常用的命令，会经常打着玩儿。

（下面回到主题）

在pythonprogramming.net的原文中，作者认为Open、High……这些feature太多了，Adj. Open、Adj. High这些指标不够好，需要做一下处理，来更好地表达股票的价格变动（至于其中缘由，也未详说。如果有时间可以深究一下，也欢迎大家来补充），我们姑且当成一种数据清洗的过程：

df = df[['Adj. Open', 'Adj. High', 'Adj. Low', 'Adj. Close', 'Adj. Volume']]

处理完以后，你可以再一次输入df看看结果：

>> df

这样就只留下了'Adj. Open', 'Adj. High', 'Adj. Low', 'Adj. Close', 'Adj. Volume'这几个column。
这里请注意一下：如果想调用dataframe中具体哪个feature的值，就用df['feature']即可，而要选择多个features，就得用df[['feature1', 'feature2']]才行，不得不顺着pandas的规矩来。

df['HL_PCT'] = (df['Adj. High'] - df['Adj. Low'])/(df['Adj. Close']) * 100.0

df['PCT_change'] = (df['Adj. Close'] - df['Adj. Open'])/(df['Adj. Open']) * 100.0

然后：

df = df[['Adj. Close', 'HL_PCT', 'PCT_change', 'Adj. Volume']]

在这里，原作者将原来12个features缩减为4个feature，降低了样本的特征维度。其中缘由，还需请高手来解释，我就按下不表了。

再来看一看结果：

            Adj. Close    HL_PCT  PCT_change  Adj. Volume
Date                                                     
2004-08-19      50.170  8.072553    0.340000     44659000
2004-08-20      54.155  7.921706    7.227007     22834300
2004-08-23      54.700  4.049360   -1.218962     18256100
2004-08-24      52.435  7.657099   -5.726357     15247300
2004-08-25      53.000  3.886792    0.990854      9188600

除了收盘价和成交量，其它两个feature都变成了百分数。这是因为上面乘了100.0的缘故。后面我们会将Adj. Close作为因变量，也就是预测的目标。将Adj. Close, HL_PCT, PCT_change, Adj. Volume作为自变量，也就是线性回归的features，建立一个对Adj. Close（这里有点怪），也就是收盘价进行预测的模型：

import math
forecast_col = 'Adj. Close'
# 将Adj. Close作为我们的预测column
df.fillna(value = -99999, inplace=True)
# df.fillna 的作用是将dataframe中的N/A值换成-99999，机器学习算法无法读取数据中的N/A值
# 且会将-99999这个值视为一个溢出值不做计算
# inplace=True是指替换原数表，相当于：df = df.fillna(value = -99999)
forecast_out = int(math.ceil(0.01 * len(df)))
# 这里我们使用0.01*len(df)来获得dataframe长度的1%，也就是取1%的sample作为预测
# 我算出来的forecast_out = 33，所以后面就用33了
# math.ceil()的作用是返回一个最大整数。ceil指天花板，floor指地板。而math.floow()指取最小整数
df['label'] = df[forecast_col].shift(-forecast_out)
# 这是很关键的一步，有助于理解如何实现“预测”
# df.shift()是对dataframe进行数据移动，df['forecast_out'].shift(-forecast_out)就是
# 将df['forecast_out']这个column的数据向上移动33个条目。这样，所有的
# sample都相当于提前“预知”了33天后的收盘价。请注意，forecast_col代表的
# 是‘Adj. Close'这个column，forecast_out代表的是我们设定预测的天数，也就是总条目
# 数的1%。

这里我要插一句，讲讲这个df[‘label’]。为什么要新建一个column，叫df[‘label’]呢？在机器学习中，分为两种数据学习方法：有监督学习(supervised learning)和无监督学习(unsupervised learning)。

有监督学习就是指用于机器学习的数据集中含有所要预测的结果。用本文第二节“什么是线性回归”中的两个例子来说：它们二者都是属于有监督学习。对于前者：回归问题，数据集给出了因变量 $y$ 的值，也就是应该预测的结果。对于后者：分类问题，数据集给出了两类点集所属的分类，也一样是应该预测的结果。

那么什么是无监督学习呢？无监督学习就是只给出数据集的features（自变量），但不给出因变量。要想实现自变量向因变量的预测，需要你指导机器学习算法自己去找规律。

这样就分清了有监督学习和无监督学习的区别，而有监督学习的那个已知的因变量，就成为标签（label），所以在这里的那个已知的提前预支的收盘价格（Adj. Close)的column，就是df[‘label’]了。

3.3 创建数据集

前一步我们对数据进行了清洗和整理，现在它完全可以用于机器学习了（去除了N/A值，也减少了features的维数）。下面我们创建一个数据集。

现在又要来科普了。对于线性回归而言，它需要有自变量，也就是features；由于它是一个有监督学习，并且我们要用它拿来预测，因此也需要有因变量，也就是label；那么还需要有什么呢？

前面已经提过，我们要用数据集对模型进行“学习”，但光学习是不够的，还需要拿出一部分数据对它进行“测试”，这样我们才知道这个模型是不是合格。评价模型是否合格的指标称为“泛化能力”，但是我们现在还不需要讨论这个概念。现在只需要知道，数据集需要按照一定比例（一般是1:9、2:8、3:7，训练样本占多数，测试样本占少数）划分成训练集合测试集：

import numpy as np
from sklearn import cross_validation
# 这里我们用到了numpy，numpy是将dataframe数据类型变成机器学习算法可以使用的数据类型的工具
X = np.array(df.drop(['label'], 1))
# 这里我们使用了np.array()，就是将pandas.dataframe数据转变为numpy.array数据类型
# 可以使用type(X)来查看获得的数据类型，这很重要。随时掌握你处理的数据类型是python编程的技巧
# 这里再一次使用了df.drop()，其中的"1"代表舍弃整个纵轴，"0"就代表舍弃横轴。
# 现在我们已经学会了三个类似的操作：df.dropna(), df.drop(), df.fillna()，它们非常好用
X = preprocessing.scale(X)
X_lately = X[-forecast_out:]
# 这一步也是需要关注的，之前我们创建了'label'的column，并把它们往上移动了33个条目
# 现在，我们将最末尾的33个sample当成预测的sample使用，将前面的sample用来训练。
# 这里的X[-forecast_out:]是python的一种语法，用来表示“倒数33个sample 至 最末尾”
X = X[:-forecast_out]
# 这里X和X_lately的顺序也不要写反了，不然会错（想想为什么会错?）
df.dropna(inplace=True)
y = np.array(df['label'])
# 这一步似乎挺取巧的，因为dataframe已经去除了N/A的值，显然df['label']的sample数与X的sample数一致
# 现在我们有了自变量X和因变量y，最后一步是分割数据集：
X_train, X_test, y_train, y_test = cross_validation.train_test_split.(X, y, test_size=0.8)
# 这里使用了sklearn.cross_validation.train_test_split()这个工具
# 注意是X_train, X_test, y_train, y_test顺序不要写错了，写错后面就都乱了
# test_size=0.8就是训练集/测试机分割比例2:8

下面对上述代码做一下注释：

在取出X后，我们需要对X中的数据做标准化预处理：

from sklearn import preprocessing
X = preprocessing.scale(X)

于是你一定会问什么是标准化？其实我也不大清楚…姑且就认为是一种对数据的预处理手段吧。
说不定你想一巴掌扇过来。

还是开一个小副本来解释一下吧：

———-以下是副本———-
标准化可以理解为将数据转化为标准分布的数据，例如：均值为0，方差为1，即标准正态分布。它将每种features转化成统一标准的数据分布。

我们来看看X的数据类型（注意，为了便于观察，我没有转化成numpy.array，而是使用了pandas.dataframe）：

            Adj. Close    HL_PCT  PCT_change  Adj. Volume
Date                                                     
2004-08-19   50.322842  8.072956    0.324968   44659000.0
2004-08-20   54.322689  7.921706    7.227007   22834300.0
2004-08-23   54.869377  4.049360   -1.227880   18256100.0
2004-08-24   52.597363  7.657099   -5.726357   15247300.0
2004-08-25   53.164113  3.886792    1.183658    9188600.0

可以发现，Adj. Close的数比HC_PCT, PCT_change的数要大好几倍，而Adj. Volume更是要超出它们很多。这在线性回归里会出现什么问题呢？就是很大的数在被分析时，其重要性可能会被看得很重要，而那些很小的数（如HL_PCT, PCT_change）的重要性却难以显现。因此我们需要将这四个features的数的数据分布整齐划一地统一起来：

x_scaled = preprocessing.scale(X)
print(x_scaled)

打印出：

array([[-1.4951861 ,  4.19528241,  0.23047628,  4.37424705],
       [-1.4767006 ,  4.08811581,  4.76739635,  1.75828185],
       [-1.47417405,  1.34439467, -0.79025756,  1.20952691],
       ..., 
       [ 2.713603  ,  0.47183083,  0.73251219, -0.65346554],
       [ 2.60642918,  0.14193002, -0.91798303, -0.59725418],
       [ 2.56880974, -0.35185968, -0.97958564, -0.71674364]])

是不是几个features的数值差不多了？

# 再试试查看均值和标准差
x_scaled.mean(axis=0)
x_scaled.std(axis=0)

好像均值并不是等于0的，尴尬了。有没有人可以解释以下子？

———-副本结束———-

这样你大概就能明白预处理的意义了：就是让各种features的数据分布方式接近一点，不会让哪一家坐得特别大。

然后我们还得再做一个步骤：

df.dropna(inplace=True)
# df.dropna()就是删除数表中N/A的sample。这个操作和df.fillna()不同，它会将出现N/A的sample整个删除

这条操作需要和上一段联系起来看。前一步操作为了创造对收盘价的“预测”，将df[‘forecast_col’]也就是df[‘Adj. Close’]整个向上移动了forecast_out条目，那么再最后的forecast_out条目上将会产生N/A数据：

            Adj. Close    HL_PCT  PCT_change  Adj. Volume  label
Date                                                            
2017-08-11      930.09  1.305250    0.690693    1589808.0    NaN
2017-08-14      938.93  0.698135   -0.014908    1140212.0    NaN
2017-08-15      938.08  0.685443   -0.313486    1006064.0    NaN
2017-08-16      944.27  1.044288    0.320850    1329301.0    NaN
2017-08-17      927.66  1.743096   -1.621507    1653779.0    NaN

这个操作就是删除这些N/A条目用的。

3.4 拟合（训练）、测试、打印结果

其实数据分析的大多数时间都在清洗、整理和研究优化方法，而最终训练、测试、打印则是水到渠成的事情，因为机器会帮你完成剩下的工作（前提是代码得调通是吧）：

from sklearn.linear_model import LinearRegression
clf = LinearRegression(n_jobs=-1) # n_jobs=-1指的是使用CPU多线程计算，这个功能不是每个算法都有
clf.fit(X_train, y_train)
confidence = clf.score(X_test, y_test) # 打印结果时使用了测试集：X_test, y_test
print(confidence)

最后一步是很简单的，我的结果大概是97%左右。那么这个97%代表什么意思呢？说了半天并没有解释现行回归是什么东西……？显然问题还有很多，我们先继续下去，对这个问题在下一节再做解释。

现在我们得到一个线性回归模型，模型的置信度（confidence）是97%，姑且认为是一个比较好的值吧（0~1，越接近1越好）。

还记得我们之前的目的是要预测最后-forecast_out天数的收盘价吧？现在我们来试试看：

forecast_set = clf.predict(X_lately)
print(forecast_set)

得到:

[  934.88742286   949.79864394   944.81590092   958.94886778   969.46911086
   971.53398824   985.88475057   986.08720328   995.40576051   993.58328027
  1005.62777954  1011.16308898  1009.56121011  1012.27153788  1015.74598931
   984.76959286   981.91499691   966.33054459   976.4747573    960.86385007
   964.22508578   964.56957141   957.21945847   963.85765919   963.54702951
   961.81564414   957.87539224   939.18417967   948.00782809   956.85739635
   955.801469     962.2408203    943.23461201]

这样一个数组。这实际上是对最后33天的价格预测。

我们来和真实的数据比较一下：

df_2 = quandl.get('WIKI/GOOGL')
df_2 = df_2['Adj. Close']
print(np.array(df_2[-33:]))

[ 919.46    932.26    927.69    940.81    951.      953.53    967.66
  968.7175  976.91    975.96    986.95    992.77    992.19    993.84
  998.31    969.03    965.31    952.51    958.33    945.5     946.56
  947.64    940.3     945.79    945.75    944.19    940.08    923.59
  930.09    938.93    938.08    944.27    927.66  ]

貌似结果差别还蛮大的……可能线性回归真的只能反映一种“趋势”吧，至少我能做股评家了~

3.5 怎么让妈妈相信你会了？作个图

好吧，我承认标题有一点雷人。但是吧，很多情况下我们码了半天，真实世界却完全不知道我们在干点什么，最后只招得一个“码农找不到女朋友”之类的尴尬名声。所以，即使我们现在还不懂线性回归、还不知道分析了点啥、结果还不怎么好，但起码还能做一张炫酷的图吧？

下面我们将使用强大的matplotlib来做一个股票收盘价格趋势预测图。先来想想我们要做一张什么图：显然是两个曲线的叠加：一条曲线的横坐标是时间，范围是[:-33]，纵坐标是样本里的收盘价（也就是Adj. Close)。另一条曲线的横坐标是时间，时间从末尾33天开始，纵坐标是收盘价，是我们预测出来的收盘价。

现在两个曲线的纵坐标已经有了，还差第二条曲线的横坐标，我们得先费尽心思生成这个横坐标：

import datetime

df['Forecast'] = np.nan # 先生成一个新的column，数值为N/A
last_date = df.iloc[-1].name # df.iloc[]是切片，后面细说
last_unix = np.int(last_date.value/10**9) # 这里涉及到关于“时间”的知识

one_day = 86400 # 一天中的秒=24*60*60=86400 seconds
next_unix = last_unix + one_day

for i in forecast_set:
        next_date = datetime.datetime.fromtimestamp(next_unix) # 这里又涉及到了时间的知识
        next_unix += 86400
        df.loc[next_date] = [np.nan for _ in range(len(df.columns)-1)]+[i]
        # 说实话我最讨厌这种“一点都不直观的”代码了，后面我们会简化它，让它说人话

说实话，上面那段代码我第一次看的时候完全懵逼了，根本不知道在说啥。这里面涉及到几个知识点需要了解：

先来说说python的切片：

>> df.iloc[-1]
Adj. Close     9.296800e+02
HL_PCT         1.655355e+00
PCT_change    -1.515906e+00
Adj. Volume    2.185444e+06
label          9.276600e+02
Name: 2017-06-30 00:00:00, dtype: float64

打印了数表的最后一个sample（注意这里的数表已经是被df.dropna()后的数表了）。

>> df.loc['2017-06-30']
Adj. Close     9.296800e+02
HL_PCT         1.655355e+00
PCT_change    -1.515906e+00
Adj. Volume    2.185444e+06
label          9.276600e+02
Name: 2017-06-30 00:00:00, dtype: float64

实际上，二者都是切片，也就是从dataframe中获取一个sample。iloc中的 $i$ 指index，因此需要输入序号索引，比如[-1]就是最后一个。loc就不能输入序号了，需要输入dataframe中的row，也就是行标签，例如最后一天就是loc[‘2017-06-30’]。这样就搞清楚了df.loc[]和df.iloc[]的用法。请注意，用的时候不要写成df.iloc()和df.loc()了，二者都必须输入[索引]。

接下来说说“timestamp”:

Unix时间戳(Unix timestamp)，或称Unix时间(Unix time)、POSIX时间(POSIX time)，是一种时间表示方式，定义为从格林威治时间1970年01月01日00时00分00秒(北京时间1970年01月01日08时00分00秒)起至现在的总秒数。Unix时间戳不仅被使用在Unix系统、类Unix系统中(比如Linux系统)，也在许多其他操作系统中被广泛采用。

在上文的代码中，使用了df.iloc[index]来提取每个sample，我们来看一下：

>> df.iloc[-1]
Adj. Close     9.296800e+02
HL_PCT         1.655355e+00
PCT_change    -1.515906e+00
Adj. Volume    2.185444e+06
label          9.276600e+02
Name: 2017-06-30 00:00:00, dtype: float64

显然它就是：

Adj. Close	HL_PCT	PCT_change	Adj. Volume	label
name

我们接着提取：

>> df.iloc[-1].name
Timestamp('2017-06-30 00:00:00')
>> type(df.iloc[-1].name)
<class 'pandas._libs.tslib.Timestamp'>

说明这是一个pandas的时间戳“timestamp”。

我们知道，timestamp是一个从1970年到现在的秒数（严格的说是纳秒），这个值怎么提取？

timestamp = df.iloc[-1].name.value//10**9
# 由于是纳秒，需要将它转换成秒

还需要补充的知识是“datetime”。事实上，python中有很多计量时间的办法，datetime和timestamp是其中的两个。具体可以参考这篇文章。你会对它们之间的转化有所了解。总之，datetime是另一种表示时间的方法，它更加直观（但会出现一些问题，这个后面会有所涉及）。为什么要提到datetime呢？因为df.loc[]不能接受数字，只能接受字符串，所以这里需要使用df.loc[datetime]。实际上对这个问题我了解不深，只是浅尝辄止地认识，希望有人可以来补充。

这样综合起来看，就可以明白上面的代码了，先来复述一下：

df['Forecast'] = np.nan # 先生成一个新的column，数值为N/A
last_date = df.iloc[-1].name # df.iloc[]是切片，后面细说
last_unix = np.int(last_date.value/10**9) # 这里涉及到关于“时间”的知识

one_day = 86400 # 一天中的秒=24*60*60=86400 seconds
next_unix = last_unix + one_day

for i in forecast_set:
        next_date = datetime.datetime.fromtimestamp(next_unix) # 这里又涉及到了时间的知识
        next_unix += 86400
        df.loc[next_date] = [np.nan for _ in range(len(df.columns)-1)]+[i]

翻译成人话就是：

在dataframe中用df.iloc[-1]切片，取最后一天的时间为timestamp()时间（此时是一个秒）
将这个时间加上一天
做一个循环：先将时间戳timestamp转换成datetime, 然后加一天， forecast_set中的数据按照每天增加的循环填入到f.loc[datetime]中去。
使用df.loc[datetime]使用一个for循环先生成一个N/A数组

注意，4当中的for循环，后面再讨论。综合来看，这样的处理虽然不易于理解，总体上还是比较精妙的，充满了程序思维。

然而，这样处理时间是有问题的，你看出来了吗？我再开一个副本来解释以下，无关线性回归的主题（感觉越跑越远了）

—————–副本2：打倒时差的魔鬼————-

其实小标题已经剧透了，实际上刚才的处理方式存在一个时差的问题，我们来看一下：

>> df.iloc[-1].name
Timestamp('2017-06-30 00:00:00')
>> df.iloc[-1].name.value//10**9
1498780800
# 说明2017-06-30 00:00:00的时间戳是1498780800

我们用datetime创造一个时间来对比一下：

>> date_time = datetime.datetime(2017,06,30,0,0)
>> date_time
datetime.datetime(2017, 6, 30, 0, 0)
# 再将datetime转换成timestamp
>> timestamp = time.mktime(date_time.timetuple())
>> timestamp
1498752000.0

两个数之间相差了28800，28800刚好是8个小时。（以下是推测，不完全肯定）由于pandas采用的是’UTC, Universial Time Coordinated’，也就是“世界统一时”的缘故。pandas使用了UTC的位于0点的时间，而我们使用datetime.datetime()创建的时间则是基于+8时区的时间，因此产生了8个小时的时差。如果你放任这个错误，那么执行代码df.loc[next_date]会发生错误，显然数表中不存在这样类似2017-06-30-08:00:00这样的时间。

我们再来印证一下：

>> date_time = pd.Timestamp('2017-06-30 00:00:00', tz='Asia/Shanghai')
>> print(date_time.value)//10**9
1498752000

我找到的解决方法是将数表的时区统一调整时间为+8区，也就是’Asia/Shanghai’这个时区，可以详见这个文档：

df_modified = df.tz_localize('Asia/Shanghai')

这样数表的时间就更正了。

———副本结束———

下面我们再来看一下

[np.nan for _ in range(len(df.columns)-1)]

这个循环，实际上就是创造一个数组，它可以改写为：

list = []
for _ in range(len(df.columns)-1):
        list.append(np.nan)

这种for循环的简写方法在后面还会用到多次，请特别注意。下面重新贴一下代码：

df['Forecast'] = np.nan # 先生成一个新的column，数值为N/A
last_date = df.iloc[-1].name # df.iloc[]是切片，后面细说
last_unix = np.int(last_date.value/10**9) # 这里涉及到关于“时间”的知识

one_day = 86400 # 一天中的秒=24*60*60=86400 seconds
next_unix = last_unix + one_day

for i in forecast_set:
        next_date = datetime.datetime.fromtimestamp(next_unix) # 这里又涉及到了时间的知识
        next_unix += 86400
        df.loc[next_date] = [np.nan for _ in range(len(df.columns)-1)]+[i]

现在我们将预测结果forecast_set补入了dataframe之中，只需要绘制df[‘Adj. Close’]就获得了股票收盘价格曲线，绘制df[‘Forecast’]就可以获得剩下的预测部分：

from matplotlib.pyplot as plt
from matplotlib import style

style.use('ggplot') # 这纯粹是为了更好看，以后可以慢慢研究
df['Adj. Close'].plot()
df['Forecast'].plot()
plt.xlabel('Date')
plt.ylabel('Price')
plt.legend(loc=2) # 这是为了曲线说明的位置，你可以改变loc的数字看看效果
plt.show()

这里我们使用了matplotlib.pyplot.plot()的用法，它的作用是将x轴、y轴的点集绘制成连续的曲线。最后要使用plt.show()将图像显示出来。

4. 线性回归的原理和实现

前面两节我们了解了线性回归的使用方法，现在来看一看它的原理。

假如我们面对一个二维点集组成的数据：
4-1

黄色曲线是线性回归曲线，代表了数据集的“趋势”，这很直观。

不过，当数据集是这样的时候：
4-2

情况就比较复杂了。我们无法直观地划出这个“回归曲线”，也没有办法定义什么样的“回归曲线”的好与坏。

所以仔细回来看看这条曲线（实际上就是直线），用初等数学的知识可以知道，它可以用一个函数来表达：

y = m x + b

$y = mx + b$

即得到一个函数关系：

y = f (x | m, b)

$y = f(x|m,b)$

在这里，x就是数据集的“features”，也是自变量。y就是数据集的“label”，也是因变量，y也是我们需要预测的结果。而 $m,b$ 是这个函数的参数，即组成模型的参数（在这个直线方程中，m就是斜率，b就是截距），带有几何含义。但这里features是一维的，如果features有很多唯，也就很难从几何上理解w和b了。

如何得到最能代表上述数据集的直线？如果feature是一维的，label也是一维的话，根据推导可以得到：

m = x ⎯ ⎯ \cdot y ⎯ ⎯ - x y ⎯ ⎯ ⎯ ⎯ ( x ⎯ ⎯ ) 2 - x 2 ⎯ ⎯ ⎯ ⎯

$m = \frac{\overline{x} \cdot \overline{y} - \overline{xy}}{(\overline{x})^2 - \overline{x^2}}$

得到了m的公式，最佳截距的公式可以反推得到：

b = y ⎯ ⎯ - m x ⎯ ⎯

$b = \overline{y} - m\overline{x}$

在这个过程我们省去了繁琐的数学证明过程（这篇文章只是为了熟练机器学习的代码部分，而非数学理论部分）。

现在我们用python写一下这个公式，并且创造几个样本，并用这个样本来生成一个拟合曲线：

import numpy as np

def best_fit_line(xs, ys):
        intercept  = (np.mean(xs)*np.mean(ys) - np.mean(xs*ys))/((np.mean(xs))**2 - np.mean(xs**2))
        slope = np.mean(ys) - intercept * np.mean(xs)

        return intercept, slope

xs = np.array([1,2,3,4,5],dtype=np.float64)
ys = np.array([5,4,6,5,6],dtype=np.float64)

x_predict = np.array([7])

import matplotlib.pyplot as plt
from matplotlib import style

intercept, slope = best_fit_line(xs,ys)

style.use('ggplot')
plt.scatter(xs, ys, color='b', label='sample')
# plt.scatter()是画散点图的意思
plt.plot(xs, intercept * xs + slope, color='r', label='linear_regression')
plt.scatter(x_predict, intercept * x_predict + slope, color='g', label='prediction')
plt.legend(loc=4)
plt.show()

那么如何评价这个曲线拟合得好还是坏呢？我们用两个统计评价指标来衡量：

SE（squared error，方差）

S E = \sum (y ˆ - y) 2

$SE = \sum(\widehat{y}-y)^2$

显然SE就是统计每个样本点的y值(y_orig)与回归曲线的纵坐标(ys_line)差值的求和（也可以求平均，效果都差不多）：
4-3

$R^{2}$ ，也称为确定系数（coefficient of determination)

R 2 = 1 - S E y ˆ S E y ⎯ ⎯

$R^2 = 1 - \frac{SE\widehat{y}}{SE\overline{y}}$

这里的 $\overline{y}$ 指拟合曲线。

显然我们也可以用python将它们写出来：

import numpy as np

def squared_error(ys_line, ys_orig):
        return np.sum((ys_line - ys_orig)**2)

def coefficient_of_determination(ys_line, ys_orig):
        squared_error_line = squared_error(ys_line, ys_orig)
        mean_line = [np.mean(ys_orig) for y in ys_orig]
        # 显然我们也可以这样写
        # mean_line = []
        # for y in ys_orig:
        #       mean_line.append(np.mean(ys_orig)
        squared_error_mean_line = squared_error(ys_orig, mean_line)
        return 1 - squared_error_line/squared_error_mean_line

通过代码我们可以进一步了解：

方差实际上就是sample原点与曲线取点的差的平方。
均线 $\overline{y}$ 是一条直线，也是一个由相同点组成的点集（这里同样使用了 $[for _ in (condition)]$ 的语法来简化一个for循环生成数组的例子。

这样我们就可以计算出上面给出例子的两个曲线评价指标：

import numpy as np

def best_fit_line(xs, ys):
        intercept  = (np.mean(xs)*np.mean(ys) - np.mean(xs*ys))/((np.mean(xs))**2 - np.mean(xs**2))
        slope = np.mean(ys) - intercept * np.mean(xs)

        return intercept, slope

xs = np.array([1,2,3,4,5],dtype=np.float64)
ys = np.array([5,4,6,5,6],dtype=np.float64)

x_predict = np.array([7])

import matplotlib.pyplot as plt
from matplotlib import style

intercept, slope = best_fit_line(xs,ys)

style.use('ggplot')
plt.scatter(xs, ys, color='b', label='sample')
# plt.scatter()是画散点图的意思
plt.plot(xs, intercept * xs + slope, color='r', label='linear_regression')
plt.scatter(x_predict, intercept * x_predict + slope, color='g', label='prediction')
plt.legend(loc=4)
plt.show()

import numpy as np

def squared_error(ys_line, ys_orig):
        return np.sum((ys_line - ys_orig)**2)

def coefficient_of_determination(ys_line, ys_orig):
        squared_error_line = squared_error(ys_line, ys_orig)
        mean_line = [np.mean(ys_orig) for y in ys_orig]
        # 显然我们也可以这样写
        # mean_line = []
        # for y in ys_orig:
        #       mean_line.append(np.mean(ys_orig)
        squared_error_mean_line = squared_error(ys_orig, mean_line)
        return 1 - squared_error_line/squared_error_mean_line

regression_line = [intercept * x + slope for x in xs]
# 这里只需要生成回归曲线的y值，存为一个list
# 显然这个for循环也是可以还原的。

SE = squared_error(regression_line, ys)
R2 = coefficient_of_determination(regression_line, ys)
print('SE = %f, R2 = %f' % (SE,R2))

下面我们上一点难度，不仅要创造几个样本，试着随机生成上百个样本，再拟合生成线性回归曲线：

import random

def create_dataset(number, variance, step, correlation=False):
        val = 1
        ys = []
        for _ in range(number): # 我觉得学会其中的编程思想比较重要
            y = val + random.randrange(-variance, variance)
            # 这一步是关键
            ys.append(y)
            if correlation == 'pos': # 生成递增的数据
                val += step
            elif correlation == 'neg': # 生成递减的数据
                val -= step

        xs = [i for i in range(len(ys))] # xs在这里就是用来打酱油的

        return np.array(xs, dtype=np.float64),np.array(ys, dtype=np.float64)

思路比较清晰，似乎没什么可补充的。现在我们学会了random.randrange()这个方法。注意最后生成出来的都是list。

现在我们将之前的代码合起来，就可以做实验了：

import numpy as np

def best_fit_line(xs, ys):
        intercept  = (np.mean(xs)*np.mean(ys) - np.mean(xs*ys))/((np.mean(xs))**2 - np.mean(xs**2))
        slope = np.mean(ys) - intercept * np.mean(xs)

        return intercept, slope

import random

def create_dataset(number, variance, step, correlation=False):
        val = 1
        ys = []
        for _ in range(number):
            y = val + random.randrange(-variance, variance)
            ys.append(y)
            if correlation == 'pos':
                val += step
            elif correlation == 'neg':
                val -= step
        xs = [i for i in range(len(ys))]

        return np.array(xs, dtype=np.float64),np.array(ys, dtype=np.float64)

xs, ys = create_dataset(100,10,1,'pos')

import matplotlib.pyplot as plt
from matplotlib import style

intercept, slope = best_fit_line(xs,ys)

style.use('ggplot')
plt.scatter(xs, ys, color='b', label='sample')
# plt.scatter()是画散点图的意思
plt.plot(xs, intercept * xs + slope, color='r', label='linear_regression')
# plt.scatter(x_predict, intercept * x_predict + slope, color='g', label='prediction')
plt.legend(loc=4)
plt.show()

import numpy as np

def squared_error(ys_line, ys_orig):
        return np.sum((ys_line - ys_orig)**2)

def coefficient_of_determination(ys_line, ys_orig):
        squared_error_line = squared_error(ys_line, ys_orig)
        mean_line = [np.mean(ys_orig) for y in ys_orig]
        # 显然我们也可以这样写
        # mean_line = []
        # for y in ys_orig:
        #       mean_line.append(np.mean(ys_orig)
        squared_error_mean_line = squared_error(ys_orig, mean_line)
        return 1 - squared_error_line/squared_error_mean_line

regression_line = [intercept * x + slope for x in xs]

SE = squared_error(regression_line, ys)
R2 = coefficient_of_determination(regression_line, ys)
print('SE = %f, R2 = %f' % (SE,R2))

怎么做呢？调整create_dataset(number, variance, step, correlation=False)中的variance就可以控制数据的“杂乱程度”，增加variance（方差），R2就越接近0，SE值增加，曲线拟合得就越不好。减小方差，R2就越接近1，SE值降低，曲线就拟合得越好。

5. 小节

由于忽略了数学推导部分，感觉有一点虎头蛇尾，争取以后再来补上（但我相信一定会损失不少读者）。

这篇文章作为机器学习系列的开篇，加入了一些我对机器学习的总体理解，导致篇幅过长，有时候也有些偏题，这将在后面修正。在这里我们大致来总结一下这篇文章：

机器学习（machine learning）总体上说就是建立一个从自变量到因变量映射的函数模型：

y = f (x | θ)

$y = f(x|\theta)$

结合股票预测的例子可以发现，影响股票收盘价的features就是这里的 $x$ ，也就是自变量。收盘价（预测的目标），就是这里的 $y$ 。机器学习会利用已有的sample组成的数据集dataset，机器学习通过对数据集的学习（这篇文章没有体现出来，因为我们直接得出了线性回归的模型参数求值方法），来得到一个良好的 $\theta$ ，从而建立起这个映射函数模型。对于一个新的sample，就可以通过输入它的feature，得到一个预测的 y <script type="math/tex" id="MathJax-Element-18">y</script>。

线性回归，就是这样一个模型，模型的参数已经定好了（对于二维数据，x为feature，y为label，我们用了intercept和slope作为它的参数）。我们使用一些方法去评价这个模型（否则就难以衡量模型的好坏），这里使用了SE（方差）和R2（确定系数）。通过生成随机数据可以发现，评价指标可以大致反映出模型的好与坏（就是模型对数据的拟合情况），这样做的好处是即使数据集无法用图像的方式可视化展现，我们也可以用评价指标来评价它。
python的库（如sklearn）可以方便地实践机器学习，而且具有很大的灵活性（相应地，难度也上升了）。我们可以使用python自己编写自己的机器学习算法，说明机器学习并不是什么神秘的东西。从实际编程中可以发现，使用python来做机器学习，需要掌握一定的编程思想。