機器學習作業---偏差和方差（線性回歸）

機器學習作業---偏差和方差（線性回歸）錯誤反例，但是理清了代碼思路，很重要

一：加載數據，顯示數據

（一）數據可視化

import numpy as np
import matplotlib.pyplot as plt
import scipy.io as sio
import scipy.optimize as opt


data = sio.loadmat("ex5data1.mat")

X = data['X']
y = data['y'].flatten()
Xval = data['Xval']
yval = data['yval'].flatten()
Xtest = data['Xtest']
ytest = data['ytest'].flatten()

m = y.size

plt.figure()
plt.scatter(X,y,c='b',marker='o')
plt.xlabel("Change in water level (x)")
plt.ylabel("Water folowing out of the dam (y)")
plt.show()

（二）數據顯示

直接從文件中讀取的數據X:

[[-15.93675813]
 [-29.15297922]
 [ 36.18954863]
 [ 37.49218733]
 [-48.05882945]
 [ -8.94145794]
 [ 15.30779289]
 [-34.70626581]
 [  1.38915437]
 [-44.38375985]
 [  7.01350208]
 [ 22.76274892]]

直接從文件中讀取的數據y:

[[ 2.13431051]
 [ 1.17325668]
 [34.35910918]
 [36.83795516]
 [ 2.80896507]
 [ 2.12107248]
 [14.71026831]
 [ 2.61418439]
 [ 3.74017167]
 [ 3.73169131]
 [ 7.62765885]
 [22.7524283 ]]

可以知道都是二維數組類型。

將y展開為一維數組后：

[ 2.13431051  1.17325668 34.35910918 36.83795516  2.80896507  2.12107248
 14.71026831  2.61418439  3.74017167  3.73169131  7.62765885 22.7524283 ]

二：代價函數實現

（一）代碼實現

def reg_cost(theta,X,y,lamda):
    # 我們在這里將輸入的X前面加上一列1，用來匹配截距θ_0
    m = y.size
    X = np.c_[np.ones(m),X]
    reg_theta = theta[1:]
    #返回正則化代價函數
    hythesis = np.dot(X,theta)-y
    return np.sum(np.power(hythesis,2))/(2*m) + np.sum(np.power(reg_theta,2))*lamda/(2*m)

（二）測試結果

theta = np.ones(X.shape[1]+1)
cost = reg_cost(theta,X,y,1)
print(cost)

三：實現梯度求解參數 

（一）代碼實現

def reg_gradient(theta,X,y,lamda=1):
    m = y.size
    X = np.c_[np.ones(m),X]
    reg_theta = theta[1:]
    hythesis = np.dot(X,theta)-y
 #可以看到和代價函數還是有好多重復，所以后面推薦一個函數實現
    theta_grad = np.dot(hythesis,X)/m
    theta_grad[1:]+=theta[1:]*lamda/m

    return theta_grad

（二）結果測試

返回的是一維數組來表示參數向量 

四：使用高級優化算法求解參數向量、擬合數據

（一）使用高級優化算法求解參數

res = opt.minimize(fun=reg_cost,x0=theta,args=(X,y,1),method="TNC",jac=reg_gradient)
print(res)

（二）可視化擬合程度

res = opt.minimize(fun=reg_cost,x0=theta,args=(X,y,1),method="TNC",jac=reg_gradient)
theta_new = res.x

plt.figure()
plt.scatter(X,y,c='b',marker='o')
plt.plot(X,theta_new[1]*X+theta_new[0])
plt.xlabel("Change in water level (x)")
plt.ylabel("Water folowing out of the dam (y)")
plt.show()

五： 繪制學習曲線（重點）

重點：

1.使用訓練集的子集來擬合數據--獲取本次循環最優參數

2.在計算訓練代價（誤差）和交叉驗證代價（誤差）時，不需要使用正則化

3.使用相同的驗證集自己來計算訓練代價

