接著上一次的一元線性迴歸http://blog.csdn.net/lulei1217/article/details/49385531往下講，這篇文章要講解的多元線性迴歸。

1、什麼是多元線性迴歸模型？

當y值的影響因素不唯一時,採用多元線性迴歸模型。

  y =y=β0+β1x1+β2x2+...+βnxn

例如商品的銷售額可能不電視廣告投入,收音機廣告投入,報紙廣告投入有關係,可以有 sales =β0+β1*TV+β2* radio+β3*newspaper.

2、使用pandas來讀取資料

pandas 是一個用於資料探索、資料分析和資料處理的python庫

import pandas as pd
 

<pre name="code" class="python"># read csv file directly from a URL and save the results  
data = pd.read_csv('/home/lulei/Advertising.csv')

# display the first 5 rows
data.head()

這裡的Advertising.csv是來自http://www-bcf.usc.edu/~gareth/ISL/Advertising.csv。大家可以自己下載。

上面程式碼的執行結果：

     TV  Radio  Newspaper  Sales
0  230.1   37.8       69.2   22.1
1   44.5   39.3       45.1   10.4
2   17.2   45.9       69.3    9.3
3  151.5   41.3       58.5   18.5
4  180.8   10.8       58.4   12.9

上面顯示的結果類似一個電子表格，這個結構稱為Pandas的資料幀(data frame)，型別全稱：pandas.core.frame.DataFrame.

pandas的兩個主要資料結構：Series和DataFrame：

• Series類似於一維陣列，它有一組資料以及一組與之相關的資料標籤(即索引)組成。
• DataFrame是一個表格型的資料結構，它含有一組有序的列，每列可以是不同的值型別。DataFrame既有行索引也有列索引，它可以被看做由Series組成的字典。

# display the last 5 rows
data.tail()

只顯示結果的末尾5行

        TV  Radio  Newspaper  Sales
195   38.2    3.7       13.8    7.6
196   94.2    4.9        8.1    9.7
197  177.0    9.3        6.4   12.8
198  283.6   42.0       66.2   25.5
199  232.1    8.6        8.7   13.4

# check the shape of the DataFrame(rows, colums)
data.shape

檢視DataFrame的形狀,注意第一列的叫索引，和資料庫某個表中的第一列類似。

(200,4) 

3、分析資料

特徵：

• TV：對於一個給定市場中單一產品，用於電視上的廣告費用（以千為單位）
• Radio：在廣播媒體上投資的廣告費用
• Newspaper：用於報紙媒體的廣告費用

響應：

• Sales：對應產品的銷量

在這個案例中，我們通過不同的廣告投入，預測產品銷量。因為響應變數是一個連續的值，所以這個問題是一個迴歸問題。資料集一共有200個觀測值，每一組觀測對應一個市場的情況。

注意：這裡推薦使用的是seaborn包。網上說這個包的資料視覺化效果比較好看。其實seaborn也應該屬於matplotlib的內部包。只是需要再次的單獨安裝。

import seaborn as sns
import matplotlib.pyplot as plt 
# visualize the relationship between the features and the response using scatterplots
sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8)
plt.show()#注意必須加上這一句，否則無法顯示。

這裡選擇TV、Radio、Newspaper 作為特徵，Sales作為觀測值

返回的結果：

seaborn的pairplot函式繪製X的每一維度和對應Y的散點圖。通過設定size和aspect引數來調節顯示的大小和比例。可以從圖中看出，TV特徵和銷量是有比較強的線性關係的，而Radio和Sales線性關係弱一些，Newspaper和Sales線性關係更弱。通過加入一個引數kind='reg'，seaborn可以新增一條最佳擬合直線和95%的置信帶。

sns.pairplot(data, x_vars=['TV','Radio','Newspaper'], y_vars='Sales', size=7, aspect=0.8, kind='reg')
plt.show()

結果顯示如下：

4、線性迴歸模型

優點：快速；沒有調節引數；可輕易解釋；可理解。

缺點：相比其他複雜一些的模型，其預測準確率不是太高，因為它假設特徵和響應之間存在確定的線性關係，這種假設對於非線性的關係，線性迴歸模型顯然不能很好的對這種資料建模。

線性模型表示式： y=β0+β1x1+β2x2+...+βnxn 其中

• y是響應
• β0是截距
• β1是x1的系數，以此類推

在這個案例中： y=β0+β1∗TV+β2∗Radio+...+βn∗Newspaper

(1)、使用pandas來構建X(特徵向量)和y(標籤列)

scikit-learn要求X是一個特徵矩陣，y是一個NumPy向量。

pandas構建在NumPy之上。

因此，X可以是pandas的DataFrame，y可以是pandas的Series，scikit-learn可以理解這種結構。

