1. 引言
這一章的主題是線性回歸。我們按照這一順序來展開教程：
- 什麼是線性回歸
- 如何用sklearn來實作線性回歸的計算
- 線性回歸的算法原理是什麼
- 如何用python建立線性回歸的算法
- 總結
2. 什麼是線性回歸
線性回歸回歸可以說是機器學習的入門算法，類似于學代碼先學會如何列印“Hello World”，是以我們從實作線性回歸入手。
先來說說線性回歸的上延概念：機器學習。機器學習分為兩種：一種是回歸（regression）學習，一種是分類（classification）學習。為什麼這裡叫做“回歸學習”和“分類學習”呢？因為機器學習算法的重點在于對資料集的學習，而不同的機器學習算法，隻是實作這種學習的一種載體，是以機器學習與以往的統計分析有所同而有所不同。
（接着回來講回歸和分類。）
-
回歸就是對一個資料集的趨勢進行拟合，進而産生一個模型。
例如，我們手裡有一段這樣的資料集：
藍色點代表資料集，它是一個二維點集，每個點代表一個樣本（sample）。紅色曲線代表線性回歸拟合曲線。線性回歸分析的結果就是創造這樣一條拟合曲線，它反映了這個樣本點集的整體”趨勢”。
Python機器學習與資料分析系列(2)-線性回歸
可以想象，當我們擁有了紅色“拟合曲線”以後，就可以用于“預測”。例如，橫坐标X軸（x-asix）給定一個數，那麼就可以根據拟合曲線，得到縱坐标Y軸（y-axis）的值、
如果我們将X軸上的數當成自變量，Y軸上的數當成因變量，那麼就相當于建立了從自變量-因變量的數學關系，上面的“預測”就可以翻譯為：給定你一個樣本的自變量的值，根據拟合曲線，就可以給出一個樣本的因變量的值，進而實作對這個樣本的“預測”。
-
分類就是對一個資料集中不同種類樣本的分類，也産生一個模型。
例如，我們手裡有這樣一段資料：
黑色的點代表一類樣本，紅色的點代表另外一類樣本，機器學習分類模型學習得到黑色、紅色兩類樣本的分類模型。然後我們就可以将一個新的樣本點（藍色點），對它的歸類進行“預測”。
Python機器學習與資料分析系列(2)-線性回歸
由此可見，機器學習也是一門用于“預測”模組化的學問，實作其“預測”能力的重點在于“學習”，而“學習”是基于資料集的。在後面的章節中，你會從代碼層面了解這種學習機制。
然後我們也從表面上了解了什麼是“線性回歸”，這一章主要圍繞它展開。
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]
注意，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. 線性回歸的原理和實作
前面兩節我們了解了線性回歸的使用方法，現在來看一看它的原理。
假如我們面對一個二維點集組成的資料：
黃色曲線是線性回歸曲線，代表了資料集的“趨勢”，這很直覺。
不過，當資料集是這樣的時候：
情況就比較複雜了。我們無法直覺地劃出這個“回歸曲線”，也沒有辦法定義什麼樣的“回歸曲線”的好與壞。
是以仔細回來看看這條曲線（實際上就是直線），用初等數學的知識可以知道，它可以用一個函數來表達：
y=mx+b
即得到一個函數關系：
y=f(x|m,b)
在這裡，x就是資料集的“features”，也是自變量。y就是資料集的“label”，也是因變量，y也是我們需要預測的結果。而 m,b 是這個函數的參數，即組成模型的參數（在這個直線方程中，m就是斜率，b就是截距），帶有幾何含義。但這裡features是一維的，如果features有很多唯，也就很難從幾何上了解w和b了。
如何得到最能代表上述資料集的直線？如果feature是一維的，label也是一維的話，根據推導可以得到：
m=x⎯⎯⋅y⎯⎯−xy⎯⎯⎯⎯(x⎯⎯)2−x2⎯⎯⎯⎯
得到了m的公式，最佳截距的公式可以反推得到：
b=y⎯⎯−mx⎯⎯
在這個過程我們省去了繁瑣的數學證明過程（這篇文章隻是為了熟練機器學習的代碼部分，而非數學理論部分）。
現在我們用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，方差）
SE=∑(yˆ−y)2
顯然SE就是統計每個樣本點的y值(y_orig)與回歸曲線的縱坐标(ys_line)內插補點的求和（也可以求平均，效果都差不多）：
- R2 ，也稱為确定系數（coefficient of determination)
R2=1−SEyˆSEy⎯⎯
這裡的 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原點與曲線取點的差的平方。
- 均線 y⎯⎯ 是一條直線，也是一個由相同點組成的點集（這裡同樣使用了 [forin(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|θ)
結合股票預測的例子可以發現，影響股票收盤價的features就是這裡的 x ，也就是自變量。收盤價（預測的目标），就是這裡的y。機器學習會利用已有的sample組成的資料集dataset，機器學習通過對資料集的學習（這篇文章沒有展現出來，因為我們直接得出了線性回歸的模型參數求值方法），來得到一個良好的 θ ，進而建立起這個映射函數模型。對于一個新的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來做機器學習，需要掌握一定的程式設計思想。