《數據分析實戰-托馬茲.卓巴斯》讀書筆記第6章-迴歸模型

python學習筆記-目錄索引

 

 第6章涵蓋了許多回歸模型，有線性的，也有非線性的。咱們還會複習隨機森林和支持向量機，它們可用來解決分類或迴歸問題。

本章會介紹一些技術，以預測發電廠生產的電量。你將學習如下主題：
·識別並解決數據中的多重共線性
·構建線性迴歸模型，預測發電廠生產的電量
·使用OLS預測生產的電量
·使用CART估算髮電廠生產的電量
·將kNN模型用於迴歸問題
·將隨機森林模型用於迴歸分析
·使用SVM預測發電廠生產的電量
·訓練神經網絡，預測發電廠生產的電量

 

6.1導論

現實世界中一個常見的問題就是預測數量，或者用更抽象的說法，是在自變量和因變量之間找到聯繫。本章中，咱們專一於預測發電廠生產的電量。
本章使用的數據集來自於美國聯邦能源信息管理局（U.S.Energy Information Adminis-tration）。從網站上下載2014年的數據：http://www.eia.gov/electricity/data/eia923/xls/f923_2014.zip。
下載頁面：https://www.eia.gov/electricity/data/eia923/ 右邊選擇ZIP下載便可。
咱們僅使用EIA923_Schedules_2_3_4_5_M_12_2014_Final_Revision.xlsx文件中Generation and Fuel Data這一sheet的數據。咱們將預測淨產量（兆瓦時）。全部數據都是類別變量（州名或燃料類型），咱們決定將它們硬編碼。
最終，咱們的數據集僅有4494行記錄，是整個數據集的一個子集。咱們僅選擇少數州2014年超出100MWh的發電廠。咱們也只選取使用特定燃料類型的發電廠——Aggregate Fuel Code（AER）：coal（COL），distillate（DFO），hydroelectric conventional（HYC），biogenic municipal solid waste and landfill gas（MLG），natural gas（NG），nuclear（NUC），solar（SUN）以及wind（WND）。另外，只選取特定發動機的發電廠：combined-cycle combustion turbine part（CT），combustion turbine（GT），hydraulic turbine（HY），internal combustion（IC），photovoltaic（PV），steam turbine（ST）以及wind turbine（WT）。
發電廠的電量數據方差很大，而且高度偏斜，中位數是18219MWh，平均數是448213MWh；25%的發電廠（2017家）在2014年的電量少於3496MWh，而美國人均用電量是12kWh左右（2010年的數據，http://energyalmanac.ca.gov/electricity/us_per_capita_electricity-2010.html），幾乎至關於291人的小鎮，以下圖所示：


大部分發電廠使用自然氣（超過1300家），將近750家使用水力或光生髮電，接近580家使用風力發電。這樣，咱們獲得了咱們數據集中使用可再生能源的發電廠比例：

 

NG自然氣，HYC水力，SUN太陽能、光生，WND風力

接近2500家發電廠使用蒸汽渦輪生產電力，緊接着的是水輪機和燃氣渦輪：

ST蒸汽渦輪，HY水輪，GT燃氣輪機

咱們數據集中加利福尼亞的發電廠最多，超過1200家。德克薩斯以超過500家的數量位居第二，排第三的是馬薩諸塞：



6.2識別並解決數據中的多重共線性

多重共線性指的是這樣的場景，一個（或多個）自變量能夠表示爲其餘自變量的線性組合。
例如，假設咱們用州中的人口、家庭數和發電廠數來預測用電量。這種狀況下，能夠想見，州中人口越多，家庭數應該也越多，也就是說，家庭數能夠表示爲州中人口的某些（接近）線性關係。
如今，若是咱們要基於多重共線的數據估算模型，頗有可能會有一個（甚至全部的共線變量）是不顯著的。將這個變量移除（並只保留和因變量最相關的變量，這樣，保留了大部分對方差的貢獻），並不會下降模型的解釋能力。
準備:需裝好pandas、NumPy和Scikit。咱們使用Scikit下降維度。
要判斷咱們的數據是否共線，咱們須要查看自變量協方差矩陣的特徵值；若是數據共線，對協方差矩陣的分析有助於辨別變量（regression_multicollinearity.py文件）：

 1 # this is needed to load helper from the parent folder
 2 import sys
 3 sys.path.append('..')
 4 
 5 # the rest of the imports
 6 import helper as hlp
 7 import pandas as pd
 8 import numpy as np
 9 import sklearn.decomposition as dc
10 
11 def reduce_PCA(x, n):
12     '''
13         Reduce the dimensions using Principal Component
14         Analysis
15     '''
16     # create the PCA object
17     pca = dc.PCA(n_components=n, whiten=True)
18 
19     # learn the principal components from all the features
20     return pca.fit(x)
21 
22 # the file name of the dataset
23 r_filename = '../../Data/Chapter06/power_plant_dataset.csv'
24 
25 # read the data
26 csv_read = pd.read_csv(r_filename)
27 
28 x = csv_read[csv_read.columns[:-1]].copy()
29 y = csv_read[csv_read.columns[-1]]
30 
31 # produce correlation matrix for the independent variables
32 #生成自變量協方差矩陣
33 corr = x.corr()
34 
35 # and check the eigenvectors and eigenvalues of the matrix
36 #檢查矩陣的特徵向量與特徵值
37 w, v = np.linalg.eig(corr)
38 print('Eigenvalues: ', w)
39 
40 # values that are close to 0 indicate multicollinearity
41 #接近0的值意味着多重共載
42 s = np.nonzero(w < 0.01)
43 # inspect which variables are collinear
44 #找出共線變量
45 print('Indices of eigenvalues close to 0:', s[0])
46 
47 all_columns = []
48 for i in s[0]:
49     print('\nIndex: {0}. '.format(i))
50 
51     t = np.nonzero(abs(v[:,i]) > 0.33)
52     all_columns += list(t[0]) + [i]
53     print('Collinear: ', t[0])

原理：首先，咱們和往常同樣載入必需的模塊並讀入數據集。而後，咱們將數據集拆分紅自變量x（咱們爲原始的數據集建立一個副本）和因變量y。
接下來，咱們生成全部自變量之間的協方差矩陣，使用NumPy在linalg模塊中的eig（...）方法找到特徵值和特徵向量。eig（...）方法要求輸入一個方陣。
只有方陣纔有特徵向量和特徵值。非方形的矩陣有奇異值（https://www.math.washington.edu/~greenbau/Math_554/Course_Notes/ch1.5.pdf）。
接近於0的特徵值意味着（多重）共線性。要找到這些值（咱們的例子中，<0.01）的索引，咱們先建立一個真值向量：每一個元素不是False（0）就是True（1）。簡單一句w<0.001就能夠完成；這個操做的結果就是到原始向量w的等長真值向量，這個向量中小於0.001的值標爲True。而後咱們可使用NumPy的.nonzero（...）方法返回非零元素的索引列表（咱們的例子中，返回的就是True元素的索引）。你頗有可能在不止一個位置得到特徵值接近0的元素。對於咱們的數據集，你可能獲得相似的結果：

/*
Indices of eigenvalues close to 0: [28 29 30 31 32]
*/

咱們遍歷索引，找到相應的特徵向量。特徵向量中顯著大於0的元素會幫咱們找到共線的變量。因爲特徵向量能夠有負值，咱們返回絕對值>0.33的元素索引：

/*
Index: 28.
Collinear:  [0 1 3]

Index: 29.
Collinear:  [ 9 11 12 13]

Index: 30.
Collinear:  [15]

Index: 31.
Collinear:  [ 2 10]

Index: 32.
Collinear:  [ 4 14]
*/

注意這些變量的重複。咱們將這些值存到all_columns列表。等咱們遍歷全部接近0的特徵值，咱們來看看哪些變量是共線的：

