機器學習_ML_迴歸分析
簡單的線性迴歸的目標是對一個單一個特徵(解釋變量)x和目標變量y之間的關係來做塑模!
y=w0+w1x
那線性迴歸就可以理解成是找通過樣本點的最適合直線!
偏移值(offsets):迴歸線到樣本點的垂直線或稱殘差，即是預測的誤差
當多個解釋變量(x)，就叫做多元線性迴歸(multiple linear regression)了!
w0:當x0=1的時候，y軸的截距
下載房屋數據集
範本為波士頓房屋資訊，506個樣本，14個特徵，其中房價是我們的目標變量，其餘13個為解釋變量。
import pandas as pd
df = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/housing/housing.data', header=None, sep='\s+')
df.columns = ['CRIM', 'ZN', 'INDUS', 'CHAS',
'NOX', 'RM', 'AGE', 'DIS', 'RAD',
'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']
df.head()
確認資料下載成功
視覺化數據集中的重要特徵
探索式數據分析(Exploratory Data Analysis,EDA)
import matplotlib.pyplot as plt
import seaborn as sns
sns.set(style='whitegrid', context='notebook')
cols = ['LSTAT', 'INDUS', 'NOX', 'RM', 'MEDV']
sns.pairplot(df[cols], size=2.5);
plt.show()
散點圖矩陣能提供我們數據中成對特徵之間關係很有用的圖形化摘要!
這是一個類似線性方程式的判讀，當資料呈現左下右上的集合時，代表兩個類型有正相關，即X大Y也大，而左上右下是即代表是負相關!
為了量化特徵間的線性關係，需要建立一個相關矩陣!
相關矩陣與PCA中共變異數矩陣有密切關係!
相關矩陣可解釋成是共變異數矩陣經過比例重新調整後的版本!
相關矩陣與利用標準化計算出來的共變異數矩陣是相同的
相關矩陣度量成對特徵之間的線性關係，r=1為正相關，r=-1為負相關，r=0為無關。
import numpy as np
cm = np.corrcoef(df[cols].values.T)
sns.set(font_scale=1.5)
hm = sns.heatmap(cm,
cbar=True,
annot=True,
fmt='.2f',
annot_kws={'size': 15},
yticklabels=cols,
xticklabels=cols)
plt.show()
為了適合出一個線性迴歸，我們的預測目標『MEDV』要找出有高度相關性的特徵!
在LSTAT(非線性關係)與RM(線性關係)有不錯的關係!
實作普通最小平方線性迴歸模型
透過OLS(普通最小平方法)來估計最小化的平方垂直距離!
以梯度下降完成迴歸參數的迴歸工作
class LinearRegressionGD(object):
def __init__(self, eta=0.001, n_iter=20):
self.eta = eta
self.n_iter = n_iter
def fit(self, X, y):
self.w_ = np.zeros(1 + X.shape[1])
self.cost_ = []
for i in range(self.n_iter):
output = self.net_input(X)
errors = (y - output)
self.w_[1:] += self.eta * X.T.dot(errors)
self.w_[0] += self.eta * errors.sum()
cost = (errors**2).sum() / 2.0
self.cost_.append(cost)
return self
def net_input(self, X):
return np.dot(X, self.w_[1:]) + self.w_[0]
def predict(self, X):
return self.net_input(X)
以RM訓練預測MEDV模型
x = df[['RM']].values
y = df['MEDV'].values
from sklearn.preprocessing import StandardScaler
sc_x = StandardScaler()
sc_y = StandardScaler()
x_std = sc_x.fit_transform(x)
y_std = sc_y.fit_transform(y[:, np.newaxis]).flatten()
lr = LinearRegressionGD()
lr.fit(x_std, y_std)
plt.plot(range(1, lr.n_iter+1), lr.cost_)
plt.ylabel('SSE')
plt.xlabel('Epoch')
plt.show()
這邊可以看到，在第五輪的時候開始收斂!
確認適合度
def lin_regplot(X, y, model):
plt.scatter(X, y, c='lightblue')
plt.plot(X, model.predict(X), color='red', linewidth=2)
return
lin_regplot(x_std, y_std, lr)
plt.show()
在房間數增加的時候，房價通常會往上，但也不一定是準的，離散率過高的情況下會有誤判情況。
預估一下五間房
num_rooms_std = sc_x.transform([5.0])
price_std = lr.predict(num_rooms_std)
price_std
print("Price in $1000's: %.3f" % sc_y.inverse_transform(price_std))
以scikit-learn實作
LIB說明
from sklearn.linear_model import LinearRegression
slr = LinearRegression()
slr.fit(x, y)
slr.coef_[0]
slr.intercept_
lin_regplot(x, y, slr)
plt.show()
使用RANSAC找出強固的迴歸模型
一般來說，線性迴歸可能會受離群值的影響，而通常都是嚴重的影響!
這時候除了刪掉這個離群值，就是用RANSAC!(RANdom SAmple Consensus)
- 使用隨機個數樣本做群內值來適合出模型
- 把剩餘的樣本拿來測試，若數據落在定義的容許範圍內，那就加入群內值
- 使用新的群內值來適合模型
- 使用群內值來預估模型誤差
- 若誤差小於定義的容許值，或迭代次數已到定義數，即終止演算法!
from sklearn.linear_model import RANSACRegressor
ransac = RANSACRegressor(LinearRegression(),
max_trials=100,
min_samples=50,
residual_metric=lambda x:np.sum(np.abs(x),axis=1),
residual_threshold=5.0,
random_state=0)
ransac.fit(x, y)
max_trials:最大迭代次數
min_samples:隨機個數
residual_metric:計算樣本與迴歸線間的垂直距離絕對值
residual_threshold:樣本與迴歸線間的絕對值定義值(如上定義，小於才計入群內值)
預設scikit-learn使用MAD來估計選擇群內值，MAD代表目標變量y的絕對中位偏差。
(這部份沒有絕對!)
圖形化RANSAC
inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)
line_x = np.arange(3, 10, 1)
line_y_ransac = ransac.predict(line_x[:, np.newaxis])
plt.scatter(x[inlier_mask], y[inlier_mask], c='blue', marker='o', label='Inliers')
plt.scatter(x[outlier_mask], y[outlier_mask], c='lightgreen', marker='s', label='Outliers')
plt.plot(line_x, line_y_ransac, color='red')
plt.xlabel('Average number of rooms[RM]')
plt.ylabel('Price in 1000[MEDV]')
plt.legend(loc='upper left')
plt.show()
ransac.estimator_.coef_[0]
ransac.estimator_.intercept_
透過RANSAC，我們減少了離群值的影響!
評估線性迴歸模型的效能
透過上面的案例我們簡單的瞭解了在二元迴歸模型下的執行狀況!
下面我們要做多元的範例!
將數據分7:3
from sklearn.model_selection import train_test_split
X = df.iloc[:, :-1].values
y = df['MEDV'].values
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
透過線性迴歸計算
slr = LinearRegression()
slr.fit(X_train, y_train)
y_train_pred = slr.predict(X_train)
y_test_pred = slr.predict(X_test)
殘差圖
因為我們考量了所有的特徵，所以無法再以二元圖來呈現，但可以透過殘差來預測迴歸模型!
殘差:實際值和預測值之間的差異或垂直距離
真實目標變量值-預測回應值
plt.scatter(y_train_pred, y_train_pred - y_train, c='blue', marker='o', label='Training data')
plt.scatter(y_test_pred, y_test_pred - y_test, c='lightgreen', marker='s', label='Test')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='red')
plt.xlim([-10, 50])
plt.show()
在完美的情況下殘差會全部是零，但是這不可能的!
對於一個良好的迴歸模型，我們會期待它的錯誤是隨機分佈的，如此殘差就會隨機分散在中心線的週圍。
離中心線愈遠，代表它是偏差越大的離群值。
均方誤差 MSE:mean Squared Error
將SSE成本函數計算平均值，用它來最佳化線性迴歸!
衡量平均誤差
from sklearn.metrics import mean_squared_error
mean_squared_error(y_train, y_train_pred)
mean_squared_error(y_test, y_test_pred)
這部份可以發現，在訓練數據的部份MSE值為19.95，而測試數據的部份是27.19。
這也代表，這個模型對訓練數據過適了!
決定系數 coefficient of determination
可理解成是標準化的MSE
R2=1−SSESST
SSE:平方誤差點
SST:總平方和(解釋變量的變異數)
決定系數是0-1的一個數值，但有可能是負值，若該模型是完美適合這些數據，此時MSE值會是0而決定系數會是1
from sklearn.metrics import r2_score
r2_score(y_train, y_train_pred)
r2_score(y_test, y_test_pred)
過適!
正規化迴歸
脊迴歸:Ridge Regression(L2模式)
最小絕對壓縮挑選機制:Least Absolute Shrinkage and Selection Operator(L1)
c:Elastic Net(三個的中間，利用L1來產生稀疏性，也用L2來克服限制)
此三項正規化皆會加入超參數(alpha)來加強正規化的強度!
脊迴歸
from sklearn.linear_model import Ridge
redge = Ridge(alpha=1.0)
最小絕對壓縮挑選機制(LASSO)
from sklearn.linear_model import Lasso
redge = Lasso(alpha=1.0)
最小絕對壓縮挑選機制
from sklearn.linear_model import ElasticNet
redge = ElasticNet(alpha=1.0, l1_ratio=0.5)
要同時設置L1與L2!
另外，若L1為0，則計算結果等於LASSO
多項式迴歸
y=w0+w1x1+w2x2+...wdxd
d為階數(degree)
二元與多元比較範例
from sklearn.preprocessing import PolynomialFeatures
X = np.array([258.0, 270.0, 294.0,
320.0, 342.0, 368.0,
396.0, 446.0, 480.0, 586.0])[:, np.newaxis]
y = np.array([236.4, 234.4, 252.8,
298.6, 314.2, 342.2,
360.8, 368.0, 391.2,
390.8])
lr = LinearRegression()
pr = LinearRegression()
quadratic = PolynomialFeatures(degree=2)
X_quad = quadratic.fit_transform(X)
# 線性特徵
lr.fit(X, y)
X_fit = np.arange(250, 600, 10)[:, np.newaxis]
y_lin_fit = lr.predict(X_fit)
# 多元特徵
pr.fit(X_quad, y)
y_quad_fit = pr.predict(quadratic.fit_transform(X_fit))
plt.scatter(X, y, label='training points')
plt.plot(X_fit, y_lin_fit, label='linear fit', linestyle='--')
plt.plot(X_fit, y_quad_fit, label='quadratic fit')
plt.legend(loc='upper left')
plt.tight_layout()
plt.show()
從圖面可以看的出來，多項式模型比線性更能描述解釋變量與目標變量的關係!
y_lin_pred = lr.predict(X)
y_quad_pred = pr.predict(X_quad)
# MSE
mean_squared_error(y, y_lin_pred)
mean_squared_error(y, y_quad_pred)
# R2
r2_score(y, y_lin_pred)
r2_score(y, y_quad_pred)
可以明顯的看的出在線性與多元(多項次)的差異!
再以房屋數據做多項式塑模
x = df[['LSTAT']].values
y = df['MEDV'].values
regr = LinearRegression()
# 建立多項式特徵
quadratic = PolynomialFeatures(degree=2)
cubic = PolynomialFeatures(degree=3)
x_quad = quadratic.fit_transform(x)
x_cubic = cubic.fit_transform(x)
# linear fit
# 成線的時候每1單位的值怎麼走
x_fit = np.arange(x.min(), x.max(), 1)[:, np.newaxis]
regr = regr.fit(x, y)
# 用每1單位的值下去計算預測走勢
y_lin_fit = regr.predict(x_fit)
linear_r2 = r2_score(y,regr.predict(x))
# quadratic fit
regr = regr.fit(x_quad, y)
y_quad_fit = regr.predict(quadratic.fit_transform(x_fit))
quadratic_r2 = r2_score(y, regr.predict(x_quad))
# cubic fit
regr = regr.fit(x_cubic, y)
y_cubic_fit = regr.predict(cubic.fit_transform(x_fit))
cubic_r2 = r2_score(y, regr.predict(x_cubic))
# 圖形化
plt.scatter(x, y, label='training points', color='lightgray')
plt.plot(x_fit, y_lin_fit,
label='linear (d=1)',
color='blue',
lw=2,
linestyle=':')
plt.plot(x_fit, y_quad_fit,
label='quad (d=2)',
color='red',
lw=2,
linestyle='-')
plt.plot(x_fit, y_cubic_fit,
label='cubic (d=3)',
color='green',
lw=2,
linestyle='--')
plt.xlabel('low status of the population[LSTAT]')
plt.ylabel('price in 1000 [MEDV]')
plt.legend(loc='upper right')
plt.show()
從圖可以看的出，立方(3次方)比平方(2次方)或是簡單線性更能說明相關性，但是也要注意過適問題!
因此在實務上會建議用一個單獨的測試數據集來評估模型。
另外，我們對特徵LSTAT做對數，並且對特徵MEDV做平方根，就能將數據投影到線性特徵空間來透過線性迴歸處理。
x_log = np.log(x)
y_sqrt = np.sqrt(y)
x_fit = np.arange(x_log.min()-1, x_log.max()+1, 1)[:, np.newaxis]
regr = regr.fit(x_log, y_sqrt)
y_lin_fit = regr.predict(x_fit)
linear_r2 = r2_score(y_sqrt, regr.predict(x_log))
plt.scatter(x_log, y_sqrt, label='training points', color='lightgray')
plt.plot(x_fit, y_lin_fit,
label='linear (d=1)',
color='blue',
lw=2)
plt.xlabel('low status of the population[LSTAT]')
plt.ylabel('price in 1000 [MEDV]')
plt.legend(loc='lower left')
plt.show()
可以發現，r2的分數比上面都還要高!
轉換解釋變量到對數空間，並將目標變量取平方根，我們似乎就能夠用一條線性迴歸來描述這兩個變量之間的關係了。
用決策樹與隨機森林處理非線性關係
與上面的迴歸模型不同，隨機森林是多個決策樹的合體，是分段的線性函數的加總!
上面的模式是全域性的函數，透過決策樹，我們將輸入空間區分成更容易處理的小區域。
決策樹
from sklearn.tree import DecisionTreeRegressor
x = df[['LSTAT']].values
y = df['MEDV'].values
tree = DecisionTreeRegressor(max_depth=3)
tree.fit(x, y)
sort_idx = x.flatten().argsort()
lin_regplot(x[sort_idx], y[sort_idx], tree)
plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('price in 1000\'s [MEDV]')
plt.show()
從圖可見，決策樹可以描述數據的整體趨勢，但是要注意，它不能清楚預測連續性與可微性，並且要注意深度造成的過適問題
隨機森林
x = df.iloc[:, :-1].values
y = df['MEDV'].values
x_train, x_test, y_train, y_test = train_test_split(x, y ,test_size=0.4, random_state=1)
from sklearn.ensemble import RandomForestRegressor
forest = RandomForestRegressor(n_estimators=1000,
criterion='mse', # 求平均
random_state=1,
n_jobs=-1)
forest.fit(x_train, y_train)
y_train_pred = forest.predict(x_train)
y_test_pred = forest.predict(x_test)
mean_squared_error(y_train, y_train_pred)
mean_squared_error(y_test, y_test_pred)
r2_score(y_train, y_train_pred)
r2_score(y_test, y_test_pred)
以結果來說，有點過適了，訓練資料集的平均錯誤很漂亮，但是測試集很高
但是對目標變量與解釋變量之間的關係還是可以清楚的說明!(r2_score)
確認殘差
plt.scatter(y_train_pred,
y_train_pred - y_train,
c='black',
marker='o',
s=35,
alpha=0.5,
label='Training data')
plt.scatter(y_test_pred,
y_test_pred - y_test,
c='lightgreen',
marker='s',
s=35,
alpha=0.7,
label='Test data')
plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend('upper left')
plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='red')
plt.xlim([-10, 50])
plt.show()
殘差分佈似乎沒有圍著中心點，這表示該模型不能夠說明所有解釋變量的資訊!
