12 回歸算法 - 手寫梯度下降代碼

首先回顧梯度下降：

當確定了初始值和步長后，根據(jù)梯度下降的算法，函數(shù)會不斷得迭代，最后收斂得到極值。那么收斂的條件是什么？

1、最終目標f_change = f(x)基本不發(fā)生變化的時候，我們認為取到了極值。

參考：10 回歸算法 - 梯度下降

## 原函數(shù)
def f(x):
    return x ** 2

## 首先要對f(x)進行求導 y'=2x
def h(x):
    return 2 * x

X=[]
Y=[]
x=2 #初始值
step = 0.8 #步長

f_change = f(x)
f_current = f(x)
X.append(x)
Y.append(f_current)
while f_change>1e-10:
    x = x-step * h(x)
    tmp = f(x)
    f_change = np.abs(f_current - tmp)
    f_current = tmp
    X.append(x)
    Y.append(f_current)
    print(u'x=',x)
    print(u'f_change:',f_change,'f_current=',f_current)
print(u'最終結果為',(x,f_current))
2、我可以給迭代的次數(shù)設置一個上限，比如迭代200次，可能這個時候函數(shù)還沒有收斂，但是我們認為已經(jīng)迭代了太多次，所以讓迭代中止。
案例：基于梯度下降法實現(xiàn)線性回歸算法

之前大家已經(jīng)看過了很多機器求解系數(shù)θ的模型，今天我們自己寫一個梯度下降的模型，親自體驗一下求解θ值的過程。

1、當我們獲取了一組樣本后，我們要將特征-X和目標-Y抽取出來并形成矩陣。在根據(jù) [Y|X] 求解之前，我們要先驗證數(shù)據(jù)的合法性。

# 數(shù)據(jù)校驗
def validate(X, Y):
    if len(X) != len(Y):
        raise Exception("參數(shù)異常")
    else:
        m = len(X[0])
        for l in X:
            if len(l) != m:
                raise Exception("參數(shù)異常")
        if len(Y[0]) != 1:
            raise Exception("參數(shù)異常")

2、計算差異值，對應的公式：

單個θ的梯度下降
# 計算差異值
def calcDiffe(x, y, a):
    lx = len(x) #特征長度
    la = len(a) #系數(shù)θ的長度
    
    # 系數(shù)向量θ中不包含常數(shù)項
    if lx == la: 
        result = 0
        for i in range(lx):
            result += x[i] * a[i]
        return y - result
    
    # 系數(shù)向量θ中包含常數(shù)項
    elif lx + 1 == la:
        result = 0
        for i in range(lx):
            result += x[i] * a[i]
        result += 1 * a[lx] # 加上常數(shù)項
        
        # 實際值-預測值
        return y - result
    else :
        raise Exception("參數(shù)異常")

3、梯度下降計算θ值的步驟：

  1. 隨機初始化0值(全部為0), a的最后一列為常數(shù)項
  2. 開始計算梯度
  3. 獲取最優(yōu)的alphas的值以及對應的θ值
  4. 返回最終的θ值

alphas ：多個學習率
threshold：收斂的值
maxIter：最多迭代次數(shù)
addConstantItem：θ是否包含常數(shù)項