for i in np.unique(all_columns):
    print('Variable {0}: {1}'.format(i, x.columns[i]))

咱們遍歷all_columns列表中去重後的值（使用NumPy的.unique（...）方法），輸出x中全部列相應的變量名：

/*
Variable 0: fuel_aer_NG
Variable 1: fuel_aer_DFO
Variable 2: fuel_aer_HYC
Variable 3: fuel_aer_SUN
Variable 4: fuel_aer_WND
Variable 9: mover_GT
Variable 10: mover_HY
Variable 11: mover_IC
Variable 12: mover_PV
Variable 13: mover_ST
Variable 14: mover_WT
Variable 15: state_CA
Variable 28: state_OH
Variable 29: state_GA
Variable 30: state_WA
Variable 31: total_fuel_cons
Variable 32: total_fuel_cons_mmbtu
*/

如今咱們能清楚看出來數據集中的共線性了：某種發動機固然與某種燃料相關。好比，風力發動機用的是風，光能發動機用的是太陽能，水力發電靠的是水。另外，加利福尼亞州也出現了；加利福尼亞是擁有最多風力渦輪機的州（第二位是德克薩斯，http://www.awea.org/resources/statefactsheets.aspx?itemnumber=890）。

（除了刪減模型中的變量而外）解決多重共線性的一個方法就是對數據集降維。咱們能夠經過PCA作到這一點（參考本書5.3節）。
咱們重用reduce_PCA（...）方法：

 1 # and reduce the data keeping only 5 principal components
 2 n_components = 5
 3 z = reduce_PCA(x, n=n_components)
 4 pc = z.transform(x)
 5 
 6 # how much variance each component explains?
 7 #每一個成分對方差貢獻多少大？
 8 print('\nVariance explained by each principal component: ',
 9     z.explained_variance_ratio_)
10 
11 # and total variance accounted for
12 #一共貢獻多少方差
13 print('Total variance explained: ',
14     np.sum(z.explained_variance_ratio_))

咱們但願獲得5個PC（principal components，主成分）。reduce_PCA（...）方法估算模型，.transform（...）方法轉換咱們的數據，返回5個主成分。咱們能夠查看每一個主成分對總體方差的貢獻：
第一個主成分大約貢獻了14.5%，第二個約13%，第三個約11%。剩下兩個的貢獻加起來大約16%，因此一共是：

/*
Variance explained by each principal component:  [0.14446578 0.13030196 0.11030824 0.0935897  0.06694611]
Total variance explained:  0.5456117974906106
*/

咱們將主成分和其餘變量一塊兒存到一個文件中：

 1 # append the reduced dimensions to the dataset
 2 for i in range(0, n_components):
 3     col_name = 'p_{0}'.format(i)
 4     x[col_name] = pd.Series(pc[:, i])
 5     
 6 x[csv_read.columns[-1]] = y
 7 csv_read = x
 8 
 9 # output to file
10 w_filename = '../../Data/Chapter06/power_plant_dataset_pc.csv'
11 with open(w_filename, 'w',newline='') as output:
12     output.write(csv_read.to_csv(index=False))

首先，咱們將主成分附加到x數據集後面；咱們但願將因變量放在最後。而後，咱們將數據集輸出到power_plant_dataset_pc.csv供後續使用。

6.3構建線性迴歸模型

線性迴歸無疑是最易於構建的模型。當你知道因變量和自變量之間是線性關係時就能夠選擇。
準備：需裝好pandas、NumPy和Scikit。
用Scikit估算迴歸模型很簡單（regression_linear.py文件）：

 1 # this is needed to load helper from the parent folder
 2 import sys
 3 sys.path.append('..')
 4 
 5 # the rest of the imports
 6 import helper as hlp
 7 import pandas as pd
 8 import numpy as np
 9 import sklearn.linear_model as lm
10 
11 @hlp.timeit
12 def regression_linear(x,y):
13     '''
14         Estimate a linear regression估算線性迴歸
15     '''
16     # create the regressor object建立迴歸對象
17     linear = lm.LinearRegression(fit_intercept=True,
18         normalize=True, copy_X=True, n_jobs=-1)
19 
20     # estimate the model估算模型
21     linear.fit(x,y)
22 
23     # return the object
24     return linear
25 
26 # the file name of the dataset
27 r_filename = '../../Data/Chapter06/power_plant_dataset_pc.csv'
28 
29 # read the data
30 csv_read = pd.read_csv(r_filename)
31 
32 # select the names of columns
33 dependent = csv_read.columns[-1]
34 independent_reduced = [
35     col
36     for col
37     in csv_read.columns
38     if col.startswith('p')
39 ]
40 
41 independent = [
42     col
43     for col
44     in csv_read.columns
45     if      col not in independent_reduced
46         and col not in dependent
47 ]
48 
49 # split into independent and dependent features
50 #拆成自變量和因變量特徵
51 x     = csv_read[independent]
52 x_red = csv_read[independent_reduced]
53 y     = csv_read[dependent]
54 
55 
56 # estimate the model using all variables (without PC)
57 #使用全部變量估算模型
58 regressor = regression_linear(x,y)
59 
60 # print model summary
61 print('\nR^2: {0}'.format(regressor.score(x,y)))
62 coeff = [(nm, coeff)
63     for nm, coeff
64     in zip(x.columns, regressor.coef_)]
65 intercept = regressor.intercept_
66 print('Coefficients: ', coeff)
67 print('Intercept', intercept)
68 print('Total number of variables: ',
69     len(coeff) + 1)

原理：按照慣例，咱們先導入必需的模塊並讀入數據集。而後將數據拆成子集。

首先，咱們取出表明因變量的列名（數據集的最後一列）、主成分（independent_reduced）以及自變量。選取主成分時，因爲全部的主成分列都以p_開頭（參見本書6.2節中的「更多」部分），咱們使用Python內建的.startswith（...）方法。選取自變量時，咱們取出全部既不是independent_reduced也不是dependent的列。而後咱們從數據集中提取這些列。
咱們調用regression_linear（...）方法估算模型。這個方法接受兩個參數：自變量x和因變量y。
regression_linear（...）方法內部，咱們先用.LinearRegression（...）方法建立模型對象。咱們指定模型要估算常數（fit_intercept=True），正規化數據，並複製自變量。
Python（默認狀況下）傳參是傳引用（不像C或C++那樣是傳值），因此方法改變了傳入的數據，這一可能性老是存在的。http://www.python-course.eu/passing_arguments.php。
咱們也將並行跑的任務數設爲計算機的核數。儘管對於線性迴歸來講，這通常不是個問題。
有了模型對象，咱們應用（.fit（...）方法）數據，即估算並返回模型。
線性迴歸模型的性能能夠R2(R平方)的形式表示。這個指標能夠理解爲模型對方差的貢獻的度量。
在本書6.4節，咱們會介紹如何手動計算R2。
LinearRegression對象提供了.score（...）方法。咱們的模型有以下的結果：

/*
R^2: 0.9965751657989851
 */

咱們使用列表表達式建立了coeff列表，以將全部變量的列表及相應的係數打印出來。咱們先（用zip（...）方法）將x的列和迴歸器的coeff_屬性壓縮在一塊兒；這建立了一個元組構成的列表，咱們將遍歷列表中的元組以建立coeff列表：

注意州的係數之間的差距有多大。另外，全部電力相關的變量（燃料和發動機）都是負的。這體現了州對咱們的模型的總體表現有一個最小的影響，特別是與截距相比時：

/*
Intercept -11531009013.842264
*/