（一）實現獲取訓練誤差、驗證誤差

def learning_curve(X,y,Xval,yval,lamda):
    m = y.size

    error_train = np.zeros(m)
    error_val = np.zeros(m)

    for i in range(m):
        # 獲取新的臨時訓練樣本集---重點事項1
        X_temp = X[:i+1]
        y_temp = y[:i+1]

        theta = np.ones(X.shape[1]+1)
        #獲取當前訓練集樣本下的最優參數θ向量
        res = opt.minimize(fun=reg_cost,x0=theta,args=(X_temp,y_temp,lamda),method="TNC",jac=reg_gradient)
        theta_new = res.x

        # 根據上面獲得的最優參數向量，我們獲取訓練集誤差---注意：訓練集誤差的數據集來自於我們新的訓練集
        error_train[i] = reg_cost(theta_new,X_temp,y_temp,lamda)
        # 驗證集誤差，數據集來自於全部驗證集
        # 每次都是將所有驗證集進行計算誤差----重點事項3
        error_val[i] = reg_cost(theta_new,Xval,yval,lamda)

    return error_train,error_val

#繪制學習曲線
##設置lamda為0，不進行正則化處理---重點事項2
error_train,error_val = learning_curve(X,y,Xval,yval,lamda=0)
print("error_train:\n",error_train)
print("error_val:\n",error_val)

（二）繪制學習曲線圖像  

plt.figure() plt.plot(np.arange(X.shape[0]),error_train,np.arange(X.shape[0]),error_val)
plt.title("Learning Curve for linear regression")
plt.xlabel("Number of Training Examples")
plt.ylabel("Error")
plt.legend(['Train','Cross Validation'])
plt.show()

（三）學習曲線得出結論---高偏差（欠擬合）

繪制非正則化學習曲線圖，隨着訓練樣本數量的增加，訓練誤差和交叉驗證逐漸接近，但誤差較大，高偏差、欠擬合---增大訓練樣本無用（高偏差）

六：改進機器學習系統性能的方法

當我們運用訓練好了的模型來預測未知數據的時候發現有較大的誤差，我們下一步可以做什么？

獲得更多的訓練實例——通常是有效的，但代價較大，下面的方法也可能有效，可考慮先采用下面的幾種方法：

1. 獲得更多的訓練實例——解決高方差 2. 嘗試減少特征的數量——解決高方差 3. 嘗試獲得更多的特征——解決高偏差 4. 嘗試增加多項式特征——解決高偏差 5. 嘗試減少正則化程度λ——解決高偏差 6. 嘗試增加正則化程度λ——解決高方差

我們不應該隨機選擇上面的某種方法來改進我們的算法，而是運用一些機器學習診斷法來幫助我們知道上面哪些方法對我們的算法是有效的。

（一）獲取更多訓練數據---解決高方差（后面再補充，可能有問題）

見：https://www.cnblogs.com/ssyfj/p/12872699.html（六：學習曲線）

前提是：

1.我們已經使用全部訓練集得出了一個模型，該模型可能存在過擬合或者欠擬合問題。下面將說明訓練集數據增大后對這個模型的影響。

例如：對於下面數據，我們提出下面線性模型：

 則無論我們如何改變參數θ，依舊存在欠擬合情況。

2.同時為了我們方便判斷這個模型的到底是過擬合還是欠擬合，我們通過繪制學習曲線來進行判斷。

1.增大訓練集數據，對高偏差無用---（已有模型，該模型欠擬合）

高偏差（欠擬合）：預測值與實際值之間的偏離程度

只是一條直線來減小擬合是不可能對這組數據進行很好得擬合。所以當我們給出交叉驗證集誤差隨着m得增大得圖像時，會發現，誤差會隨着m增大逐漸減小，但是會存在一個下界。

同樣訓練誤差一開始也是很小的，而在高偏差的情況下，你會發現訓練集誤差會逐漸增大，最后接近交叉驗證誤差：