## 要求X必須是List集合，Y也必須是List集合
def fit(X, Y, alphas, threshold=1e-6, maxIter=200, addConstantItem=True):
    import math
    import numpy as np
    ## 校驗
    validate(X, Y)
    ## 開始模型構建
    l = len(alphas)
    m = len(Y)
    n = len(X[0]) + 1 if addConstantItem else len(X[0])#樣本的個數(shù)
    B = [True for i in range(l)]#模型的格式：控制最優(yōu)模型
    
    ## 差異性(損失值)
    ## J是不同的alphas學習率(步長)下的損失值
    J = [np.nan for i in range(l)]#loss函數(shù)的值
    
    # 1. 隨機初始化0值(全部為0), a的最后一列為常數(shù)項
    ## a是不同的alphas學習率(步長)下的定義個多個[0,0, ... ,0]
    ## 比如有3個步長參數(shù)，那么 a= [[0,0, ... ,0]，[0,0, ... ,0]，[0,0, ... ,0]]
    a = [[0 for j in range(n)] for i in range(l)]#theta，是模型的系數(shù)
    
    # 2. 開始計算
    # 迭代的次數(shù)
    for times in range(maxIter):
        # l是alphas學習率的長度，根據(jù)每個學習率依次迭代
        for i in range(l):
            if not B[i]:
                # 如果當前alpha的值已經(jīng)計算到最優(yōu)解了，那么不進行繼續(xù)計算
                continue
            
            # 初始a=[[0,0, ... ,0]，[0,0, ... ,0]，[0,0, ... ,0]] 
            # a[i] = [0,0, ... ,0]
            # ta[i] 即系數(shù)θi  
            ta = a[i]
            # 更新每一個θ的取值
            for j in range(n):
                # 從第i個學習率開始更新
                alpha = alphas[i]
                # 設初始值為0
                ts = 0
                # m對應Y的長度，m對應第幾個觀測值
                for k in range(m):
                    # 如果增加常數(shù)項為 TRUE
                    if j == n - 1 and addConstantItem:
                        # calcDiffe 計算差異，最終的體現(xiàn)是 真實值-預測值
                        # 不斷的迭代是一個累加的過程
                        # 由于常數(shù)項為TRUE，[h(x)^(i) -y^(i)]xj 中的xj為1
                        
                        # X[k]=[1.17640523]  Y[k][0]= -18.55967185 
                        ts += alpha*calcDiffe(X[k], Y[k][0], a[i]) * 1
                    else:
                        # 由于常數(shù)項為FALSE，[h(x)^(i) -y^(i)]xj 
                        ts =ts+alpha*calcDiffe(X[k], Y[k][0], a[i]) * X[k][j]
                
                # 初始a=[[0,0, ... ,0]，[0,0, ... ,0]，[0,0, ... ,0]] 
                # a[i] = [0,0, ... ,0]
                # ta[i] 即系數(shù)θi  
                # ts:梯度
                t = ta[j] + ts
                ta[j] = t
            ## 計算完一個alpha值的0的損失函數(shù)
            flag = True
            js = 0
            for k in range(m):
                js += math.pow(calcDiffe(X[k], Y[k][0], a[i]),2)+a[i][j]
                if js > J[i]:
                    flag = False
                    break;
            if flag:
                J[i] = js
                for j in range(n):
                    a[i][j] = ta[j]
            else:
                # 標記當前alpha的值不需要再計算了
                B[i] = False        
        ## 計算完一個迭代，當目標函數(shù)/損失函數(shù)值有一個小于threshold的結束循環(huán)
        r = [0 for j in J if j <= threshold]
        if len(r) > 0:
            break
        # 如果全部alphas的值都結算到最后解了，那么不進行繼續(xù)計算
        r = [0 for b in B if not b]
        if len(r) > 0:
            break;
    # 3. 獲取最優(yōu)的alphas的值以及對應的0值
    min_a = a[0]
    min_j = J[0]
    min_alpha = alphas[0]
    for i in range(l):
        if J[i] < min_j:
            min_j = J[i]
            min_a = a[i]
            min_alpha = alphas[i]
    
    print("最優(yōu)的alpha值為:",min_alpha)
    
    # 4. 返回最終的0值
    return min_a

4、 輸入樣本X和參數(shù)θ，并預測結果Y

def predict(X,a):
    Y = []
    n = len(a) - 1
    for x in X:
        result = 0
        for i in range(n):
            result += x[i] * a[i]
        result += a[n]
        Y.append(result)
    return Y

5、 計算實際值和預測值之間的相關性

def calcRScore(y,py):
    if len(y) != len(py):
        raise Exception("參數(shù)異常")
    import math 
    import numpy as np
    avgy = np.average(y)
    m = len(y)
    rss = 0.0
    tss = 0
    for i in range(m):
        rss += math.pow(y[i] - py[i], 2)
        tss += math.pow(y[i] - avgy, 2)
    r = 1.0 - 1.0 * rss / tss
    return r
測試自己寫的梯度下降：
import numpy as np
import matplotlib as mpl
import matplotlib.pyplot as plt
import pandas as pd
import warnings
import sklearn
from sklearn.linear_model import LinearRegression,Ridge, LassoCV, RidgeCV, ElasticNetCV
from sklearn.preprocessing import PolynomialFeatures
from sklearn.pipeline import Pipeline
from sklearn.linear_model.coordinate_descent import ConvergenceWarning

# inline 在行內(nèi)顯示
# plt.show() 在行內(nèi)顯示
%matplotlib inline 

## 設置字符集，防止中文亂碼
mpl.rcParams['font.sans-serif']=[u'simHei']
mpl.rcParams['axes.unicode_minus']=False