咱們已由前一技巧知道，燃料總量是相關變量，因此咱們能夠安心地將其移除。
別忘了數據集中有多重共線性，咱們只用主成分估算模型：

 1 # estimate the model
 2 regressor_nm = regression_linear(x_no_state,y)
 3 
 4 # print model summary
 5 print('\nR^2: {0}'.format(regressor_nm.score(x_no_state,y)))
 6 coeff = [(nm, coeff)
 7     for nm, coeff
 8     in zip(x_no_state.columns, regressor_nm.coef_)]
 9 intercept = regressor_nm.intercept_
10 print('Coefficients: ', coeff)
11 print('Intercept', intercept)
12 print('Total number of variables: ',
13     len(coeff) + 1)

這個模型的R2結果更差：

/*
R^2: 0.12740548278756392
 */

因爲州的影響不大，咱們決定試試不考慮這個因素。咱們也不考慮total_fuel_consumption：

1 # removing the state variables and keeping only fuel and state
2 columns = [col for col in independent if 'state' not in col and col != 'total_fuel_cons']
3 x_no_state = x[columns]

R2有降低，但不顯著：

/*
R^2: 0.9957289744500172
*/

不過，如今咱們模型的自變量由35個（包括截距）變成了18個。有個精神要把握住，相似的表現下，你應該選解釋變量更少的模型。



要找到最佳直線（而且你只有一個解釋變量），你可使用Seaborn。儘管結果並不直觀，但咱們可使用主成分來展現。
首先，因爲咱們將主成分存在列中，咱們將它們放在棧頂，這樣數據集中只有三列——PC是主成分，x是主成分的值，y是正規化的發電量：

 1 # stack up the principal components
 2 #將主成分放在棧頂
 3 pc_stack = pd.DataFrame()
 4 
 5 # stack up the principal components
 6 #將主成分放在棧頂
 7 for col in x_red.columns:
 8     series = pd.DataFrame()
 9     series['x'] = x_red[col]
10     series['y'] = y
11     series['PC'] = col
12     pc_stack = pc_stack.append(series)

咱們先建立一個空的pandas DataFrame。而後，咱們遍歷x_red中的列（只有主成分）。DataFrame的內部序列將主成分的值放在x，將因變量的值放在y，將主成分的名字放在PC。而後咱們將序列附加到pc_stack。
如今繪製圖線：

1 # Show the results of a linear regression within each
2 # principal component
3 sns.lmplot(x='x', y='y', col='PC', hue='PC', data=pc_stack,
4            col_wrap=2, size=5)
5 
6 pl.savefig('../../Data/Chapter06/charts/regression_linear.png',
7     dpi=300)

.lmplot（...）方法將線性迴歸用到咱們的數據上。col參數指定了圖表的分組，咱們將爲5個主成分生成圖表。hue參數將爲每一個主成分更換顏色。col_wrap=2意味着咱們在一行中放兩個圖表，每一個都是5英寸高。
結果以下：

Tips:

/*
D:\tools\Python37\lib\site-packages\seaborn\regression.py:546: UserWarning: The `size` paramter has been renamed to `height`; please update your code.
  warnings.warn(msg, UserWarning)
*/

解決方案：size=5改成height=5

/*
sns.lmplot(x='x', y='y', col='fuel', hue='fuel',
    data=fuel_stack, col_wrap=2, height=5)
 */

咱們也查看兩個燃料變量，看是否存在線性關係：

# select only the fel consumption
fuel_cons = ['total_fuel_cons','total_fuel_cons_mmbtu']
x = csv_read[fuel_cons]

下面的圖顯示了，total_fuel_cons_mmbtu存在線性關係，total_fuel_cons也存在接近線性的關係：

/* 
MMBTU，million British Thermal Units表明百萬英熱單位，百萬英制熱單位。
1mmBtu=2.52X10^8cal(卡) =2.52X10^5 kcal(千卡)
1桶原油=5.8 MMBTU
 */

燃料消耗總量不考慮存儲量；只考慮使用的資源整體積（或重量）。可是，一磅的鈾或鈈比起一磅的煤，產生的能量要多得多。因此total_fuel_cons在x的低值一端有很大的偏離。當咱們考慮特定的能量（或對液體燃料來講的能量密度），total_fuel_cons_mmbtu，咱們在能量的消耗和產出之間看到了線性關係。
你能夠從圖中看出來，這條線的趨勢<1，表明着系統的損失。注意，有些燃料是可再生的，風力或太陽能發電時咱們幾乎不用承擔輸入相關（燃料）的成本。

6.4使用OLS預測生產的電量

OLS（Ordinary Least Squares，最小二乘法）也是一個線性模型。使用最小二乘法實際上也估算了一個線性迴歸。然而，即便自變量與因變量之間不是線性關係，只要參數之間是線性關係，OLS就能夠估算出一個模型。
準備:需裝好pandas和Statsmodels。
按照慣例，爲模型估算包一層函數（regression_ols.py文件）：

 1 import statsmodels.api as sm
 2 
 3 @hlp.timeit
 4 def regression_ols(x,y):
 5     '''
 6         Estimate a linear regression
 7     '''
 8     # add a constant to the data
 9     x = sm.add_constant(x)
10 
11     # create the model object
12     model = sm.OLS(y, x)
13 
14     # and return the fit model
15     return model.fit()
16 
17 # the file name of the dataset
18 r_filename = '../../Data/Chapter06/power_plant_dataset_pc.csv'
19 
20 # read the data
21 csv_read = pd.read_csv(r_filename)
22 
23 # select the names of columns
24 dependent = csv_read.columns[-1]
25 independent_reduced = [
26     col
27     for col
28     in csv_read.columns
29     if col.startswith('p')
30 ]
31 
32 independent = [
33     col
34     for col
35     in csv_read.columns
36     if      col not in independent_reduced
37         and col not in dependent
38 ]
39 
40 # split into independent and dependent features
41 x = csv_read[independent]
42 y = csv_read[dependent]
43 
44 # estimate the model using all variables (without PC)
45 regressor = regression_ols(x,y)
46 print(regressor.summary()) 

原理：首先，咱們按照慣例讀入必需的數據並選取自變量和因變量。
而後咱們估算模型。先往數據中加一個常量；這僅僅是加一列1。而後用.OLS（...）方法建立迴歸器模型；這個方法與Scikit中的對應方法不一樣，接受自變量和因變量。最後，咱們應用並返回模型。接着，打印出模型總結（這個在Scikit的線性迴歸中沒有實現）。咱們簡化了表以節省空間：

Tips:

/* 
D:\tools\Python37\lib\site-packages\numpy\core\fromnumeric.py:2495: FutureWarning: Method .ptp is deprecated and will be removed in a future version. 
Use numpy.ptp instead.
  return ptp(axis=axis, out=out, **kwargs)
*/

解決方案：無

可見，燃料類型、發動機和州名這幾個變量在統計上是不顯著的：p-value大於0.05。這些變量能夠放心地從模型中移除，不會損失太多精度。注意，數字2告訴咱們數據集中有很嚴重的多重共線性問題——這一點咱們已經知道了。
R2以及total_fuel_cons與total_fuel_cons_mmbtu統計上的顯著性證明了咱們以前的發現：僅靠這兩個變量就決定了發電廠的發電量。使用什麼燃料並無多大差異：一樣的輸入能量會獲得一樣的電量，無論是什麼燃料。其實發電廠本質上就是個能量轉換器：它將一種形式的能量轉化爲另外一種。問題來了——爲何使用化石燃料，而不是可再生能源？
咱們看下若是僅僅使用這兩個變量，模型是否是一樣有效：

1 # remove insignificant variables
2 #移除不顯著的變量
3 significant = ['total_fuel_cons', 'total_fuel_cons_mmbtu']
4 x_red = x[significant]
5 
6 # estimate the model with limited number of variables
7 # 用僅有的變量估算模型
8 regressor = regression_ols(x_red,y)
9 print(regressor.summary()) 

模型的表現並無變差：