2.增大數據集數據，對高方差有用 ---（已有模型，該模型過擬合）

高方差（過擬合）：預測值變化波動情況

訓練誤差---當訓練樣本越多的時候，就越難把訓練集數據擬合得更好。但總的來說訓練集誤差還是很小得

其學習曲線如下：

交叉驗證誤差---過擬合下（泛化能力很差），驗證誤差將會一直很大，盡管隨着樣本增大，但是總的來說還是很大。

其特點是：在訓練誤差和交叉驗證誤差之間有一段很大的差距。

當我們增大樣本數量時，訓練誤差會增大，而驗證誤差會減少。所以高方差下，增大樣本數量還是有用的。 

（二）多項式特征

對於訓練集，當d較小時，模型擬合程度更低，誤差較大；隨着d的增長，擬合程度提高，誤差減小。

對於交叉驗證集，當d較小時，模型擬合程度低，誤差較大；但是隨着d的增長，誤差呈現先減小后增大的趨勢，轉折點是我們的模型開始過擬合訓練數據集的時候。

 訓練集誤差和交叉驗證集誤差近似時：偏差/欠擬合---多項式次數d太小，我們應該適當增大

 交叉驗證集誤差遠大於訓練集誤差時：方差/過擬合---多項式次數d太大，我們應該適當減小

（三）正則化程度入---線性回歸舉例

1.當我們設置入較大時，比如入=1000，這時θ_1,...,θ_m都將受到很大的懲罰。所以θ_1,...,θ_m幾乎都等於0.這時H_θ(x)≈θ_0

因此這個假設處於高偏差，對數據集嚴重欠擬合。

2.與之對應的另一種情況是，如果我們的lambda值很小，比如說lambda的值等於0的時候，在這種情況下，如果我們要擬合一個高階多項式的話，那么此時我們通常會處於過擬合的情況。

3.只有當我們取一個中間大小的，既不大也不小的lambda值時，我們才會得到一組合理的，對數據剛好擬合的theta參數值

學習曲線如下：

當入較小時，訓練集誤差較小（過擬合）而交叉驗證集誤差較大--高方差

隨着入的增加，訓練集誤差不斷增加（欠擬合），而交叉驗證集誤差則是先減小后增加---高偏差

七：增加多項式特征解決欠擬合問題 

從我們在五中學習曲線中，可以看出，使用線性模型模型存在欠擬合（高偏差）狀態，所以這里我們為了解決欠擬合問題，應該嘗試增加多項式特征。

（一） 增加特征，擴展特征（從1到10階）

def poly_feature(X,p):  #傳入訓練集數據和多項式次數---因為原來只有一個x_1特征項，所以這里我們增加特征，是x_1^n,沒有中間項
    X_poly = np.zeros((X.shape[0],p))   #由原來m行，1列，變為m行p列
    p = np.arange(1,p+1)
    X_poly = X**p   #獲取多項式特征
    return X_poly

p = 8　　#這里用8階來測試我們的函數

X_poly = poly_feature(X,p)
print(X_poly.shape)
print(X_poly)

補充：兩個向量之間進行平方運算

a**b,b會對a進行列擴展。

（二）獲取擴展多項式后訓練集的均值和方差

在進行歸一化處理之前，我們需要獲取均值和方差（和我們五中，獲取theta一樣），是從訓練集中獲取數據信息，然后運用到后面訓練集、驗證集、測試集中進行歸一化處理。並且在繪制擬合圖形和學習曲線時，都會受到我們從訓練集中獲取的均值和方差影響。

def get_mean_std(X):
    fea_mean = np.mean(X,0) #對特征進行均值求解（按列）
    sigma = np.std(X,0,ddof=1)  #ddof=0求解總體標准差，ddof=1求解樣本標准差
    return fea_mean,sigma

p = 6　　#這里用6階來進行測試
lmd = 0