# warnings.filterwarnings(action = 'ignore', category=ConvergenceWarning)
## 創(chuàng)建模擬數(shù)據(jù)
np.random.seed(0)
np.set_printoptions(linewidth=1000, suppress=True)
N = 10
x = np.linspace(0, 6, N) + np.random.randn(N)
y = 1.8*x**3 + x**2 - 14*x - 7 + np.random.randn(N)
x.shape = -1, 1
y.shape = -1, 1
x

array([[ 1.76405235],
[ 1.06682388],
[ 2.31207132],
[ 4.2408932 ],
[ 4.53422466],
[ 2.35605545],
[ 4.95008842],
[ 4.51530946],
[ 5.23011448],
[ 6.4105985 ]])

plt.figure(figsize=(12,6), facecolor='w')

## 模擬數(shù)據(jù)產(chǎn)生
x_hat = np.linspace(x.min(), x.max(), num=100)
x_hat.shape = -1,1

## 線性模型
model = LinearRegression()
model.fit(x,y)
y_hat = model.predict(x_hat)
s1 = calcRScore(y, model.predict(x))
print(model.score(x,y)) ## 自帶R^2輸出
print ("模塊自帶實現(xiàn)===============")
print ("參數(shù)列表:", model.coef_)
print ("截距:", model.intercept_)

## 自模型
ma = fit(x,y,np.logspace(-4,-2,100), addConstantItem=True)
y_hat2 = predict(x_hat, ma)
s2 = calcRScore(y, predict(x,ma))
print ("自定義實現(xiàn)模型=============")
print ("參數(shù)列表:", ma)

## 開始畫圖
plt.plot(x, y, 'ro', ms=10, zorder=3)
plt.plot(x_hat, y_hat, color='#b624db', 
  lw=2, alpha=0.75, label=u'Python模型，
    $R^2$:%.3f' % s1, zorder=2)

plt.plot(x_hat, y_hat2, color='#6d49b6', 
   lw=2, alpha=0.75, label=u'自己實現(xiàn)模型，
    $R^2$:%.3f' % s2, zorder=1)

plt.legend(loc = 'upper left')
plt.grid(True)
plt.xlabel('X', fontsize=16)
plt.ylabel('Y', fontsize=16)

plt.suptitle(u'自定義的線性模型和模塊中的線性模型比較', fontsize=22)
plt.show()

0.837437698825
模塊自帶實現(xiàn)===============
參數(shù)列表: [[ 72.0576022]]
截距: [-163.71132966]
最優(yōu)的alpha值為: 0.01
自定義實現(xiàn)模型=============
參數(shù)列表: [70.879363936338876, -158.49974583659909]


結論：比較LinearRegression()和自己寫的梯度下降求θ的模型，最后根據(jù)兩種模型對測試集進行預測，我們發(fā)現(xiàn)兩條模型幾乎重合。說明自己寫的模型還是不錯的。
最后看看sklearn自帶的梯度下降法：
from sklearn.ensemble import GradientBoostingRegressor
clf = GradientBoostingRegressor()
y1 = y.ravel()
clf.fit(x,y1)
print ("自帶梯度下降法R方:", clf.score(x,y1))
y_hat3=clf.predict(x_hat)
s3=calcRScore(y, clf.predict(x))

## 開始畫圖
plt.plot(x, y, 'ro', ms=10, zorder=3)

plt.plot(x_hat, y_hat, color='#b624db', 
  lw=2, alpha=0.75, label=u'Python模型，
    $R^2$:%.3f' % s1, zorder=2)

plt.plot(x_hat, y_hat2, color='#6d49b6', 
  lw=2, alpha=0.75, label=u'自己實現(xiàn)模型，
    $R^2$:%.3f' % s2, zorder=1)

plt.plot(x_hat, y_hat3, color='#6daaba', 
  lw=2, alpha=0.75, label=u'自帶梯度下降方法，
    $R^2$:%.3f' % s3, zorder=1)
plt.legend(loc = 'upper left')
plt.grid(True)
plt.xlabel('X', fontsize=16)
plt.ylabel('Y', fontsize=16)

plt.suptitle(u'自定義的線性模型和模塊中的線性模型比較', fontsize=22)
plt.show()

自帶梯度下降法R方: 0.999999998927

最后編輯于 ：
?著作權歸作者所有,轉載或內(nèi)容合作請聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時請結合常識與多方信息審慎甄別。
平臺聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點，簡書系信息發(fā)布平臺，僅提供信息存儲服務。

相關閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容