如你所見，R2從0.997降到了0.996，差別可忽略不計。注意這個常量，儘管統計上不顯著，對最後的結果幾乎沒有影響。實際上，模型中值得保留的變量只有total_fuel_cons_mmbtu，由於它與發電量有最強的（線性）關係。
咱們也可使用MLPY估算OLS模型（regression_ols_alternative.py文件）：

 1 import mlpy as ml
 2 
 3 @hlp.timeit
 4 def regression_linear(x,y):
 5     '''
 6         Estimate a linear regression
 7     '''
 8     # create the model object
 9     ols = ml.OLS()
10 
11     # estimate the model
12     ols.learn(x, y)
13 
14     # and return the fit model
15     return ols
16 
17 # the file name of the dataset
18 r_filename = '../../Data/Chapter06/power_plant_dataset_pc.csv'
19 
20 # read the data
21 csv_read = pd.read_csv(r_filename)
22 
23 # remove insignificant variables
24 significant = ['total_fuel_cons_mmbtu']
25 x = csv_read[significant]
26 
27 # x = np.array(csv_read[independent_reduced[0]]).reshape(-1,1)
28 y = csv_read[csv_read.columns[-1]]
29 
30 # estimate the model using all variables (without PC)
31 regressor = regression_linear(x,y)
32 
33 # predict the output
34 predicted = regressor.pred(x)
35 
36 # and calculate the R^2
37 score = hlp.get_score(y, predicted)
38 print('R2: ', score)

MLPY的.OLS（）方法只能構建簡單的線性模型。咱們使用最顯著的變量，total_fuel_cons_mmbtu。因爲MLPY不提供關於模型表現的信息，咱們決定手寫一個評分函數（helper.py文件）：

 1 def get_score(y, predicted):
 2     '''
 3         Method to calculate R^2
 4     '''
 5     # calculate the mean of actuals
 6     mean_y = y.mean()
 7 
 8     # calculate the total sum of squares and residual
 9     # sum of squares
10     sum_of_square_total = np.sum((y - mean_y)**2)
11     sum_of_square_resid = np.sum((y - predicted)**2)
12 
13     return 1 - sum_of_square_resid / sum_of_square_total

這個方法遵循了決定係數的標準定義：用1減去殘差平方和與總平方和的比值。總平方和是因變量中每一個觀測值與平均值的差的平方和。殘差平方和是因變量每一個實際觀測值與模型輸出值的差的平方和。
get_score（...）方法要求兩個輸入參數：因變量的實際觀測值和模型預測的值。
對於本模型，咱們算出下述值：

/*
The method regression_linear took 0.00 sec to run.
R2:  0.9951670803634637
 */

參考:訪問這個網址，看看OLS的另外一種作法：http://www.datarobot.com/blog/ordinary-least-squares-in-python/。

6.5使用CART估算髮電廠生產的電量

以前，咱們用決策樹給銀行電話分類過（參考本書3.6節）。對於迴歸問題來講，分類樹和迴歸樹是一碼事。
準備：需裝好pandas和Scikit。
用Scikit估算CART很簡單（regression_cart.py文件）：

 1 import sklearn.tree as sk
 2 
 3 @hlp.timeit
 4 def regression_cart(x,y):
 5     '''
 6         Estimate a CART regressor
 7     '''
 8     # create the regressor object
 9     cart = sk.DecisionTreeRegressor(min_samples_split=80,
10         max_features="auto", random_state=66666,
11         max_depth=5)
12 
13     # estimate the model
14     cart.fit(x,y)
15 
16     # return the object
17     return cart

原理：按照慣例，咱們先載入數據，提取出因變量y和自變量x。一樣以主成分的形式將自變量存到x_red。
估算模型時，和以前全部用Scikit估算的模型同樣，咱們建立模型對象。.DecisionTree Regressor（...）方法接受一系列參數。min_samples_split是節點中爲了實現拆分至少要有的樣本數。max_features是樹可使用的最多變量數；auto參數是說，若是須要的話，模型可使用全部的自變量。random_state決定了模型的初始隨機狀態，max_depth決定了樹最多能有多少層。
而後應用數據，返回估算的模型。咱們看看模型表現如何：

# print out the results
print('R2: ', regressor.score(x,y))

結果看上去和以前差很少：

/*
The method regression_cart took 0.01 sec to run.
R2:  0.9595455136505057

The method regression_cart took 0.01 sec to run.
R:  0.7035525409649943
0. p_0: 0.17273577812298352
1. p_1: 0.08688236875805344
2. p_2: 0.0
3. p_3: 0.5429677683710962
4. p_4: 0.19741408474786684
*/

決策樹的美在於，模型根本不會使用不顯著的自變量。咱們看看哪些是顯著的：

1 for counter, (nm, label) \
2     in enumerate(
3         zip(x.columns, regressor.feature_importances_)
4     ):
5     print("{0}. {1}: {2}".format(counter, nm,label))

/*

0. fuel_aer_NG: 0.0
1. fuel_aer_DFO: 0.0
2. fuel_aer_HYC: 0.0
3. fuel_aer_SUN: 0.0
4. fuel_aer_WND: 0.0
5. fuel_aer_COL: 0.0
6. fuel_aer_MLG: 0.0
7. fuel_aer_NUC: 0.0
8. mover_CT: 0.0
9. mover_GT: 0.0
10. mover_HY: 0.0
11. mover_IC: 0.0
12. mover_PV: 0.0
13. mover_ST: 0.0
14. mover_WT: 0.0
15. state_CA: 0.0
16. state_TX: 0.0
17. state_NY: 0.0
18. state_FL: 0.0
19. state_MN: 0.0
20. state_MI: 0.0
21. state_NC: 0.0
22. state_PA: 0.0
23. state_MA: 0.0
24. state_WI: 0.0
25. state_NJ: 0.0
26. state_IA: 0.0
27. state_IL: 0.0
28. state_OH: 0.0
29. state_GA: 0.0
30. state_WA: 0.0
31. total_fuel_cons: 0.0
32. total_fuel_cons_mmbtu: 1.0

*/

咱們遍歷regressor（.DecisionTreeRegressor（...））的feature_importances_屬性，查看使用zip（...）的相應名字。下面是簡略的列表：

能夠看到，模型只認爲total_fuel_cons_mmbtu是顯著的，具備100%的預測能力；其餘變量都對模型的表現沒有影響。
能夠經過最後的決策樹確認這一點：

# and export to a .dot file
sk.export_graphviz(regressor,
    out_file='../../Data/Chapter06/CART/tree.dot')

要將.dot文件轉換爲易讀格式，能夠在命令行使用dot命令。切到路徑下（假設你在Codes/Chapter6文件夾）：

要爲咱們的樹生成一個PDF文件，只需執行下面這行命令：

/*
D:\PythonLearn\數據分析實戰-托馬茲·卓巴斯\tools\GraphViz\release\bin\dot -Tpdf D:\Java2018\practicalDataAnalysis\Data\Chapter06\CART\tree.dot -o D:\Java2018\practicalDataAnalysis\Data\Chapter06\CART\tree.pdf

D:\PythonLearn\數據分析實戰-托馬茲·卓巴斯\tools\GraphViz\release\bin\dot -Tpdf D:\Java2018\practicalDataAnalysis\Data\Chapter06\CART\tree_red.dot -o D:\Java2018\practicalDataAnalysis\Data\Chapter06\CART\tree_red.pdf

*/

獲得的樹應該只有一個決定變量：


X[32]，看一下列名列表，果真就是total_fuel_cons_mmbtu啊。
如今看看若是模型只有主成分，表現如何：

1 # estimate the model using Principal Components only
2 regressor_red = regression_cart(x_red,y)
3 
4 # print out the results
5 print('R: ', regressor_red.score(x_red,y)) 

R2的值也可接受：

/*
R^2:  0.7035525409649943
*/