X_poly = poly_feature(X,p)

mean,sigma = get_mean_std(X_poly)
print("mu:\n",mean)
print("sigma:\n",sigma)

X_norm = feature_normalize(X_poly,mean,sigma)

（三）歸一化數據 （6階來測試）

對於我們進行擴展后的多項式，各個特征中的數據大小差異可能過大，我們要進行歸一化處理，防止某一個特征影響太大，從而影響其他特征的作用。這里我們使用標准差標准化。

標准差標准化： x =(x - u)/σ u是均值 σ是標准差

def feature_normalize(X,mean,sigma):
    X_norm = (X-mean)/sigma #對每一行都減去這個均值
    return X_norm

（四）對訓練數據、驗證樣本、測試樣本進行多項式處理、歸一化處理（6階測試）

X_poly = poly_feature(X,p)

mean,sigma = get_mean_std(X_poly)
print("mu:\n",mean)
print("sigma:\n",sigma)

X_norm = feature_normalize(X_poly,mean,sigma)

X_val_poly = poly_feature(Xval,p)
X_val_norm = feature_normalize(X_val_poly,mean,sigma)

X_test_poly = poly_feature(Xtest,p)
X_test_norm = feature_normalize(X_test_poly,mean,sigma)

print("X_norm:\n",X_norm[:1])
print("X_val_norm:\n",X_val_norm[:1])

（五）高級優化算法獲取多項式參數向量

theta = np.ones(p+1)　　#p=6,lmd=0

#繪制原始數據集和我們新的多項式模型---查看擬合狀態
res = opt.minimize(fun=reg_cost,x0=theta,args=(X_norm,y,lmd),method="TNC",jac=reg_gradient)
theta_new = res.x
print(theta_new)

（六）獲取訓練誤差和驗證誤差

def poly_learning_curve(X,y,Xval,yval,p_f,lamda):
    error_train = np.zeros(p_f)
    error_val = np.zeros(p_f)

    lmd = 0

    for p in range(1,p_f+1):
        # 獲取新的臨時訓練樣本集---重點事項1
        X_poly = poly_feature(X, p)
        mean, sigma = get_mean_std(X_poly)
        X_norm = feature_normalize(X_poly, mean, sigma)

        X_val_poly = poly_feature(Xval, p)
        X_val_norm = feature_normalize(X_val_poly, mean, sigma)

        theta = np.ones(p + 1)

        # 繪制原始數據集和我們新的多項式模型---查看擬合狀態
        res = opt.minimize(fun=reg_cost, x0=theta, args=(X_norm, y, lmd), method="TNC", jac=reg_gradient)
        theta_new = res.x

        # 根據上面獲得的最優參數向量，我們獲取訓練集誤差---注意：訓練集誤差的數據集來自於我們新的訓練集
        error_train[p-1] = reg_cost(theta_new,X_norm,y,lamda)
        # 驗證集誤差，數據集來自於全部驗證集
        # 每次都是將所有驗證集進行計算誤差----重點事項3
        error_val[p-1] = reg_cost(theta_new,X_val_norm,yval,lamda)

    print("error_train:\n",error_train)
    print("error_val:\n",error_val)
    return error_train,error_val

data = sio.loadmat("ex5data1.mat")

X = data['X']
y = data['y'].flatten()
Xval = data['Xval']
yval = data['yval'].flatten()
Xtest = data['Xtest']
ytest = data['ytest'].flatten()

m = y.size
theta = np.ones(X.shape[1]+1)

lmd = 0
p = 10

#繪制學習曲線
##設置lamda為0，不進行正則化處理---重點事項2
error_train,error_val = poly_learning_curve(X,y,Xval,yval,p,lmd)

（七）繪制學習曲線---用於判斷求取合適的多項式項數

#繪制學習曲線
##設置lamda為0，不進行正則化處理---重點事項2
error_train,error_val = poly_learning_curve(X,y,Xval,yval,p,lmd)

