阿曼多·凡丹戈（Armando Fandango）的《python資料分析（第2版）》，2018年6月出版，本系列為讀書筆記。主要是為了系統整理，加深記憶。第10章“預測性分析與機器學習”通過一個例子來說明人工智能在天氣預報上的應用，這主要借助于scikit-learn。不過，有些機器學習算法在scikit-learn中尚未實作，是以有時還要求助其他API。
python資料分析個人學習讀書筆記-目錄索引
第10章預測性分析與機器學習
最近，預測性分析與機器學習已經納入許多行業的主流資料科學和資料分析的行列。相對于其他領域而言，我們對這兩個領域取得突飛猛進發展的期待要更熱切一些。甚至有人預言，機器學習的發展速度将日益加快，是以，幾十年内，人工就會被智能機器所替代。當然，就目前人工通用智能（AGI）的發展現狀來看，這隻是一個遙遠的烏托邦而已；但是機器學習已經取得了長足的進步，可用于自駕車、聊天機器人和AI助手，如亞馬遜的Alexa、蘋果的Siri以及Ok Google。然而，目前即使是進行非常簡單的判斷，如判斷網絡圖檔中是否含有貓或狗等，都需要大量的運算和資料作為支撐。預測性分析則需要借助各種各樣的技術，包括機器學習，才能做出有用的判斷。例如，某客戶是否有能力償還其貸款，或者某位女性客戶是否有孕在身。
為了完成這些預測，需要從海量資料中提取特征。關于特征，我們之前也曾經提到過，它們又被稱為預測變量。進行預測時，特征通常用于輸入變量。實質上，特征可以從資料中提取，然後要做的是，找到一個函數，将特征映射到目标上。當然，這個目标可能是已知的，也可能是未知的。尋找合适的函數很難，為此，通常需要把多種不同的算法和模型組合在一起，也就是所謂的內建。內建的輸出結果可以是一組模型投票決出的結果，也可以是所有結果的一個折中。但是，我們還可以使用另外一種更加進階的算法來獲得最終結果。雖然我們不會在本章使用內建技術，但是大家還是有必要記住這種技術。
實際上，在前面的章節中我們已經接觸過機器學習算法了，樸素貝葉斯分類算法便是其中之一。我們可以将機器學習分為下列幾種類型。
- 監督學習：要求為訓練資料提供标簽，也就是說需要給算法提供已經分好類的樣本。利用帶标簽的訓練資料，我們可以建立一個函數，将輸入變量映射為相應的輸出變量。比如，如果想對垃圾郵件進行分類，那麼我們必須提供垃圾郵件和正常電子郵件的相應樣本。監督學習算法的例子包括線性回歸、邏輯回歸、狀态向量機、随機森林、K最近鄰算法等。
- 無監督學習：這類機器學習算法不需要人工輸入。它能夠自行發現資料中存在的模式，如大型資料集中的聚類等。無監督學習算法的例子包括K均值和分層聚類。
- 強化學習：這種類型的學習技術無需進行輔導，但是需要我們提供一些回報資訊。例如一台電腦可以跟它自己下棋，或者玩1983的《戰争遊戲》電影中的井字遊戲和熱核戰争遊戲。
下面我們以天氣預報為例進行說明。本章會大量用到一個Python程式庫，即scikit- learn。雖然這個程式庫提供了許多聚類分析算法、回歸算法和分類算法，但是，還有一些機器學習算法未包含在scikit-learn中，是以，我們還需要用到其他一些API。本章涉及的主題如下。
- 預處理
- 基于邏輯回歸的分類
- 基于支援向量機的分類
- 基于ElasticNetCV的回歸分析
- 支援向量回歸
- 基于相似性傳播（affinity propagation）的聚類分析
- 均值漂移算法
- 遺傳算法
- 神經網絡
- 決策樹算法
10.1預處理
在scikit-learn庫中提供了一個預處理子產品，下面我們将詳細介紹該子產品。在第9章中，我們已經安裝了scikit-learn并且已經做過一次資料預處理了，即過濾掉剔除字（Stopwords）。一些機器學習算法對某些資料比較頭疼，因為這些資料不服從高斯分布，即不滿足數學期望為0、标準方差為1的條件。子產品sklearn.preprocessing進而應運而生，本節将詳細介紹這個子產品的使用方法，我們會針對來自荷蘭皇家氣象學會的氣象資料進行預處理。比爾特（DE BILT）氣象研究中心的原始資料較多，我們所要的資料隻是原始資料檔案中的一列而已，這一列記錄的是日降雨量。需要的資料将使用第5章中介紹的.npy格式來存放。下面把資料加載到一個NumPy數組中。這些整數值必須全部乘以0.1，進而得到以毫米為機關的日降水量。
該資料有一個很詭異的地方，那就是凡是小于 0.05mm的數值，都将用作−1。是以，我們把這些數值全部設為0.025，即0.05的一半。在原始資料中，如果某一天的降雨量低于這些數值，那麼這一天将不會記錄在案。對于這樣的缺失資料，我們将完全忽略。我們可以那麼做，因為本來缺失的資料點就非常多，是以也不差這一點。例如，本世紀初大約缺失一年的資料，同時後來許多天的資料也是缺失的。在preprocessing子產品中，它的Imputer類提供了許多處理缺失資料的預設政策。不過，這些政策對于本節讨論的内容而言，好像不太合适。資料分析如同一扇透視資料的視窗，一扇通向知識的視窗。資料的清洗和填補可以使我們的視窗看起來更加清爽。可是，我們應該注意不要過于扭曲原始資料。
這裡的機器學習樣本的主要特征是一個數組，用來存放表示一年中各天（1～366）日期的數值。這些可以幫助解釋季節影響因素。
期望值、标準差和安德森-達林（Anderson-Darling）檢驗結果（詳情請參閱第4章）如下
Rain mean 2.17919594267
Rain variance 18.803443919
Anderson rain (inf, array([ 0.576, 0.656, 0.787, 0.918,1.092]), array([ 15. , 10. , 5. , 2.5, 1. ])) 透過以上輸出内容，可以得出可靠的結論，該資料不滿足數學期望為0标準方差為1的要求，是以它不符合正态分布。在這個資料集中，0值占比很大，說明在這些天沒有下雨。大的降雨量已經越來越罕見（這倒是一件好事）。資料分布情況完全不對稱，是以，這不是一個高斯分布。下面通過scale()函數對資料進行縮放處理，代碼如下。
1 scaled = preprocessing.scale(rain) 如今，數學期望和标準方差已經滿足要求，隻是資料的分布還是不對稱。
Scaled mean 3.4130160280768244e-17
Scaled variance 1.0
Anderson scaled AndersonResult(statistic=5731.267850009339, critical_values=array([0.576, 0.656, 0.787, 0.918, 1.092]),
significance_level=array([15. , 10. , 5. , 2.5, 1. ])) 有時，我們需要把特征值由數值型轉換為布爾型。例如，進行文本分析時就經常需要進行此類變換，以簡化計算。這時，可以借助binarize()函數進行類型轉換，代碼如下。
1 binarized = preprocessing.binarize(rain.reshape(-1,1))
2 print(np.unique(binarized), binarized.sum()) 預設情況下，這将生成一個新數組；此外，我們還讓它進行了某些運算。預設的門檻值為0。也就是說，正值用1替代，負值用0替代。
[0. 1.] 24594.0 進行分類時，類LabelBinarizer可以用整數來标注類别。
1 lb = preprocessing.LabelBinarizer()
2 lb.fit(rain.astype(int))
3 print(lb.classes_) 上述代碼的輸出0～62的一組整數。相關代碼摘自本書代碼包中的ch-10.ipynb檔案。
1 import numpy as np
2 from sklearn import preprocessing
3 from scipy.stats import anderson
4
5 rain = np.load('rain.npy')
6 rain = .1 * rain
7 rain[rain < 0] = .05/2
8 print("Rain mean", rain.mean())
9 print("Rain variance", rain.var())
10 print("Anderson rain", anderson(rain))
11
12 scaled = preprocessing.scale(rain)
13 print("Scaled mean", scaled.mean())
14 print("Scaled variance", scaled.var())
15 print("Anderson scaled", anderson(scaled))
16
17 binarized = preprocessing.binarize(rain.reshape(-1,1))
18 print(np.unique(binarized), binarized.sum())
19
20 lb = preprocessing.LabelBinarizer()
21 lb.fit(rain.astype(int))
22 print(lb.classes_) 10.2基于邏輯回歸的分類
邏輯回歸是一種分類算法，這種算法可以用于預測事件發生的機率或者某事物屬于某一類别的機率。對于多元分類問題，可以将其簡化為二進制分類問題。在最簡單的情形下，某一類别的機率高，則意味着另一個類别的機率低。邏輯回歸是以logistic函數為基礎的，該函數的取值介于0～1中間，這與機率值正好吻合。是以，可以利用logistic函數将任意值轉換為一個機率。
為了使用邏輯回歸進行分類，需要定義一個相應的函數。下面開始建立分類器對象，代碼如下。
1 clf = LogisticRegression(random_state=12) 參數random_state的作用，類似為随機數生成器指定的種子。之前已經講過，交叉驗證在避免過拟合方面有着非常重要的作用。K-折交叉驗證是一種交叉驗證技術，它會把資料集随機分為k（一個小整數）份，每一份稱為一個包。在這k次疊代過程中，每個包會有1次被用于驗證，其餘9次用于訓練。對于scikit-learn來說，它每個類的預設k值都是3，但是，通常需要将這個值設定得更大一些，如5或10。疊代的結果可以在最後進行合并。對于k-折交叉驗證，scikit-learn專門提供了一個KFold類。為了建立一個具有10個包的KFold對象，可以使用下列代碼。
1 kf = KFold(n_splits=10) 然後，使用fit()函數訓練資料，代碼如下。
1 clf.fit(x[train], y[train]) 為了衡量分類的準确性，可以使用score()方法，具體如下。
1 scores.append(clf.score(x[test], y[test])) 在這個例子中，可以把日期和日降雨量作為特征，下面使用這些特征來構造一個數組，代碼如下。
1 x = np.vstack((dates[:-1], rain[:-1])) 因為要進行分類，是以首先要定義無降雨天，即降雨量為0；然後用−1表示有微弱降雨的小雨天；最後剩下的就是雨天。我們把這3種類别與資料值的符号關聯起來。
1 y = np.sign(rain[1:]) 據此，得到的平均準确率為57%。在scikit-learn的資料集上，我們的平均準确率為41%，代碼具體參考本書代碼包中的log_regress.py檔案。
1 from sklearn.linear_model import LogisticRegression
2 from sklearn.model_selection import KFold
3 from sklearn import datasets
4 import numpy as np
5
6 def classify(x, y):
7 clf = LogisticRegression(random_state=12)
8 scores = []
9 kf = KFold(n_splits=10)
10 for train,test in kf.split(x):
11 clf.fit(x[train], y[train])
12 scores.append(clf.score(x[test], y[test]))
13
14 print("Accuracy: ",np.mean(scores))
15
16 rain = np.load('rain.npy')
17 dates = np.load('doy.npy')
18
19 x = np.vstack((dates[:-1], rain[:-1]))
20 y = np.sign(rain[1:])
21 classify(x.T, y)
22
23 #iris example
24 iris = datasets.load_iris()
25 x = iris.data[:, :2]
26 y = iris.target
27 classify(x, y) Accuracy: 0.577211706428828
Accuracy: 0.7066666666666667 10.3基于支援向量機的分類
支援向量機（Support Vector Machines，SVM）不僅可以用來進行回歸分析，如支援向量回歸（Support Vector Regression，SVR），而且還可以用來進行分類，支援向量機算法是Vladimir Vapnik于1993年發明出來的。
SVM可以把資料點映射到多元空間的資料點，這種映射是通過核函數來完成的。這裡的核函數既可以是線性的，也可以是非線性的。這樣，分類問題就簡化為尋找一個将空間一分為二的超平面，或者是能夠将資料點恰如其分地劃分到不同空間（類别）的多個超平面。利用超平面進行分類是一件非常困難的事情，因為這會引出軟間隔的概念。軟間隔用來表示對錯誤分類的容忍度，通常由一個用C表示的常量給出。另一個重要的參數就是核函數的類型，有如下幾種。
- 線性函數
- 多項式函數
- 徑向基函數
- Sigmoid函數
網格搜尋法可以用來尋找合适的參數，這是一種系統性的方法，即嘗試所有可能的參數組合。為了進行網格搜尋，可以借助scikit-learn提供的GridSearchCV類。使用這個類時，我們可以通過字典來提供分類器或回歸器的類型對象。字典的鍵就是我們将要調整的參數，而字典的值就是需要嘗試的參數值的相應清單。對于用來進行分類和回歸分析的那些類，scikit-learn API都有對應的、為其提供了添加交叉驗證功能的類。不過，預設情況下，交叉驗證功能都未啟用。下面建立一個GridSearchCV對象，代碼如下。
1 clf = GridSearchCV(SVC(random_state=42, max_iter=100), {'kernel': ['linear', 'poly', 'rbf'], 'C':[1, 10]}) 上面的代碼規定了最多疊代次數，這個數字不要太大，否則對我們的耐心絕對是一個嚴峻的考驗。此外，這裡沒有啟用交叉驗證，為的是加速處理過程。而且，我們還指定了核方法的類型以及軟間隔的相應參數。
上述代碼為可能的參數變異建立了一個2×3的網格。如果時間寬裕，我們可以建立一個更大的網格來包含更多的可能值。另外，我們可以把GridSearchCV的cv參數的值設為預期的包總量，如5或10。實際上，最多疊代次數設定得越高越好。至于不同的核方法，可以根據拟合的實際情況擇機而用。如果把參數verbose設定為一個非零整數，那就可以得到每個參數值組合更為詳盡的輸出資訊，如執行時間等。通常情況下，我們希望按照幅值大小的順序，如1～10 000，依次變換 soft-margin 參數的取值，這時可以借助 NumPy 的logspace()函數來完成。
使用這個分類器時，在天氣資料上得到的準确性是56%，在鸢尾花樣本資料集上得到的準确性是82%。GridSearchCV的gridscores區儲存有網格搜尋的成效資料。對于天氣資料，成效資料如下。
Accuracy: 0.4752925545505647
{'mean_fit_time': array([0.10729136, 0.15199084, 0.27729039, 0.07938633, 0.15557127,
0.26728883]),
'mean_score_time': array([0.02516398, 0.03653703, 0.07361264, 0.01555986, 0.03631983,
0.070818 ]),
'mean_test_score': array([0.42610763, 0.53735336, 0.28585569, 0.29270417, 0.53349511,
0.28585569]),
'param_C': masked_array(data=[1, 1, 1, 10, 10, 10],
mask=[False, False, False, False, False, False],
fill_value='?',
dtype=object),
'param_kernel': masked_array(data=['linear', 'poly', 'rbf', 'linear', 'poly', 'rbf'],
mask=[False, False, False, False, False, False],
fill_value='?',
dtype=object),
'params': [{'C': 1, 'kernel': 'linear'},
{'C': 1, 'kernel': 'poly'},
{'C': 1, 'kernel': 'rbf'},
{'C': 10, 'kernel': 'linear'},
{'C': 10, 'kernel': 'poly'},
{'C': 10, 'kernel': 'rbf'}],
'rank_test_score': array([3, 1, 5, 4, 2, 5]),
'split0_test_score': array([0.40252938, 0.55952989, 0.23441492, 0.5167348 , 0.55556975,
0.23441492]),
'split1_test_score': array([0.31732243, 0.55838017, 0.38924374, 0.21525294, 0.55620848,
0.38924374]),
'split2_test_score': array([0.55799693, 0.55825243, 0.41377108, 0.29036791, 0.55697496,
0.41377108]),
'split3_test_score': array([0.55960138, 0.45189728, 0.16877475, 0.30586432, 0.45189728,
0.16877475]),
'split4_test_score': array([0.29308803, 0.55870704, 0.22307397, 0.13530088, 0.54682509,
0.22307397]),
'std_fit_time': array([0.00754749, 0.00337823, 0.01541952, 0.00506452, 0.00372173,
0.01022942]),
'std_score_time': array([0.00146306, 0.00118295, 0.00172066, 0.00135341, 0.00120141,
0.00209969]),
'std_test_score': array([0.11428013, 0.04273036, 0.09731092, 0.12739141, 0.04096434,
0.09731092])} View Code
對于鸢尾花樣本資料來說，我們的成績如下。
Accuracy: 0.82
{'mean_fit_time': array([0.0009974 , 0.00099716, 0.00099726, 0.00079803, 0.0009974 ,
0.00119658]),
'mean_score_time': array([0.00039902, 0. , 0.00019932, 0.00039868, 0.00019941,
0.00019956]),
'mean_test_score': array([0.80666667, 0.73333333, 0.82 , 0.76666667, 0.59333333,
0.81333333]),
'param_C': masked_array(data=[1, 1, 1, 10, 10, 10],
mask=[False, False, False, False, False, False],
fill_value='?',
dtype=object),
'param_kernel': masked_array(data=['linear', 'poly', 'rbf', 'linear', 'poly', 'rbf'],
mask=[False, False, False, False, False, False],
fill_value='?',
dtype=object),
'params': [{'C': 1, 'kernel': 'linear'},
{'C': 1, 'kernel': 'poly'},
{'C': 1, 'kernel': 'rbf'},
{'C': 10, 'kernel': 'linear'},
{'C': 10, 'kernel': 'poly'},
{'C': 10, 'kernel': 'rbf'}],
'rank_test_score': array([3, 5, 1, 4, 6, 2]),
'split0_test_score': array([0.73333333, 0.7 , 0.73333333, 0.73333333, 0.56666667,
0.73333333]),
'split1_test_score': array([0.83333333, 0.86666667, 0.86666667, 0.83333333, 0.5 ,
0.83333333]),
'split2_test_score': array([0.76666667, 0.66666667, 0.76666667, 0.76666667, 0.6 ,
0.8 ]),
'split3_test_score': array([0.86666667, 0.56666667, 0.86666667, 0.86666667, 0.9 ,
0.83333333]),
'split4_test_score': array([0.83333333, 0.86666667, 0.86666667, 0.63333333, 0.4 ,
0.86666667]),
'std_fit_time': array([3.23406696e-07, 4.67203091e-07, 2.78041453e-07, 3.99017505e-04,
1.90734863e-07, 3.98779158e-04]),
'std_score_time': array([0.00048869, 0. , 0.00039864, 0.00048829, 0.00039883,
0.00039911]),
'std_test_score': array([0.04898979, 0.11737878, 0.05811865, 0.08164966, 0.16786238,
0.04521553])} 相關代碼取自本書代碼包中的ch-10.ipynb檔案。
1 from sklearn.svm import SVC
2 from sklearn.model_selection import GridSearchCV
3 from sklearn import datasets
4 from sklearn.preprocessing import StandardScaler
5 import numpy as np
6 from pprint import PrettyPrinter
7
8 def classify(x, y):
9 clf = GridSearchCV(SVC(random_state=42, max_iter=100), {'kernel': ['linear', 'poly', 'rbf'], 'C':[1, 10]})
10 clf.fit(x, y)
11 print("Accuracy: ", clf.score(x, y))
12 PrettyPrinter().pprint(clf.cv_results_)
13
14 rain = np.load('rain.npy')
15 dates = np.load('doy.npy')
16
17 x = np.vstack((dates[:-1], rain[:-1]))
18 y = np.sign(rain[1:])
19
20 classify(x.T, y)
21
22 #iris example
23 iris = datasets.load_iris()
24 x = iris.data[:, :2]
25 y = iris.target
26 classify(x, y) 10.4基于ElasticNetCV的回歸分析
彈性網絡正則化（Elastic Net Regularization）是一種降低回歸分析的過拟合風險的方法。彈性網絡正則化實際上就是LASSO（The Least Absolute Shrinkage and Selection Operator，LASSO）算法和嶺回歸方法的線性組合。LASSO能夠有效限制L1範數或曼哈頓距離。對于兩點來說，L1範數表示的是它們坐标值之差的絕對值之和。嶺回歸方法使用L1範數的平方作為懲罰項。處理回歸問題時，拟合優度通常是由所謂的R平方的判定系數來判斷的。令人遺憾的是，R平方的定義并不統一。此外，這個名稱還有誤導之嫌，因為實際上它可以為負數。就拟合優度而言，最好的情況下其判定系數為1。根據判定系數的定義看，其允許的取值範圍還是不小的，但是我們的目标很明确，就是設法讓它向1靠攏。
下面我們使用10-折交叉驗證并定義一個ElasticNetCV對象，代碼如下。
clf = ElasticNetCV(max_iter=200, cv=10, l1_ratio = [.1, .5, .7, .9, .95, .99, 1]) 這裡，ElasticNetCV類使用了一個名為l1_ratio的參數，其取值為0～1。如果該參數的值為0，隻使用嶺回歸算法；如果該參數的值為1，隻使用LASSO回歸算法；否則，就使用混合算法。對于這個參數，可以給它指定單個數值，也可以給它提供一個數值清單。對于降雨資料，我們得分情況如下。
Score 0.052783876094198205 這個得分表明，我們對資料訓練不足，或者說産生欠拟合（underftting）。出現這種情況的原因有很多，如特征數量不足或者選用的模型不合适等。對于波士頓房價資料，針對現有的特征得分如下。
Score 0.6831524128217573 Predict()方法可以針對新資料進行預測。對于預測的結果，我們可以使用散點圖來可視化其效果。對于降雨資料，得到的散點圖如圖10-1所示。
觀察圖10-1可以發現，這次的拟合效果不佳，即欠拟合。穿過原點的對角線，才是最理想的拟合曲線，就像圖10-2所示的波士頓房價資料拟合情況。如上，兩圖合并
相關代碼摘自本書代碼包中的ch-10.ipynb檔案。
1 from sklearn.linear_model import ElasticNetCV
2 import numpy as np
3 from sklearn import datasets
4 import matplotlib.pyplot as plt
5
6 def regress(x, y, title):
7 clf = ElasticNetCV(max_iter=200, cv=10, l1_ratio = [.1, .5,
8 .7, .9, .95, .99, 1])
9
10 clf.fit(x, y)
11 print("Score", clf.score(x, y))
12
13 pred = clf.predict(x)
14 plt.title("Scatter plot of prediction and " + title)
15 plt.xlabel("Prediction")
16 plt.ylabel("Target")
17 plt.scatter(y, pred)
18 # Show perfect fit line
19 if "Boston" in title:
20 plt.plot(y, y, label="Perfect Fit")
21 plt.legend()
22
23 plt.grid(True)
24 plt.show()
25
26 rain = .1 * np.load('rain.npy')
27 rain[rain < 0] = .05/2
28 dates = np.load('doy.npy')
29
30 x = np.vstack((dates[:-1], rain[:-1]))
31 y = rain[1:]
32 regress(x.T, y, "rain data")
33
34 boston = datasets.load_boston()
35 x = boston.data
36 y = boston.target
37 regress(x, y, "Boston house prices") 10.5支援向量回歸
如前所述，支援向量機也可以用于回歸分析。就回歸來說，我們是通過超平面而非單獨的點來拟合資料的。學習曲線是一種将學習算法行為特點可視化的好辦法。對于不同的訓練資料量來說，學習曲線都是訓練成效和測試成效所考量的一部分。建立學習曲線可能迫使我們重複訓練估計器，是以會導緻整體過程變慢。對此，我們可以通過使用多個并行估計器作業進行補償。
支援向量回歸是需要進行标定處理（scaling）的算法之一。
我們得到的兩個高分如下。
Max test score Rain 0.0161004084576
Max test score Boston 0.662188537037 這與ElasticNetCV類得到的結果相仿。為此，許多scikit-learn類都提供了一個n_jobs參數。根據經驗，通常系統中有幾個CPU，我們就建立幾個作業。建立作業時，可以借助Python标準的多程序API。為了進行訓練和測試，可以調用learning_curve()函數，代碼 如下。
1 train_sizes, train_scores, test_scores = learning_curve(clf, X, Y, n_jobs=ncpus) 求平均數，然後畫出相應的得分情況。
1 plt.plot(train_sizes, train_scores.mean(axis=1), label="Train score")
2 plt.plot(train_sizes, test_scores.mean(axis=1), '--', label="Test score") 降雨資料的學習曲線大緻如圖10-3所示。
學習曲線在日常生活中也較常見，即練習的次數越多，學到的就越多。按照資料分析術語來講，資料越多，學習的效果越好。如果訓練成績不錯，但是測試效果不好，這就說明出現了過拟合現象。也就是說，我們的模型隻适用于訓練資料。圖10-4所示是波士頓房價資料的學習曲線，它看起來更好一些。兩圖已合并。
代碼取自本書代碼包中的sv_regress.py檔案。
1 import numpy as np
2 from sklearn import datasets
3 from sklearn.model_selection import learning_curve
4 from sklearn.svm import SVR
5 from sklearn import preprocessing
6 import multiprocessing
7 import matplotlib.pyplot as plt
8
9 def regress(x, y, ncpus, title):
10 X = preprocessing.scale(x)
11 Y = preprocessing.scale(y)
12 clf = SVR(max_iter=ncpus * 200)
13
14 train_sizes, train_scores, test_scores = learning_curve(clf, X, Y, n_jobs=ncpus)
15
16 plt.figure()
17 plt.title(title)
18 plt.plot(train_sizes, train_scores.mean(axis=1), label="Train score")
19 plt.plot(train_sizes, test_scores.mean(axis=1), '--', label="Test score")
20 print("Max test score " + title, test_scores.max())
21 plt.grid(True)
22 plt.legend(loc='best')
23 plt.show()
24
25 rain = .1 * np.load('rain.npy')
26 rain[rain < 0] = .05/2
27 dates = np.load('doy.npy')
28
29 x = np.vstack((dates[:-1], rain[:-1]))
30 y = rain[1:]
31 ncpus = multiprocessing.cpu_count()
32 regress(x.T, y, ncpus, "Rain")
33
34 boston = datasets.load_boston()
35 x = boston.data
36 y = boston.target
37 regress(x, y, ncpus, "Boston") 10.6基于相似性傳播算法的聚類分析
聚類分析，就是把資料分成一些組，這些組即所謂的聚類。聚類分析通常無需提供目标資料，從這個意義上來說，它屬于無監督學習方法。一些聚類算法需要對聚類數進行推測，而另一些聚類算法則不必如此，相似性傳播（Affinity Propagation，AP）算法就屬于後一類算法。資料集中的各個元素都會通過特征值被映射到歐氏空間，然後計算資料點之間的歐氏距離，以此建構一個矩陣，這個矩陣就是AP算法的基礎。因為這個矩陣可能會迅速膨脹，對此我們一定要當心，不要讓它耗盡記憶體。實際上，scikit-learn程式庫提供了一個實用程式，來幫我們生成結構化資料。下面我們通過代碼生成3個資料塊，具體如下。
1 x, _ = datasets.make_blobs(n_samples=100, centers=3, n_features=2,
2 random_state=10) 調用euclidean_distances()函數，建立前面提到的矩陣。
1 S = euclidean_distances(x) 利用這個矩陣，我們就可以給資料标注其所屬的聚類了，代碼如下。
1 aff_pro = cluster.AffinityPropagation().fit(S)
2 labels = aff_pro.labels_ 繪制聚類後，得到的圖像如圖10-5所示。
以上代碼摘自本書代碼包中的ch-10.ipynb.py檔案。
1 from sklearn import datasets
2 from sklearn import cluster
3 import numpy as np
4 import matplotlib.pyplot as plt
5 from sklearn.metrics import euclidean_distances
6
7 x, _ = datasets.make_blobs(n_samples=100, centers=3, n_features=2,
8 random_state=10)
9 S = euclidean_distances(x)
10
11 aff_pro = cluster.AffinityPropagation().fit(S)
12 labels = aff_pro.labels_
13
14 styles = ['o', 'x', '^']
15
16 for style, label in zip(styles, np.unique(labels)):
17 print(label)
18 plt.plot(x[labels == label], style, label=label)
19 plt.title("Clustering Blobs")
20 plt.grid(True)
21 plt.legend(loc='best')
22 plt.show() 10.7均值漂移算法
均值漂移算法是另外一種不需要估算聚類數的聚類算法。目前，這個算法已經成功應用于圖像處理。該算法通過疊代來尋找一個密度函數的最大值。展示均值漂移算法前，我們首先需要使用pandas DataFrame計算日降雨量的平均值。下面來建立DataFrame并計算其資料的平均值，代碼如下。
1 df = pd.DataFrame.from_records(x.T, columns=['dates', 'rain'])
2 df = df.groupby('dates').mean()
3 df.plot() 結果如圖10-6所示。
基于均值漂移算法的聚類代碼如下所示。
1 x = np.vstack((np.arange(1, len(df) + 1) , df.values.ravel()))
2 x = x.T
3 ms = cluster.MeanShift()
4 ms.fit(x) 如果使用不同的線寬和陰影，那麼繪制出的3個聚類如圖10-7所示。
兩圖合并如上。
如你所見，我們根據年内日平均雨量（以mm為機關）劃分出3個聚類。完整的代碼摘自本書代碼包中的ch-10.ipynb檔案。
1 import numpy as np
2 from sklearn import cluster
3 import matplotlib.pyplot as plt
4 import pandas as pd
5
6 rain = .1 * np.load('rain.npy')
7 rain[rain < 0] = .05/2
8 dates = np.load('doy.npy')
9 x = np.vstack((dates, rain))
10 df = pd.DataFrame.from_records(x.T, columns=['dates', 'rain'])
11 df = df.groupby('dates').mean()
12 df.plot()
13 # x = np.vstack((np.arange(1, len(df) + 1) , df.as_matrix().ravel()))
14 #'DataFrame' object has no attribute 'as_matrix'，改為df.values即可，邀月注
15 x = np.vstack((np.arange(1, len(df) + 1) , df.values.ravel()))
16 x = x.T
17 ms = cluster.MeanShift()
18 ms.fit(x)
19 labels = ms.predict(x)
20
21 plt.figure()
22 grays = ['0', '0.5', '0.75']
23
24 for gray, label in zip(grays, np.unique(labels)):
25 match = labels == label
26 x0 = x[:, 0]
27 x1 = x[:, 1]
28 plt.plot(x0[match], x1[match], lw=label+1, label=label)
29 plt.fill_between(x0, x1, where=match, color=gray)
30
31 plt.grid(True)
32 plt.legend()
33 plt.show() 10.8遺傳算法
遺傳算法是本書最具争議的部分，這種算法基于生物學領域的進化論。這種類型的算法在搜尋和優化方面用途廣泛，例如，可以使用遺傳算法來搜尋回歸問題或分類問題的最佳參數。
在地球上，人類及其他生命形式都通過染色體來攜帶遺傳資訊。通常，我們使用字元串對染色體進行模組化。在遺傳算法中，我們也沿用了這種類似的遺傳資訊表示方法。算法的第一步是，利用随機個體初始化種群并使用相應的表示方法表達遺傳資訊。此外，我們還可以使用待解決問題的已知候選解來進行初始化。之後，我們需要進入一個疊代處理過程，即所謂的逐代（generation）演化。對于每一代個體，算法都會根據預定義的适應度函數來選擇一些個體，作為交配個體供下一步雜交用。适應度函數的作用是評估個體與滿意解的接近程度。
下面介紹兩種遺傳算子（genetic operators），它們可以用來生成新的遺傳資訊。
- 交叉（Crossover）：這個遺傳操作就是通過交配的方式來産生新個體。這裡将介紹單點交叉（one-point crossover）操作，具體操作是，親本的遺傳資訊互相交換，進而産生兩個新的個體，第1個個體前半段遺傳資訊取自父本，後半段取自母本，第2個個體正好相反。例如，假設遺傳資訊使用100個清單元素來表示，那麼交叉操作可以從親本中的一方摘取前80個元素，然後從另一方取得後20個元素。當然，對于遺傳算法來說，它能夠通過交叉兩個以上的父本和母本來繁殖新個體，這一技術正處于研究階段（詳情參見Eiben, A. E. et al. Genetic algorithms with multi-parent recombination, Proceedings of the International Conference on Evolutionary Computation- PPSN III. The Third Conference on Parallel Problem Solving from Nature: 78-87. ISBN 3-540-58484-6, 1994）。
- 突變（Mutation）：這個遺傳算子受固定突變率的限制。這個概念在好萊塢電影和大衆文化中屢見不鮮。我們知道，突變雖然極為罕見，但是卻能帶來緻命危害。不過，有時候突變也能帶來我們夢寐以求的特性，而且在某些情況下，這些特性還能夠遺傳到下一代。
最終，新個體會取代舊個體，這時，我們就可以進入下一輪疊代了。本例中，我們将使用Python的DEAP程式庫來示範遺傳算法。我們安裝DEAP，具體方法如下。
1 $ pip3 install deap 首先，定義最大化适應度的Fitness子類，具體如下。
1 creator.create("FitnessMax", base.Fitness, weights=(1.0,)) 然後，為種群中的個體定義一個模闆。
1 creator.create("Individual", array.array, typecode='d',
2 fitness=creator.FitnessMax) DEAP中的工具箱用來注冊必需的函數。下面建立一個工具箱并注冊初始化函數，代碼如下。
1 toolbox = base.Toolbox()
2 toolbox.register("attr_float", random.random)
3 toolbox.register("individual", tools.initRepeat,
4 creator.Individual, toolbox.attr_float, 200)
5 toolbox.register("populate", tools.initRepeat, list,
6 toolbox.individual) 第1個函數的作用是産生浮點數，這些數的取值範圍為0～1。第2個函數的作用是生成一個個體，它實際上是一個由200個浮點數組成的清單。第3個函數的作用是建立由個體構成的清單，這個清單代表的是一個種群，也就是搜尋或優化問題的一組可能解。
現實社會中，大部分人都很“普通”，或者說很“正常”，但是也有極少數人很“另類”，如愛因斯坦。第3章曾經介紹了一個shapiro()函數，它可以用來進行正态檢驗。對于一個個體是否為普通個體，我們需要通過其正态檢驗p值來衡量它是普通個體的可能性有多大。下面的代碼将定義一個适應度函數。
1 def eval(individual):
2 return shapiro(individual)[1], 下面定義遺傳算子，代碼如下。
1 toolbox.register("evaluate", eval)
2 toolbox.register("mate", tools.cxTwoPoint)
3 toolbox.register("mutate", tools.mutFlipBit, indpb=0.1)
4 toolbox.register("select", tools.selTournament, tournsize=4) 下面對這些遺傳算子進行簡單解釋。
- Evaluate：即評估算子。這項操作用來度量各個個體的适應度。本例中，正态檢驗的p值便是适應度的衡量名額。
- Mate：即交配算子。這項操作用來産生子代。本例中采用的是兩點交叉。
- Mutate：即突變算子。這項操作會随機修改個體。對于由布爾值構成的清單而言，這就意味着某些值會被反轉，即由True變為False，或者正好相反。
- Mutate：即選擇算子。這項操作用來選擇可以進行交配的個體。
在上面的代碼中，我們規定使用兩點交叉，同時指定了屬性被翻轉的機率，下面生成一個由400個個體組成的原始群體，代碼如下。
1 pop = toolbox.populate(n=400) 現在，啟動進化過程，代碼如下。
1 hof = tools.HallOfFame(1)
2 stats = tools.Statistics(key=lambda ind: ind.fitness.values)
3 stats.register("max", np.max)
4
5 algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=80,
6 stats=stats, halloffame=hof) 這個程式将提供包括每代最大适應度在内的各種統計資訊。我們給它規定了交叉機率、突變率和最大代數（maximum generations），即超過這個代數後就會停止運作。下面的資料就是摘自這個程式輸出的統計報告。
gen nevals max
0 400 0.000484774
1 222 0.000656187
2 246 0.00745961
3 239 0.00745961
4 240 0.0184182
5 216 0.0309736
6 237 0.06957
7 243 0.06957
8 231 0.224381
9 226 0.224381
10 247 0.224381
11 228 0.247313
12 241 0.28318
13 242 0.354144
14 246 0.46282
15 239 0.46282
16 266 0.480937
17 233 0.648529
18 266 0.703452
19 229 0.742139
20 249 0.791523
21 227 0.795744
22 241 0.856505
23 237 0.898701
24 242 0.949595
25 237 0.949595
26 249 0.950638
27 245 0.961397
28 226 0.974227
29 235 0.974227
30 239 0.99192
31 225 0.99192
32 252 0.99192
33 255 0.992183
34 234 0.995053
35 250 0.995053
36 228 0.995053
37 252 0.995466
38 228 0.996523
39 265 0.996992
40 221 0.99771
41 238 0.998339
42 248 0.997815
43 226 0.998156
44 236 0.998442
45 257 0.998442
46 247 0.998485
47 227 0.998564
48 243 0.998564
49 242 0.998692
50 237 0.998692
51 233 0.998692
52 240 0.998692
53 252 0.998692
54 245 0.998713
55 245 0.998713
56 239 0.998713
57 220 0.998713
58 239 0.998782
59 247 0.998838
60 238 0.998838
61 235 0.998838
62 248 0.998861
63 236 0.998861
64 240 0.998861
65 210 0.998861
66 240 0.998861
67 217 0.998861
68 239 0.998861
69 255 0.998861
70 228 0.998861
71 254 0.998861
72 211 0.998861
73 252 0.998861
74 250 0.998861
75 231 0.998861
76 230 0.998861
77 252 0.998861
78 232 0.998861
79 243 0.998861
80 235 0.998861
0.9988605380058289 可見，最初的分布與正态分布相去甚遠，但是，最終得到的個體的分布情況如圖10-8所示。
相關代碼摘自本書代碼包中的ch-10.ipynb檔案：
1 import array
2 import random
3 import numpy as np
4 from deap import algorithms
5 from deap import base
6 from deap import creator
7 from deap import tools
8 from scipy.stats import shapiro
9 import matplotlib.pyplot as plt
10
11 creator.create("FitnessMax", base.Fitness, weights=(1.0,))
12 creator.create("Individual", array.array, typecode='d',
13 fitness=creator.FitnessMax)
14
15 toolbox = base.Toolbox()
16 toolbox.register("attr_float", random.random)
17 toolbox.register("individual", tools.initRepeat,
18 creator.Individual, toolbox.attr_float, 200)
19 toolbox.register("populate", tools.initRepeat, list,
20 toolbox.individual)
21
22 def eval(individual):
23 return shapiro(individual)[1],
24
25 toolbox.register("evaluate", eval)
26 toolbox.register("mate", tools.cxTwoPoint)
27 toolbox.register("mutate", tools.mutFlipBit, indpb=0.1)
28 toolbox.register("select", tools.selTournament, tournsize=4)
29
30 random.seed(42)
31
32 pop = toolbox.populate(n=400)
33 hof = tools.HallOfFame(1)
34 stats = tools.Statistics(key=lambda ind: ind.fitness.values)
35 stats.register("max", np.max)
36
37 algorithms.eaSimple(pop, toolbox, cxpb=0.5, mutpb=0.2, ngen=80,
38 stats=stats, halloffame=hof)
39
40 print(shapiro(hof[0])[1])
41 plt.hist(hof[0])
42 plt.grid(True)
43 plt.show() 10.9神經網絡
人工神經網絡（ANN）的計算模型靈感來自高等動物的大腦。所謂神經網絡，實際上就是由神經元組成的網絡，這些神經元都具有輸入和輸出。例如，可以将圖檔像素相關的數值作為神經網絡的輸入，然後将神經元的輸出傳遞給下一個神經元，依此類推，最終得到一個多層網絡。神經網絡蘊涵了許多自适應元素，這使得它非常适合用于處理非線性模型和模式識别問題。接下來，我們将再次根據前些年對應某日降雨量和前一日降雨量來預測是否會下雨。因為後面要用到Python庫theanets，是以這裡先行安裝，指令如下。
1 $ pip3 install theanets nose_parameterized 一位技術評審在安裝過程中曾經遇到問題，不過更新NumPy和SciPy後，問題就迎刃而解了。首先建立一個Experiment，也就是建立了一個神經網絡，然後再對這個網絡進行訓練。下面的代碼将建立一個具有兩個輸入神經元和一個輸出神經元的神經網絡。
1 net = theanets.Regressor(layers=[2,3,1]) 該網絡擁有一個隐藏層，層内含有3個神經元；另外，這個網絡還使用了Python标準的多程序API來加速運算。訓練網絡所需的資料包括訓練樣本和驗證樣本。
1 train = [x[:N], y[:N]]
2 valid = [x[N:], y[N:]]
3
4 net.train(train,valid,learning_rate=0.1,momentum=0.5) 如果要對驗證資料進行預測，可以使用以下代碼。
1 pred = net.predict(x[N:]).ravel() Scikit-learn程式庫提供了一個實用函數，可以用來計算分類器的準确性。下面是計算準确性的具體代碼。
1 print("Pred Min", pred.min(), "Max", pred.max())
2 print("Y Min", y.min(), "Max", y.max())
3 print("Accuracy", accuracy_score(y[N:], pred >= .5)) 注意，輸出值可能每次都不一樣，這是神經網絡的特性所導緻的。這裡的輸出情況如下。
Pred Min 0.615606596762 Max 0.615606596762
Y Min 0.0 Max 1.0
Accuracy 0.634133878385 相關代碼摘自本書代碼包中的ch-10.ipynb檔案。
1 import numpy as np
2 import theanets
3 import multiprocessing
4 from sklearn import datasets
5 from sklearn.metrics import accuracy_score
6
7 rain = .1 * np.load('rain.npy')
8 rain[rain < 0] = .05/2
9 dates = np.load('doy.npy')
10 x = np.vstack((dates[:-1], np.sign(rain[:-1])))
11 x = x.T
12
13 y = np.vstack(np.sign(rain[1:]),)
14 N = int(.9 * len(x))
15
16 train = [x[:N], y[:N]]
17 valid = [x[N:], y[N:]]
18
19 net = theanets.Regressor(layers=[2,3,1])
20
21 net.train(train,valid,learning_rate=0.1,momentum=0.5)
22
23 pred = net.predict(x[N:]).ravel()
24 print("Pred Min", pred.min(), "Max", pred.max())
25 print("Y Min", y.min(), "Max", y.max())
26 print("Accuracy", accuracy_score(y[N:], pred >= .5)) 所有相關的類都報這個錯，涉及到引用climate，看來這個包有問題。有知道的朋友請告知邀月，謝謝！
d:\tools\python37\lib\site-packages\theanets\layers\base.py in <module>
14 from .. import util
15
---> 16 logging = climate.get_logger(__name__)
17
18 __all__ = [
AttributeError: module 'climate' has no attribute 'get_logger' 10.10決策樹
形如if a: else b這樣的判斷語句，恐怕是Python程式中最常見的語句了。通過嵌套群組合這些語句，就能夠建立所謂的決策樹。決策樹跟老式流程圖非常類似，隻不過流程圖允許循環而已。機器學習領域中，應用決策樹的過程通常被稱為決策樹學習。進行決策樹學習時，決策樹的末端節點通常又叫作葉節點，其中存放着分類問題的類标簽。每個非葉節點都對應特征值之間的一個布爾條件判斷。Scikit-learn使用基尼不純度（Gini impurity）和熵作為資訊的衡量名額。實際上，這兩種名額衡量的是資料項被錯誤分類的機率。決策樹非常易于了解、使用、可觀化和驗證。為了形象展示決策樹，可以借助 Graphviz，該軟體可以從Graphviz官網下載下傳。此外，我們還需要安裝pydot2，指令如下。
1 $ pip3 install pydot2
2 $ pip3 install pydotplus #這一行必須執行，邀月注 下面把降雨資料分為訓練集和測試集兩部分，這需要用到scikit-learn提供的traintest split()函數，代碼如下。
1 x_train, x_test, y_train, y_test = train_test_split(x, y,
2 random_state=37) 建立DecisionTreeClassifier，代碼如下。
1 clf = tree.DecisionTreeClassifier(random_state=37) 我們可以使用scikit-learn的RandomSearchCV類來試驗各種參數的取值範圍，代碼如下。
1 params = {"max_depth": [2, None],
2 "min_samples_leaf": sp_randint(1, 5),
3 "criterion": ["gini", "entropy"]}
4 rscv = RandomizedSearchCV(clf, params)
5 rscv.fit(x_train,y_train) 經過一番搜尋後，得到的最佳成效和參數如下。
Best Train Score 0.703164923517
Test Score 0.705058763413
Best params {‘criterion’: ‘gini’, ‘max_depth’: 2, ‘min_samples_leaf’:2} 在任何情況下，即使隻是想驗證我們的設想，都不妨将決策樹可觀化，為此，可以使用下以下代碼來繪制決策樹的圖像。
1 sio = io.StringIO()
2 tree.export_graphviz(rscv.best_estimator_, out_file=sio,
3 feature_names=['day-of-year','yest'])
4 dec_tree = pydot.graph_from_dot_data(sio.getvalue())
5
6 print("Best Train Score", rscv.best_score_)
7 print("Test Score", rscv.score(x_test, y_test))
8 print("Best params", rscv.best_params_) 最終結果如圖10-9所示。
在上面的非葉節點中，各種判斷條件都被放置在最上一行，如果條件成立，就進入左側的子節點；否則，就進入右側的子節點。當進入一個葉節點時，節點内最下面一行中具有最大值的類别獲勝。也就是說，決策樹預測輸入樣本屬于該類别。
1 from sklearn.model_selection import train_test_split
2 from sklearn import tree
3 from sklearn.model_selection import RandomizedSearchCV
4 from scipy.stats import randint as sp_randint
5 import pydotplus as pydot
6 import io
7 import numpy as np
8 # from tempfile import NamedTemporaryFile
9
10 rain = .1 * np.load('rain.npy')
11 rain[rain < 0] = .05/2
12 dates = np.load('doy.npy').astype(int)
13 x = np.vstack((dates[:-1], np.sign(rain[:-1])))
14 x = x.T
15
16 y = np.sign(rain[1:])
17
18 x_train, x_test, y_train, y_test = train_test_split(x, y,
19 random_state=37)
20
21 clf = tree.DecisionTreeClassifier(random_state=37)
22 params = {"max_depth": [2, None],
23 "min_samples_leaf": sp_randint(1, 5),
24 "criterion": ["gini", "entropy"]}
25 rscv = RandomizedSearchCV(clf, params)
26 rscv.fit(x_train,y_train)
27
28 sio = io.StringIO()
29 tree.export_graphviz(rscv.best_estimator_, out_file=sio,
30 feature_names=['day-of-year','yest'])
31 dec_tree = pydot.graph_from_dot_data(sio.getvalue())
32
33 print("Best Train Score", rscv.best_score_)
34 print("Test Score", rscv.score(x_test, y_test))
35 print("Best params", rscv.best_params_)
36
37 # print(help(dec_tree.create_png))
38
39 #以下解決路徑問題，邀月注：InvocationException: GraphViz’s executables not found
40 import os
41 os.environ["PATH"] += os.pathsep + 'D:/tools/Graphviz/bin/' #注意修改你的路徑
42
43 # from IPython.display import Image
44 display(Image(dec_tree.create_png())) 10.11小結
本章緻力于介紹預測模組化和機器學習。對于本章涉及的這些主題，讀者如果有興趣，可以進一步參考Packt出版的相關圖書。預測性分析通常利用多種技術（其中包括機器學習），才能做出有效的預測，如判斷明天是否有雨等。
SVM能夠将資料映射到多元空間，是以，我們通過它可以把分類問題簡化為尋找最佳的一個或多個超平面來區隔資料點，進而達到分類目的。
彈性網絡正則化則是LASSO方法和嶺回歸方法的線性組合。處理回歸問題時，拟合優度通常是通過判定系數，即所謂的R平方來确定的。一些聚類算法需要推測聚類數，而另一些則沒有這樣的要求。
對于遺傳算法，我們要做的第一步是根據随機個體和遺傳資訊的表示方法來初始化種群。對于每一代個體，都會根據預定義的适應度函數選出一些個體用于配種。機器學習領域中，運用決策樹通常叫作決策樹學習。
第11章将向大家介紹互用性和雲計算。
第10章完。
随書源碼官方下載下傳：
https://www.ptpress.com.cn/shopping/buy?bookId=bae24ecb-a1a1-41c7-be7c-d913b163c111
需要登入後免費下載下傳。
邀月注：本文版權由邀月和部落格園共同所有，轉載請注明出處。
助人等于自助! [email protected]