驚奇的是，（若是你還記得每一個主成分對方差的貢獻）第三個主成分是最重要的，而第二個是徹底沒用的：

/*
0. p_0: 0.17273577812298352
1. p_1: 0.08688236875805344
2. p_2: 0.0
3. p_3: 0.5429677683710962
4. p_4: 0.19741408474786684
 */

（一般）使用主成分的問題在於沒法直接看懂結果，要得出真實意義就有點繞（https://onlinecourses.science.psu.edu/stat505/node/54）。
獲得的樹比前一棵更復雜：

參考:這裏對這裏對CART的介紹挺好：http://www.stat.wisc.edu/~loh/treeprogs/guide/wires11.pdf。

6.6將kNN模型用於迴歸問題

第3章中用到的k鄰近模型，儘管多用於解決分類問題，但其實也能夠用於解決迴歸問題。本技巧中來談談怎麼用。
準備：需裝好pandas和Scikit。

用Scikit估算模型很簡單（regression_knn.py文件）：

 1 import sklearn.neighbors as nb
 2 import sklearn.cross_validation as cv
 3 
 4 @hlp.timeit
 5 def regression_kNN(x,y):
 6     '''
 7         Build the kNN classifier
 8     '''
 9     # create the classifier object
10     knn = nb.KNeighborsRegressor(n_neighbors=80,
11         algorithm='kd_tree', n_jobs=-1)
12 
13     # fit the data
14     knn.fit(x,y)
15 
16     #return the classifier
17     return knn

原理：首先，咱們讀入數據，並拆成因變量y和自變量x_sig；咱們只選取以前發現是顯著的變量total_fuel_cons和total_fuel_cons_mmbtu。咱們也會測試只使用主成分構建的模型；這種狀況下，咱們使用x_principal。
咱們使用Scikit的.KNeighborsRegressor（...）方法構建kNN迴歸模型。咱們將n_neighbors設爲80，指定kd_tree做爲使用的算法。
查看kd-樹是如何造成的：http://www.alglib.net/other/nearestneighbors.php。
指定並行任務數爲-1，模型會使用與處理器核數相同的任務數。
估算的模型R2分值爲0.94。
可是問題來了：模型的表現對數據來講有多敏感？咱們用Scikit的cross_validation模塊來測一下這個：

1 import sklearn.cross_validation as cv
2 # test the sensitivity of R2
3 scores = cv.cross_val_score(regressor, x_sig, y, cv=100)
4 print('Expected R2: {0:.2f} (+/- {1:.2f})'\
5     .format(scores.mean(), scores.std()**2))

首先，咱們加載模塊。而後，咱們使用.cross_val_score（...）方法測試模型。第一個參數是要測試的模型，第二個參數是自變量的集合，第三個參數是因變量的集合。最後一個參數cv，指定了摺疊的數量：原始數據集隨機劃分紅的分塊數。一般狀況下，交叉驗證只用到一個分塊，剩下的用做訓練數據集。
咱們的結果：

/*
The method regression_kNN took 0.00 sec to run.
R^2:  0.9433911310360819
Expected R^2: 0.91 (+/- 0.03)
 */

Tips:

/*
ModuleNotFoundError: No module named 'sklearn.cross_validation'
*/

解決方案：改成從 sklearn.model_selection

/*
import sklearn.model_selection as cv
 */

咱們運行獲得的R^2分數是0.94。不過，交叉驗證獲得的指望值是0.91，標準差是0.03（因此咱們運行獲得的結果在這個範圍以內）。
咱們測試只使用主成分的模型：

1 # estimate the model using Principal Components only
2 #只使用主成分估算模型
3 regressor_principal = regression_kNN(x_principal,y)
4 
5 print('R^2: ', regressor_principal.score(x_principal,y))

獲得的結果你也許會意外（特別是若是你學到的說法是R^2的取值範圍在0到1之間）：

/*
The method regression_kNN took 0.01 sec to run.
R^2:  0.4602300520714968
Expected R^2: -18.55 (+/- 14427.43)
 */

因此，咱們運行獲得的R^2分數是0.46。然而，指望值是負的（標準差極高）。
咱們怎麼會獲得負值？考慮這個例子：咱們有一組數據，差很少遵循這個關係y=x^2。如今，咱們假設有些模型（只能估算線性關係）給咱們一個遵循y_pred=3*x函數的預測。
這樣的預測，R^2的值就在-0.4左右。若是你分析下R^2的計算過程，你就會明白某些狀況下R^2的值會變成負的。
參考StackOverflow上這個回答：http://stats.stackexchange.com/questions/12900/when-is-r-squared-negative。


6.7將隨機森林模型用於迴歸分析

相似於決策樹，隨機森林也能夠用於解決迴歸問題。咱們以前用它歸類電話（參考本書3.7節）。這裏，咱們使用隨機森林預測電廠的發電量。
準備：需裝好pandas、NumPy和Scikit。
隨機森林是一組模型。這個例子的代碼來自第3章（regression_randomForest.py文件）：

 1 import sys
 2 sys.path.append('..')
 3 
 4 # the rest of the imports
 5 import helper as hlp
 6 import pandas as pd
 7 import numpy as np
 8 import sklearn.ensemble as en
 9 import sklearn.model_selection as cv
10 
11 @hlp.timeit
12 def regression_rf(x,y):
13     '''
14         Estimate a random forest regressor
15         估算隨機森林迴歸器
16     '''
17     # create the regressor object
18     random_forest = en.RandomForestRegressor(
19         min_samples_split=80, random_state=666,
20         max_depth=5, n_estimators=10)
21 
22     # estimate the model
23     random_forest.fit(x,y)
24 
25     # return the object
26     return random_forest

原理：按照慣例，咱們讀入數據並拆成自變量和因變量。後面會用於估算模型。
.RandomForestRegressor（...）方法建立了模型對象。和決策樹同樣，min_samples_split決定拆分時節點最少的樣本數。random_state決定了模型的初始隨機狀態，max_depth決定了樹最多能有多少層。最後一個參數n_estimators，告訴模型森林中估算多少樹。
估算後，咱們打印出R^2評分，測試它對輸入的敏感度，以列表形式打印出變量及其重要性：

 1 # print out the results
 2 print('R^2: ', regressor_red.score(x_red,y))
 3 
 4 # test the sensitivity of R2 測試r^2敏感度
 5 scores = cv.cross_val_score(regressor_red, x_red, y, cv=100)
 6 print('Expected R^2: {0:.2f} (+/- {1:.2f})'\
 7     .format(scores.mean(), scores.std()**2))
 8 
 9 # print features importance
10 for counter, (nm, label) \
11     in enumerate(
12         zip(x_red.columns, regressor_red.feature_importances_)
13     ):
14     print("{0}. {1}: {2}".format(counter, nm,label))    

 Tips:

/*     
Traceback (most recent call last):
  File "D:\Java2018\practicalDataAnalysis\Codes\Chapter06\regression_randomForest.py", line 79, in <module>
    x_red = csv_read[features[0]]
    
    KeyError: "None of [Int64Index([32], dtype='int64')] are in the [columns]"
*/

解決方案：
改成：x_red = csv_read.iloc[:,features[0]]

代碼的結果以下（已簡化）：

/*

The method regression_rf took 0.06 sec to run.
R^2:  0.9704592485243615
Expected R^2: 0.82 (+/- 0.21)
0. fuel_aer_NG: 0.0
1. fuel_aer_DFO: 0.0
2. fuel_aer_HYC: 0.0
......
8. mover_CT: 0.0
9. mover_GT: 0.0
10. mover_HY: 0.0
11. mover_IC: 0.0
12. mover_PV: 0.0
13. mover_ST: 1.9909417055476134e-05
14. mover_WT: 0.0
15. state_CA: 0.0
16. state_TX: 0.0
......
28. state_OH: 0.0
29. state_GA: 0.0
30. state_WA: 0.0
31. total_fuel_cons: 0.0
32. total_fuel_cons_mmbtu: 0.9999800905829446
 */