plt.figure()
plt.plot(np.arange(1,p+1),error_train,np.arange(1,p+1),error_val)
plt.title("Learning Curve for linear regression")
plt.xlabel("Number of Poly Feature")
plt.ylabel("Error")
plt.legend(['Train','Cross Validation'])
plt.show()

由六（二）可以知道，我們選取項數為3時，可以很好的擬合數據集！！！

（八）繪制擬合圖形

#產生一系列新的x坐標值，一會用於繪制擬合曲線
p = 3
theta = np.ones(p + 1)

X_poly = poly_feature(X, p)
mean, sigma = get_mean_std(X_poly)
X_norm = feature_normalize(X_poly, mean, sigma)

# 繪制原始數據集和我們新的多項式模型---查看擬合狀態
res = opt.minimize(fun=reg_cost, x0=theta, args=(X_norm, y, lmd), method="TNC", jac=reg_gradient)
theta_new = res.x

# print(np.min(X),np.max(X))  #-48.058829452570066 37.49218733199513 所以我們取-60到60
x_plot = np.linspace(np.min(X),np.max(X),100).reshape(100,1) #轉為列向量
x_plot_poly = poly_feature(x_plot,p)
x_plot_norm = feature_normalize(x_plot_poly,mean,sigma)

#獲取y值
#print(x_plot_norm.shape,theta_new.shape)    #(100, 8) (9,)
#先向我們坐標序列加一列
x_plot_norm = np.c_[np.ones(100),x_plot_norm]   #(100, 9)
y_plot = np.dot(x_plot_norm,theta_new)  #(100,)

plt.figure()
plt.scatter(X,y,c='b',marker='o')   #繪制原始數據集
#開始繪制擬合曲線
plt.plot(x_plot,y_plot,'r--')
plt.xlabel("Change in water level (x)")
plt.ylabel("Water folowing out of the dam (y)")
plt.show()

八：討論正則化程度入對數據的影響

由於我們上面使用了多項式解決了欠擬合問題，這里我們使用正則化參數入來解決一個過擬合問題！！！

從七中，我們可以知道，當我們項數小於3時，趨於欠擬合（高偏差），當我們選取的項數大於3時，趨於過擬合（高方差）。

所以我們這里選取多項式特征為6，以此為基礎來解決過擬合問題！！！

（一）設置多項式特征為6，入=0，進行過擬合測試，繪制擬合圖形

#產生一系列新的x坐標值，一會用於繪制擬合曲線
p = 6 lmd = 0

theta = np.ones(p + 1)

X_poly = poly_feature(X, p)
mean, sigma = get_mean_std(X_poly)
X_norm = feature_normalize(X_poly, mean, sigma)

# 繪制原始數據集和我們新的多項式模型---查看擬合狀態
res = opt.minimize(fun=reg_cost, x0=theta, args=(X_norm, y, lmd), method="TNC", jac=reg_gradient)
theta_new = res.x

# print(np.min(X),np.max(X))  #-48.058829452570066 37.49218733199513 所以我們取-60到60
x_plot = np.linspace(np.min(X),np.max(X),100).reshape(100,1) #轉為列向量
x_plot_poly = poly_feature(x_plot,p)
x_plot_norm = feature_normalize(x_plot_poly,mean,sigma)

#獲取y值
#print(x_plot_norm.shape,theta_new.shape)    #(100, 8) (9,)
#先向我們坐標序列加一列
x_plot_norm = np.c_[np.ones(100),x_plot_norm]   #(100, 9)
y_plot = np.dot(x_plot_norm,theta_new)  #(100,)

plt.figure()
plt.scatter(X,y,c='b',marker='o')   #繪制原始數據集
#開始繪制擬合曲線
plt.plot(x_plot,y_plot,'r--')
plt.xlabel("Change in water level (x)")
plt.ylabel("Water folowing out of the dam (y)")
plt.show()