#create a python list of feature names
feature_cols = ['TV', 'Radio', 'Newspaper']
# use the list to select a subset of the original DataFrame
X = data[feature_cols]
# equivalent command to do this in one line
X = data[['TV', 'Radio', 'Newspaper']]
# print the first 5 rows
print X.head()
# check the type and shape of X
print type(X)
print X.shape

輸出結果如下：

      TV  Radio  Newspaper
0  230.1   37.8       69.2
1   44.5   39.3       45.1
2   17.2   45.9       69.3
3  151.5   41.3       58.5
4  180.8   10.8       58.4
<class 'pandas.core.frame.DataFrame'>
(200, 3)

# select a Series from the DataFrame
y = data['Sales']
# equivalent command that works if there are no spaces in the column name
y = data.Sales
# print the first 5 values
print y.head()

輸出的結果如下：

0    22.1
1    10.4
2     9.3
3    18.5
4    12.9
Name: Sales

（2）、構建訓練集與測試集

<pre name="code" class="python"><span style="font-size:14px;">##構造訓練集和測試集
from sklearn.cross_validation import train_test_split  #這裡是引用了交叉驗證
X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)

#default split is 75% for training and 25% for testing

print X_train.shape
print y_train.shape
print X_test.shape
print y_test.shape

輸出結果如下：

(150, 3)
(150,)
(50, 3)
(50,)

注：上面的結果是由train_test_spilit()得到的，但是我不知道為什麼我的版本的sklearn包中居然報錯：

ImportError                               Traceback (most recent call last)
<ipython-input-182-3eee51fcba5a> in <module>()
      1 ###構造訓練集和測試集
----> 2 from sklearn.cross_validation import train_test_split
      3 #import sklearn.cross_validation
      4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)
      5 # default split is 75% for training and 25% for testing

ImportError: cannot import name train_test_split

處理方法：1、我後來重新安裝sklearn包。再一次呼叫時就沒有錯誤了。                   2、自己寫函式來認為的隨機構造訓練集和測試集。(這個程式碼我會在最後附上。) （3）sklearn的線性迴歸

from sklearn.linear_model import LinearRegression
linreg = LinearRegression()
model=linreg.fit(X_train, y_train)
print model
print linreg.intercept_
print linreg.coef_

輸出的結果如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)
2.66816623043
[ 0.04641001  0.19272538 -0.00349015]

# pair the feature names with the coefficients
zip(feature_cols, linreg.coef_)

輸出如下：

[('TV', 0.046410010869663267),
 ('Radio', 0.19272538367491721),
 ('Newspaper', -0.0034901506098328305)]

y=2.668+0.0464∗TV+0.192∗Radio-0.00349∗Newspaper
如何解釋各個特徵對應的係數的意義？

對於給定了Radio和Newspaper的廣告投入，如果在TV廣告上每多投入1個單位，對應銷量將增加0.0466個單位。就是加入其它兩個媒體投入固定，在TV廣告上每增加1000美元（因為單位是1000美元），銷量將增加46.6（因為單位是1000）。但是大家注意這裡的newspaper的係數居然是負數，所以我們可以考慮不使用newspaper這個特徵。這是後話，後面會提到的。

（4）、預測

y_pred = linreg.predict(X_test)
print y_pred

print type(y_pred)

輸出結果如下：

[ 14.58678373   7.92397999  16.9497993   19.35791038   7.36360284
   7.35359269  16.08342325   9.16533046  20.35507374  12.63160058
  22.83356472   9.66291461   4.18055603  13.70368584  11.4533557
   4.16940565  10.31271413  23.06786868  17.80464565  14.53070132
  15.19656684  14.22969609   7.54691167  13.47210324  15.00625898
  19.28532444  20.7319878   19.70408833  18.21640853   8.50112687
   9.8493781    9.51425763   9.73270043  18.13782015  15.41731544
   5.07416787  12.20575251  14.05507493  10.6699926    7.16006245
  11.80728836  24.79748121  10.40809168  24.05228404  18.44737314
  20.80572631   9.45424805  17.00481708   5.78634105   5.10594849]
<type 'numpy.ndarray'>

5、迴歸問題的評價測度

(1) 評價測度

對於分類問題，評價測度是準確率，但這種方法不適用於迴歸問題。我們使用針對連續數值的評價測度(evaluation metrics)。
這裡介紹3種常用的針對線性迴歸的測度。

1)平均絕對誤差(Mean Absolute Error, MAE)

(2)均方誤差(Mean Squared Error, MSE)

(3)均方根誤差(Root Mean Squared Error, RMSE)

這裡我使用RMES。

<pre name="code" class="python">#計算Sales預測的RMSE
print type(y_pred),type(y_test)
print len(y_pred),len(y_test)
print y_pred.shape,y_test.shape
from sklearn import metrics
import numpy as np
sum_mean=0
for i in range(len(y_pred)):
    sum_mean+=(y_pred[i]-y_test.values[i])**2