又一次，咱們（果真）看到total_fuel_cons_mmbtu起決定做用。mover_ST也出現了，但它的重要性幾乎等於沒有，因此咱們能夠移除它。
只用total_fuel_cons_mmbtu構建模型的結果以下：

/*
The method regression_rf took 0.03 sec to run.
R^2:  0.9704326205779759
Expected R^2: 0.81 (+/- 0.22)
0. total_fuel_cons_mmbtu: 1.0
 */

如你所見，R2和敏感度的差異能夠忽略不計。

6.8使用SVM預測發電廠生產的電量

SVM（Support Vector Machine，支持向量機）已流行多年。模型使用多種核技巧，自變量和因變量之間最複雜的關係也能夠建模。
本技巧中，使用人工生成的數據，咱們展現SVM的威力。
準備：需裝好pandas、NumPy、Scikit和Matplotlib。
本技巧中，咱們用4種核測試用於迴歸的SVM（regression_svm.py文件）：

 1 import sys
 2 sys.path.append('..')
 3 
 4 # the rest of the imports
 5 import helper as hlp
 6 import pandas as pd
 7 import numpy as np
 8 import sklearn.svm as sv
 9 import matplotlib.pyplot as plt
10 
11 @hlp.timeit
12 def regression_svm(x, y, **kw_params):
13     '''
14         Estimate a SVM regressor
15     '''
16     # create the regressor object
17     svm = sv.SVR(**kw_params)
18 
19     # estimate the model
20     svm.fit(x,y)
21 
22     # return the object
23     return svm
24 
25 # simulated dataset
26 x = np.arange(-2, 2, 0.004)
27 errors = np.random.normal(0, 0.5, size=len(x))
28 y = 0.8 * x**4 - 2 * x**2 +  errors
29 
30 # reshape the x array so its in a column form
31 x_reg = x.reshape(-1, 1)
32 
33 models_to_test = [
34     {'kernel': 'linear'},
35     {'kernel': 'poly','gamma': 0.5, 'C': 0.5, 'degree': 4},
36     {'kernel': 'poly','gamma': 0.5, 'C': 0.5, 'degree': 6},
37     {'kernel': 'rbf','gamma': 0.5, 'C': 0.5}
38 ]

原理：首先，咱們導入全部必需模塊。而後，咱們模擬數據。咱們先用NumPy的.arange（...）方法生成x向量；數據集從-2到2，步長爲0.004。而後引入一些偏差：.random.normal（...）方法返回和x長度相等的樣本數，x是從一個均值爲0標準差爲0.5的正態分佈中取出的。咱們的因變量這樣指定：y=0.8x4-2x2+errors。
爲了估算模型，咱們要將x轉爲列的形式；咱們使用.reshape（...）方法重塑x。
models_to_test列表包括了4個關鍵詞參數，用於regression_svm（...）方法；咱們估算一個線性核SVM，兩個多項式核SVM（一個度數爲4，一個度數爲6），以及一個RBF核SVM。
regression_svm（...）方法和全部模型估算方法的邏輯一致：先建立模型對象，而後估算模型，最後返回估算的模型。
本技巧中，咱們生成一個圖表，以具象化模型的預測。
建立圖，調整框：

1 # create chart frame
2 plt.figure(figsize=(len(models_to_test) * 2 + 3, 9.5))
3 plt.subplots_adjust(left=.05, right=.95,
4     bottom=.05, top=.96, wspace=.1, hspace=.15)

建立指定大小的圖，調整子圖：爲每一個子圖調整left和right，bottom和top邊距，以及子圖間的空間，wspace和hspace。
如今遍歷要構建的模型，加到圖中：

 1 for i, model in enumerate(models_to_test):
 2     # estimate the model
 3     regressor = regression_svm(x_reg, y, **model)
 4 
 5     # score
 6     score = regressor.score(x_reg, y)
 7 
 8     # plot the chart
 9     plt.subplot(2, 2, i + 1)
10     if model['kernel'] == 'poly':
11         plt.title('Kernel: {0}, deg: {1}'\
12             .format(model['kernel'], model['degree']))
13     else:
14         plt.title('Kernel: {0}'.format(model['kernel']))
15     plt.ylim([-4, 8])
16     plt.scatter(x, y)
17     plt.plot(x, regressor.predict(x_reg), color='r')
18     plt.text(.9, .9, ('R^2: {0:.2f}'.format(score)),
19                  transform=plt.gca().transAxes, size=15,
20                  horizontalalignment='right')

先用regression_svm（...）方法估算模型。估算好後，計算R^2分值。
而後用.subplot（...）方法指定（2乘2的網格上）子圖的位置，將當前的循環次數i加1；圖表中，子圖從1開始計數，而列表索引從0開始。若是是多項式核，圖表的標題中應體現度數；不然給出核的名稱就足夠了。
而後，咱們用原始數據x和y建立散佈圖。以後，咱們用迴歸器的預測值繪圖；爲了醒目，咱們使用紅色。在各個子圖中，y軸的範圍保持在-4到8。
而後，咱們在右上角放上R^2的值；指定座標爲（.9，.9）並使用轉換能夠作到這一點。字體設爲15，右對齊。
最後保存到文件：

plt.savefig('../../data/Chapter06/charts/regression_svm.png',
    dpi= 300)

從獲得的圖中能夠看出，RBF覈對於高度非線性數據表現有多好。咱們又一次看到了線性核的負的R2值。能夠注意到，多項式核的表現也並無足夠好：

更多：使用RBF核估算SVM的主要難點在於決定爲C（偏差懲罰參數）和gamma（控制對特徵向量差別的敏感度的RBF核參數）參數選擇何值。
最後，咱們可使用格點搜索法，一個計算密集型的作法。不過，咱們將使用Optunity，一個多種環境下各類優化器的庫。
關於Optunity能夠參考：http://optunity.readthedocs.org/en/latest/。
Anaconda不直接支持Optunity，咱們須要自行安裝。

獨立安裝 optunity
>pip install optunity

/*
Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
Collecting optunity
  Downloading https://pypi.tuna.tsinghua.edu.cn/packages/32/4d/d49876a49e105b56755eb5ba06a4848ee8010f7ff9e0f11a13aefed12063/Optunity-1.1.1.tar.gz (4.6MB)
.......
Successfully built optunity
Installing collected packages: optunity
Successfully installed optunity-1.1.1
FINISHED

 */

既然有了Optunity，咱們來看看代碼（regression_svm_alternative.py文件）。這裏我跳過了各類引入；完整的引入列表能夠查看文件：

 1 import optunity as opt
 2 import optunity.metrics as mt
 3 
 4 # simulated dataset
 5 x = np.arange(-2, 2, 0.002)
 6 errors = np.random.normal(0, 0.5, size=len(x))
 7 y = 0.8 * x**4 - 2 * x**2 +  errors
 8 
 9 # reshape the x array so its in a column form
10 x_reg = x.reshape(-1, 1)
11 
12 @opt.cross_validated(x=x_reg, y=y, num_folds=10, num_iter=5)
13 def regression_svm(
14     x_train, y_train, x_test, y_test, logC, logGamma):
15     '''
16         Estimate a SVM regressor
17     '''
18     # create the regressor object
19     svm = sv.SVR(kernel='rbf',
20         C=0.1 * logC, gamma=0.1 * logGamma)
21 
22     # estimate the model
23     svm.fit(x_train,y_train)
24 
25     # decision function
26     decision_values = svm.decision_function(x_test)
27 
28     # return the object
29     return mt.roc_auc(y_test, decision_values)