與我們上面獲得的，多項式特征為3比較，發現在兩端有明顯過擬合現象。所以下面我們通過修改正則化程度入來進行改善！！！

（二） 獲取訓練誤差和驗證誤差

def reg_learning_curve(X,y,Xval,yval,lmds):
    error_train = np.zeros(len(lmds))
    error_val = np.zeros(len(lmds))

    for i in range(len(lmds)):
        theta = np.ones(X.shape[1]+1)
        lmd = lmds[i]
        # 繪制原始數據集和我們新的多項式模型---查看擬合狀態
        res = opt.minimize(fun=reg_cost, x0=theta, args=(X_norm, y, lmd), method="TNC", jac=reg_gradient)
        theta_new = res.x

        # 根據上面獲得的最優參數向量，我們獲取訓練集誤差---注意：訓練集誤差的數據集來自於我們新的訓練集
        error_train[i] = reg_cost(theta_new,X,y,0)
        # 驗證集誤差，數據集來自於全部驗證集
        # 每次都是將所有驗證集進行計算誤差----重點事項3
        error_val[i] = reg_cost(theta_new,Xval,yval,0)
    return error_train,error_val

重點注意：五中提及的，對於我們后面的訓練集、驗證集求解誤差時，不需要進行正則化，所以我們需要將其入設置為0！！！

#繪制學習曲線
p = 6
lmds = [0, 0.001, 0.003, 0.01, 0.03, 0.1, 0.3, 1, 3, 10]

X_poly = poly_feature(X, p)
mean, sigma = get_mean_std(X_poly)
X_norm = feature_normalize(X_poly, mean, sigma)

X_val_poly = poly_feature(Xval, p)
X_val_norm = feature_normalize(X_val_poly, mean, sigma)

# lam_error(X_norm,y,X_val_norm,yval)

error_train,error_val = reg_learning_curve(X_norm,y,X_val_norm,yval,lmds)

（三）繪制學習曲線---用於判斷求取合適的正則化程度入值

plt.figure()
plt.plot(lmds,error_train,lmds,error_val)
plt.title("Learning Curve for linear regression")
plt.xlabel("Lamda")
plt.ylabel("Error")
plt.legend(['Train','Cross Validation'])
plt.show()

由六（三）可以知道，當我們選擇正則化程度入值為3時，可以很好的防止過擬合！！！

（四）繪制擬合圖形

#產生一系列新的x坐標值，一會用於繪制擬合曲線
p = 6
lmd = 3

theta = np.ones(p + 1)

X_poly = poly_feature(X, p)
mean, sigma = get_mean_std(X_poly)
X_norm = feature_normalize(X_poly, mean, sigma)

# 繪制原始數據集和我們新的多項式模型---查看擬合狀態
res = opt.minimize(fun=reg_cost, x0=theta, args=(X_norm, y, lmd), method="TNC", jac=reg_gradient)
theta_new = res.x

# print(np.min(X),np.max(X))  #-48.058829452570066 37.49218733199513 所以我們取-60到60
x_plot = np.linspace(np.min(X),np.max(X),100).reshape(100,1) #轉為列向量
x_plot_poly = poly_feature(x_plot,p)
x_plot_norm = feature_normalize(x_plot_poly,mean,sigma)

#獲取y值
#print(x_plot_norm.shape,theta_new.shape)    #(100, 8) (9,)
#先向我們坐標序列加一列
x_plot_norm = np.c_[np.ones(100),x_plot_norm]   #(100, 9)
y_plot = np.dot(x_plot_norm,theta_new)  #(100,)

plt.figure()
plt.scatter(X,y,c='b',marker='o')   #繪制原始數據集
#開始繪制擬合曲線
plt.plot(x_plot,y_plot,'r--')
plt.xlabel("Change in water level (x)")
plt.ylabel("Water folowing out of the dam (y)")
plt.show()