sum_erro=np.sqrt(sum_mean/50)
# calculate RMSE by hand
print "RMSE by hand:",sum_erro

最後的結果如下：

<type 'numpy.ndarray'> <class 'pandas.core.series.Series'>
50 50
(50,) (50,)
RMSE by hand: 1.42998147691

（2）做ROC曲線

import matplotlib.pyplot as plt
plt.figure()
plt.plot(range(len(y_pred)),y_pred,'b',label="predict")
plt.plot(range(len(y_pred)),y_test,'r',label="test")
plt.legend(loc="upper right") #顯示圖中的標籤
plt.xlabel("the number of sales")
plt.ylabel('value of sales')
plt.show()

顯示結果如下：（紅色的線是真實的值曲線，藍色的是預測值曲線）

直到這裡整個的一次多元線性迴歸的預測就結束了。

6、改進特徵的選擇
在之前展示的資料中，我們看到Newspaper和銷量之間的線性關係竟是負關係（不用驚訝，這是隨機特徵抽樣的結果。換一批抽樣的資料就可能為正了），現在我們移除這個特徵，看看線性迴歸預測的結果的RMSE如何？

依然使用我上面的程式碼，但只需修改下面程式碼中的一句即可：

#create a python list of feature names
feature_cols = ['TV', 'Radio', 'Newspaper']
# use the list to select a subset of the original DataFrame
X = data[feature_cols]
# equivalent command to do this in one line
#X = data[['TV', 'Radio', 'Newspaper']]#只需修改這裡即可<pre name="code" class="python" style="font-size: 15px; line-height: 35px;">X = data[['TV', 'Radio']]  #去掉newspaper其他的程式碼不變

# print the first 5 rowsprint X.head()# check the type and shape of Xprint type(X)print X.shape

最後的到的係數與測度如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)

2.81843904823
[ 0.04588771  0.18721008]
RMSE by hand: 1.28208957507

然後再次使用ROC曲線來觀測曲線的整體情況。我們在將Newspaper這個特徵移除之後，得到RMSE變小了，說明Newspaper特徵可能不適合作為預測銷量的特徵，於是，我們得到了新的模型。我們還可以通過不同的特徵組合得到新的模型，看看最終的誤差是如何的。

備註：

之前我提到了這種錯誤： 注：上面的結果是由train_test_spilit()得到的，但是我不知道為什麼我的版本的sklearn包中居然報錯：

ImportError                               Traceback (most recent call last)
<ipython-input-182-3eee51fcba5a> in <module>()
      1 ###構造訓練集和測試集
----> 2 from sklearn.cross_validation import train_test_split
      3 #import sklearn.cross_validation
      4 X_train,X_test, y_train, y_test = train_test_split(X, y, random_state=1)
      5 # default split is 75% for training and 25% for testing

ImportError: cannot import name train_test_split

處理方法：1、我後來重新安裝sklearn包。再一次呼叫時就沒有錯誤了。                   2、自己寫函式來認為的隨機構造訓練集和測試集。(這個程式碼我會在最後附上。) 這裡我給出我自己寫的函式：

import random 

<span style="font-family:microsoft yahei;">######自己寫一個隨機分配數的函式，分成兩份，並將數值一次儲存在對應的list中##########
def train_test_split(ylabel, random_state=1):
    import random 
    index=random.sample(range(len(ylabel)),50*random_state)
    list_train=[]
    list_test=[]
    i=0
    for s in range(len(ylabel)):
        if i in index:
            list_test.append(i)
        else:
            list_train.append(i)
        i+=1
    return list_train,list_test
 ###############對特徵進行分割#############################
feature_cols = ['TV', 'Radio','Newspaper']
X1 = data[feature_cols]
y1 = data.Sales
list_train,list_test=train_test_split(y1)#random_state的預設值是1

X1_train=X1.ix[list_train] #這裡使用來DataFrame的ix（）函式，可以將指定list中的索引的記錄全部放在一起
X1_test=X1.ix[list_test]
y1_train=y1.ix[list_train]
y1_test=y1.ix[list_test]
#######################開始進行模型的訓練########################################
linreg.fit(X1_train, y1_train)</span>

<span style="font-family:microsoft yahei;">######################預測#############</span>

<span style="font-family:microsoft yahei;">y1_pred = linreg.predict(X1_test)
print model
print linreg.intercept_
print linreg.coef_
#################評價測度##############</span>

<span style="font-family:microsoft yahei;">sum_mean1=0
for i in range(len(y1_pred)):
    sum_mean1+=(y1_pred[i]-y1_test.values[i])**2
sum_erro1=np.sqrt(sum_mean1/50)
# calculate RMSE by hand
print "RMSE by hand:",sum_erro1</span>

運算結果如下：

LinearRegression(copy_X=True, fit_intercept=True, normalize=False)
3.1066750253
[ 0.04588016  0.18078772 -0.00187699]
RMSE by hand: 1.39068687332