首先，相似前一個例子，咱們建立模擬數據。
而後用一個Optunity中的方法裝飾regression_svm（...）方法。.cross_validated（...）方法接受一系列參數：x是自變量，y是因變量。
num_folds是用於交叉驗證的塊數。咱們將1/10的樣本用於交叉驗證，9/10用於訓練。num_iter是交叉驗證的循環次數。
注意regression_svm（...）方法的定義變了。cross_validated（...）方法先注入訓練的x和y樣本，而後是測試的x和y，以及指定的其餘參數：咱們的例子中，logC和logGamma。對於後兩個參數咱們要爲其尋找最優值。
方法從估算模型開始：先用指定的參數建立模型，應用數據。decision_function（...）方法返回一個列表，放的是x中元素到超平面的距離。咱們返回的是在ROC曲線之下的區域。（ROC參考本書3.2節）
使用Optunity提供的maximize（...）方法找到C和gamma的最佳值。第一個參數是咱們要最大化的函數，即regression_svm（...）。num_evals指定了容許的函數評估的次數。最後兩個指定了logC和logGamma搜索空間的邊界。
優化以後，咱們就能夠用找到的值估算最終的模型並保存圖片：

 1 # and the values are...
 2 print('The optimal values are: C - {0:.2f}, gamma - {1:.2f}'\
 3     .format(0.1 * hps['logC'], 0.1 * hps['logGamma']))
 4 
 5 # estimate the model with optimal values
 6 regressor = sv.SVR(kernel='rbf',
 7             C=0.1 * hps['logC'],
 8             gamma=0.1 * hps['logGamma'])\
 9         .fit(x_reg, y)
10 
11 # predict the output
12 predicted = regressor.predict(x_reg)
13 
14 # and calculate the R^2
15 score = hlp.get_score(y, predicted)
16 print('R2: ', score)
17 
18 # plot the chart
19 plt.scatter(x, y)
20 plt.plot(x, predicted, color='r')
21 plt.title('Kernel: RBF')
22 plt.ylim([-4, 8])
23 plt.text(.9, .9, ('R^2: {0:.2f}'.format(score)),
24              transform=plt.gca().transAxes, size=15,
25              horizontalalignment='right')
26 
27 plt.savefig(
28     '../../data/Chapter06/charts/regression_svm_alt.png',
29     dpi=300
30 )

大部分代碼以前已經解釋過了，這裏再也不重複。

獲得的值應該是這樣：

/* D:\tools\Python37\lib\site-packages\optunity\metrics.py:204: RuntimeWarning: invalid value encountered in true_divide return float(TP) / (TP + FN) ... D:\tools\Python37\lib\site-packages\optunity\metrics.py:204: RuntimeWarning: invalid value encountered in true_divide return float(TP) / (TP + FN)
The optimal values are: C - 0.38, gamma - 1.44
R^2:  0.8662968996734814

 */

和前一個例子相比，R^2值上升了；如今迴歸更匹配數據了：

參考：強烈推薦熟悉一下Optunity：http://optunity.readthedocs.org/en/latest/。

6.9訓練神經網絡，預測發電廠生產的電量

神經網絡是強有力的迴歸器，能夠應用到任何數據上。然而，本例中，咱們已經知道數據遵循線性關係了，因此就使用一個簡單的模型來應用數據。
準備：需裝好pandas和PyBrain。
相似第3章中的方式，咱們指定fitANN（...）方法（regression_ann.py文件）。有些引入語句已略去：

 1 import sys
 2 sys.path.append('..')
 3 
 4 # the rest of the imports
 5 import helper as hlp
 6 import pandas as pd
 7 import pybrain.structure as st
 8 import pybrain.supervised.trainers as tr
 9 import pybrain.tools.shortcuts as pb
10 
11 @hlp.timeit
12 def fitANN(data):
13     '''
14         Build a neural network regressor
15     '''
16     # determine the number of inputs and outputs
17     inputs_cnt = data['input'].shape[1]
18     target_cnt = data['target'].shape[1]
19 
20     # create the regressor object
21     ann = pb.buildNetwork(inputs_cnt,
22         inputs_cnt * 3,
23         target_cnt,
24         hiddenclass=st.TanhLayer,
25         outclass=st.LinearLayer,
26         bias=True
27     )
28 
29     # create the trainer object
30     trainer = tr.BackpropTrainer(ann, data,
31         verbose=True, batchlearning=False)
32 
33     # and train the network
34     trainer.trainUntilConvergence(maxEpochs=50, verbose=True,
35         continueEpochs=2, validationProportion=0.25)
36 
37     # and return the regressor
38     return ann

原理：咱們先讀入數據，拆分數據，並建立ANN訓練數據：

 1 # split the data into training and testing
 2 train_x, train_y, \
 3 test_x,  test_y, \
 4 labels = hlp.split_data(csv_read,
 5     y='net_generation_MWh', x=['total_fuel_cons_mmbtu'])
 6 
 7 # create the ANN training and testing datasets
 8 training = hlp.prepareANNDataset((train_x, train_y),
 9     prob='regression')
10 testing  = hlp.prepareANNDataset((test_x, test_y),
11     prob='regression')

自變量只使用total_fuel_cons_mmbtu。注意.prepareANNDataset（...）方法使用了新參數。
prob（表明problem，問題）指定了要用數據集解決一個迴歸問題：對於分類問題，咱們須要有類別指標及其到1的補充abs（item-1），參考本書3.8節。
如今能夠應用模型了：

1 # train the model
2 regressor = fitANN(training)
3 
4 # predict the output from the unseen data
5 predicted = regressor.activateOnDataset(testing)
6 
7 # and calculate the R^2
8 score = hlp.get_score(test_y, predicted[:, 0])
9 print('R^2: ', score)

fitANN（...）方法先建立神經網絡。咱們的神經網絡有一個輸入神經元，三個（帶tanh激活函數的）隱藏層神經元，一個（帶線性激活函數的）輸出神經元。
和第3章同樣使用在線訓練（batchlearning=False）。咱們指定maxEpochs爲50，即便到時候網絡沒有達到最小的偏差，咱們也終止訓練。咱們將四分之一的數據用於驗證。
訓練以後，咱們經過預測發電量來測試網絡，計算R^2：

/*
The method fitANN took 68.70 sec to run.
R^2:  0.9843443717345242 
*/

能夠看出，咱們神經網絡的表現不輸其餘模型。
參考
Tips:

/* Traceback (most recent call last):
  File "D:\Java2018\practicalDataAnalysis\Codes\Chapter06\regression_ann.py", line 8, in <module>
    import pybrain.structure as st
  File "D:\tools\Python37\lib\site-packages\pybrain\__init__.py", line 1, in <module>
    from structure.__init__ import *
ModuleNotFoundError: No module named 'structure'
 */

解決方案（均失敗）：

/*
pip install git+https://github.com/pybrain/pybrain.git
This seems to be an issue with the library when used from Python 3 and installed using pip. The latest version listed is 0.3.3, but I'm still getting v0.3.0 when installing.

Try installing numpy and scipy in Anaconda, then (after activating the conda environment) install directly from the git repo:

pip install git+https://github.com/pybrain/pybrain.git@0.3.3


嘗試方法1：pip install git+https://github.com/pybrain/pybrain.git@0.3.3

  Running command git clone -q https://github.com/pybrain/pybrain.git 'd:\Temp\pip-req-build-0soc1yjx'
Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
Collecting git+https://github.com/pybrain/pybrain.git@0.3.3
  Cloning https://github.com/pybrain/pybrain.git (to revision 0.3.3) to d:\temp\pip-req-build-0soc1yjx
  Running command git checkout -q 882cf66a89b6b6bb8112f41d200c55461464aa06
Requirement already satisfied: scipy in d:\tools\python37\lib\site-packages (from PyBrain==0.3.1) (1.3.3)
Requirement already satisfied: numpy>=1.13.3 in d:\tools\python37\lib\site-packages (from scipy->PyBrain==0.3.1) (1.17.4)
Building wheels for collected packages: PyBrain
  Building wheel for PyBrain (setup.py): started
  Building wheel for PyBrain (setup.py): finished with status 'done'
  Created wheel for PyBrain: filename=PyBrain-0.3.1-cp37-none-any.whl size=472111 sha256=07f80767b5065f40565e595826df759d2dc66740d3fb6db709d7343cbdc2a8c0
  Stored in directory: d:\Temp\pip-ephem-wheel-cache-h4bvm49q\wheels\a3\1c\d1\9cfc0e65ef0ed5559fd4f2819e4423e9fa4cfe02ff3a5c3b3e
Successfully built PyBrain
Installing collected packages: PyBrain
  Found existing installation: PyBrain 0.3
    Uninstalling PyBrain-0.3:
      Successfully uninstalled PyBrain-0.3
Successfully installed PyBrain-0.3.1
FINISHED

嘗試方法2：
D:\tools\Python37\Scripts>pip3  install https://github.com/pybrain/pybrain/archive/0.3.3.zip
Looking in indexes: https://pypi.tuna.tsinghua.edu.cn/simple
Collecting https://github.com/pybrain/pybrain/archive/0.3.3.zip
  Downloading https://github.com/pybrain/pybrain/archive/0.3.3.zip
     / 1.9MB 163kB/s
Requirement already satisfied: scipy in d:\tools\python37\lib\site-packages (from PyBrain==0.3.1) (1.3.3)
Requirement already satisfied: numpy>=1.13.3 in d:\tools\python37\lib\site-packages (from scipy->PyBrain==0.3.1) (1.17.4)
Building wheels for collected packages: PyBrain
  Building wheel for PyBrain (setup.py) ... done
  Created wheel for PyBrain: filename=PyBrain-0.3.1-cp37-none-any.whl size=468236 sha256=ce9591da503451919d71cd322324930d3d2259c7fdbe7c3041a8309bcf1ff9b2
  Stored in directory: d:\Temp\pip-ephem-wheel-cache-ouwcf150\wheels\0b\04\38\2f174aa3c578350870947ca6ab12e0eb89aef3478c9610eb0a
Successfully built PyBrain
Installing collected packages: PyBrain
Successfully installed PyBrain-0.3.1

D:\tools\Python37\Scripts>


*/

最終反覆檢查：

pip install https://github.com/pybrain/pybrain/archive/0.3.3.zip
安裝0.3.3，其實仍是0.3.1

/*
 File "D:\tools\Python37\lib\site-packages\pybrain\tools\functions.py", line 4, in <module>
    from scipy.linalg import inv, det, svd, logm, expm2
ImportError: cannot import name 'expm2' from 'scipy.linalg' (D:\tools\Python37\lib\site-packages\scipy\linalg\__init__.py)
 */

解決方案補充：

根據提示錯誤是「D:\tools\Python37\lib\site-packages\pybrain\tools\functions.py」文件中的import scipy.linalg import inv, det, svd, logm, expm2出錯，查找到一些英文資料說明SciPy中expm2和expm3魯棒性不好，官方建議再也不使用這兩個，而是使用expm，所以，在前的所提到的路徑中找到functions.py文件，找到語句「from scipy.linalg import inv, det, svd, logm, expm2」，將expm2更改爲expm，保存，退出，再次嘗試，引用成功。

/* Total error:  0.0234644034879
Total error:  0.00448033176755
Total error:  0.00205475261566
Total error:  0.000943766244883
Total error:  0.000441922091808
Total error:  0.00023137206216
Total error:  0.00014716265305
Total error:  0.000113546124273
Total error:  0.000100591874178
Total error:  9.48731641822e-05
Total error:  9.17187548879e-05
Total error:  9.00545474786e-05
Total error:  8.8754314954e-05
Total error:  8.78622211966e-05
Total error:  8.64148298491e-05
Total error:  8.52580031627e-05
Total error:  8.45496751335e-05
Total error:  8.31773350186e-05
Total error:  8.22848381841e-05
Total error:  8.13736136075e-05
Total error:  8.01448506931e-05
Total error:  7.91777696701e-05
Total error:  7.82493713543e-05
Total error:  7.75330633562e-05
Total error:  7.63823501292e-05
Total error:  7.55286037883e-05
Total error:  7.44680561538e-05
Total error:  7.36680418929e-05
Total error:  7.27012120487e-05
Total error:  7.18980979088e-05
Total error:  7.09005144214e-05
Total error:  6.99157006989e-05
Total error:  6.9222957236e-05
Total error:  6.82360145737e-05
Total error:  6.71353795523e-05
Total error:  6.63631396761e-05
Total error:  6.56685823709e-05
Total error:  6.53158010387e-05
Total error:  6.4518956904e-05
Total error:  6.35146379129e-05
Total error:  6.28304190101e-05
Total error:  6.22125569291e-05
Total error:  6.13769678787e-05
Total error:  6.09232354752e-05
Total error:  6.0155987076e-05
Total error:  5.94907225251e-05
Total error:  5.88562745717e-05
Total error:  5.80850744217e-05
Total error:  5.75905974532e-05
Total error:  5.65098746457e-05
Total error:  5.60159718614e-05
('train-errors:', '[0.023464 , 0.00448  , 0.002055 , 0.000944 , 0.000442 , 0.000231 , 0.000147 , 0.000114 , 0.000101 , 9.5e-05  , 9.2e-05  , 9e-05    , 8.9e-05  , 8.8e-05  , 8.6e-05  , 8.5e-05  , 8.5e-05  , 8.3e-05  , 8.2e-05  , 8.1e-05  , 8e-05    , 7.9e-05  , 7.8e-05  , 7.8e-05  , 7.6e-05  , 7.6e-05  , 7.4e-05  , 7.4e-05  , 7.3e-05  , 7.2e-05  , 7.1e-05  , 7e-05    , 6.9e-05  , 6.8e-05  , 6.7e-05  , 6.6e-05  , 6.6e-05  , 6.5e-05  , 6.5e-05  , 6.4e-05  , 6.3e-05  , 6.2e-05  , 6.1e-05  , 6.1e-05  , 6e-05    , 5.9e-05  , 5.9e-05  , 5.8e-05  , 5.8e-05  , 5.7e-05  , 5.6e-05  ]')
('valid-errors:', '[2.98916  , 0.006323 , 0.003087 , 0.001464 , 0.000735 , 0.00041  , 0.000263 , 0.000197 , 0.000168 , 0.000154 , 0.000147 , 0.000144 , 0.00014  , 0.000143 , 0.000136 , 0.000132 , 0.000131 , 0.000131 , 0.000136 , 0.000127 , 0.000126 , 0.000124 , 0.000124 , 0.000122 , 0.000123 , 0.000121 , 0.00012  , 0.000117 , 0.000116 , 0.000118 , 0.000114 , 0.000113 , 0.000112 , 0.000117 , 0.00011  , 0.000108 , 0.000107 , 0.000107 , 0.000105 , 0.000105 , 0.000105 , 0.000103 , 0.000101 , 0.000102 , 0.000103 , 9.8e-05  , 9.7e-05  , 9.6e-05  , 9.5e-05  , 9.6e-05  , 9.7e-05  , 9.4e-05  ]')
The method fitANN took 68.70 sec to run.
R^2:  0.9843443717345242
 */

查看使用PyBrain訓練神經網絡的例子：http://fastml.com/pybrain-a-simple-neural-networks-library-in-python/。

 

 第6章完。

杭州這個季節的窗外細雨打在窗上，感受有點陰冷。東方既白，還要早起找警察叔叔。

 python學習筆記-目錄索引

 

隨書源碼官方下載：
http://www.hzcourse.com/web/refbook/detail/7821/92