線性回歸的推導

什么是線性回歸?

回歸算法是一種有監督算法

回歸算法是一種比較常用的機器學習算法，用於構建一個模型來做特征向量到標簽的映射。，用來建立“解釋”變量(自變量X)和觀測值(因變量Y)之間的關系。在算法的學習過程中，試圖尋找一個模型，最大程度擬合訓練數據。

回歸算法在使用時，接收一個n維度特征向量，輸出一個連續的數據值。

下面通過一則案例引出線性回歸:

房屋面積(m^2)	租賃價格(1000￥)
10	0.8
15	1
20	1.8
30	2
50	3.2
60	3
60	3.1
70	3.5

目的:假如有一個新的房屋面積比如35平方米,那我們應該怎么預測價格?

假設

​ 輸入：x特征向量

​ 輸出：hθ(x)即預測值

假設擬合函數:

 

解決方案:

最小二乘法（又稱最小平方法）是一種數學優化技術。它由兩部分組成：

​ 一、計算所有樣本誤差的平均（代價函數）

​ 二、使用最優化方法尋找數據的最佳函數匹配（抽象的）

最小二乘法是抽象的，具體的最優化方法有很多，比如正規方程法、梯度下降法、牛頓法、擬牛頓法等等。

 

線性回歸的表達式

選擇最優的θ值構成算法公式,並最終要求是計算出θ的值:

 

預備知識

大數定理

意義

隨着樣本容量n的增加，樣本平均數將接近於總體平均數(期望μ)，所以在統計推斷中，一般都會使用樣本平均數估計總體平均數的值。也就是我們會使用一部分樣本的平均值來代替整體樣本的期望/均值，出現偏差的可能是存在的，但是當n足夠大的時候，偏差的可能性是非常小的，當n無限大的時候，這種可能性的概率基本為0。

作用

就是為使用頻率來估計概率提供了理論支持；為使用部分數據來近似的模擬構建全部數據的特征提供了理論支持。

如下圖的例子:

從1到10中進行隨機抽取數字,經過非常多次循環后,我們發現均值趨向5

import matplotlib.pyplot as plt
import random
import numpy as np
​
random.seed(28)
​
​
def generate_random_int(n):
    return [random.randint(1, 9) for i in range(n)]
​
​
if __name__ == '__main__':
    number = 8000
    x = [i for i in range(number + 1) if i != 0]
    total_random_int = generate_random_int(number)
    y = [np.mean(total_random_int[0:i + 1]) for i in range(number)]
    plt.plot(x, y, "b-")
    plt.xlim(0, number)
    plt.grid(True)
    plt.show()

 

 

中心極限定理

意義

中心極限定理就是一般在同分布的情況下，抽樣樣本值的規范和在總體數量趨於無窮時的極限分布近似於正態分布。

作用

假設{Xn}為獨立同分布的隨機變量序列，並具有相同的期望μ和方差為σ2，則{Xn}服從中心極限定理，且Zn為隨機序列{Xn}的規范和：

 

 

隨機的拋六面的骰子，計算三次的點數的和, 三次點數的和其實就是一個事件A，現在問事件A發生的概率以及事件A所屬的分布是什么?

import matplotlib.pyplot as plt
import random
import numpy as np
​
random.seed(28)
def generate_random_int():
    # 模擬產生篩子
    return random.randint(1,6)
​
def generate_mean(n):
    # 隨機選擇三個篩子的和
    return np.sum([generate_random_int() for i in range(n)])
​
if __name__ == '__main__':
    number1 = 1000000
    number2 = 3
    keys = [generate_mean(number2) for i in range(number1)]
    result = {}
    for key in keys:
        count = 1
        if key in result:
            count+=result[key]
        result[key] = count
    x = sorted(np.unique(list(result.keys())))
    y = []
    for key in x:
        y.append(result[key]/number1)
    plt.plot(x,y,"b-")
    plt.xlim(x[0]-1,x[-1]+1)
    plt.grid(True)
    plt.show()

 

 

最大似然函數

最大概似估計、極大似然估計，是一種具有理論性的參數估計方法。基本思想是：當從模型總體隨機抽取n組樣本觀測值后，最合理的參數估計量應該使得從模型中抽取該n組樣本觀測值的概率最大；一般步驟如下：

1. 寫出似然函數；

2. 對似然函數取對數，並整理；

3. 求導數；

4. 解似然方程

梯度下降

優化思想

用當前位置負梯度方向作為搜索方向，因為該方向為當前位置的最快下降方向，所以梯度下降法也被稱為“最速下降法”。梯度下降法中越接近目標值，變量變化越小。計算公式如下：

 

α被稱為步長或者學習率(learningrate)，表示自變量θ每次迭代變化的大小。

收斂條件：當目標函數的函數值變化非常小的時候或者達到最大迭代次數的時候，就結束循環。

分類

批量梯度下降法（BatchGradient Descent）

•有N個樣本，求梯度的時候就用了N個樣本的梯度數據

•優點：准確 缺點：速度慢

 

隨機梯度下降法（StochasticGradient Descent，SGD）

•和批量梯度下降法原理類似，區別在於求梯度時沒有用所有的N個樣本的數據，而是僅僅選取一個樣本來求梯度

•優點：速度快 缺點：准確度低

 

小批量梯度下降法（Mini-batchGradient Descent）

•批量梯度下降法和隨機梯度下降法的折衷，比如N=100

•spark中使用的此方法

 

例子

import numpy as np
import matplotlib.pyplot as plt
def f(x):
    return 0.5*(x-0.25)**2
def h(x):
    return x-1/4
​
X = []
Y = []
x =2
# 此時的學習率,是由一定區間的,過小會導致性能變慢,過大可能發散
step = 0.8
f_change = f(x)
f_current = f(x)
X.append(x)
Y.append(f_current)
# 當變化值小於1e-10的時候停止,也就是梯度
while f_change > 1e-10 and len(X)<100:
    # 這里是梯度下降的變化程度x = x - k*f(x)'
    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)
fig = plt.figure()
X2 = np.arange(-2.1,2.65,0.05)
Y2 = 0.5*(X2-0.25)**2
plt.plot(X2,Y2,'-',color="#666666",linewidth = 2)
plt.plot(X,Y,'bo--')
plt.title("$y= 0.5*(x-0.25)^2$通過梯度下降法得到目標值x=%.2f,y=%.2f迭代次數%d"%(x,f_current,len(X)))
plt.show()

 

 

三維數據中的梯度下降

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
def f(x, y):
    return x ** 2 + y ** 2
​
​
def h(t):
    return 2 * t
​
​
X = []
Y = []
Z = []
​
x = 2
y = 2
f_change = x ** 2 + y ** 2
f_current = f(x, y)
step = 0.1
X.append(x)
Y.append(y)
Z.append(f_current)
while f_change > 1e-10:
    # 對於各自未知數求偏導
    x = x - step * h(x)
    y = y - step * h(y)
    f_change = np.abs(f_current - f(x,y))
    f_current = f(x,y)
    X.append(x)
    Y.append(y)
    Z.append(f_current)
fig = plt.figure()
ax = Axes3D(fig)
X2 = np.arange(-2,2,0.2)
Y2 = np.arange(-2,2,0.2)
X2,Y2 = np.meshgrid(X2,Y2)
Z2 = X2**2 + Y2**2
ax.plot_surface(X2,Y2,Z2,rstride=1,cstride=1,cmap='rainbow')
ax.plot(X,Y,Z,'ro--')
​
ax.set_title("")
plt.show()

 

 

推導線性回歸

確定目標函數:

 

 

首先我們找到的數據,或者在測量的過程中肯定會產生或多或少的數據,所以這里就利用中心極限定理,我們可以將誤差epsilon認為是獨立分布的服從均值為0方差為σ^2的高斯分布。

 

 

利用最大似然函數來求解:

目標函數為:

 

 

因為累乘求導比較復雜,這里利用對數化,形成對數似然函數:

 

 

這里我們主要是求似然函數的最大值,前面是常量忽略不計,所以主要是后面的內容。

 

 

這里的目標函數是平方和損失函數,除了它還有0-1損失函數,感知損失函數,絕對值損失函數,對數損失函數。

利用梯度下降對目標函數θ求解

 

 

初始化θ(隨機初始化，可以初始為0)沿着負梯度方向迭代，更新后的θ使J(θ)更小

對損失函數進行求導:

 

 

實現代碼:

批量下降線性回歸:

import numpy as np
​
# 構造訓練數據集
x_train = np.array([[2, 0., 3], [3, 1., 3], [0, 2., 3], [4, 3., 2], [1, 4., 4]])
m = len(x_train)
print(m)
​
x0 = np.full((m, 1), 1)
print(x0)
​
# 構造一個每個數據第一維特征都是1的矩陣
input_data = np.hstack([x0, x_train])
m, n = input_data.shape
theta1 = np.array([[2, 3, 4]]).T
# 構建標簽數據集，后面的np.random.randn是將數據加一點噪聲，以便模擬數據集。
# y_train = (input_data.dot(np.array([1, 2, 3, 4]).T)).T
y_train = x_train.dot(theta1) + np.array([[2], [2], [2], [2], [2]])
​
# 設置兩個終止條件
loop_max = 1000000
epsilon = 1e-5
​
# 初始theta
np.random.seed(0)  # 設置隨機種子
theta = np.random.randn(n, 1)  # 隨機取一個1維列向量初始化theta
​
# 初始化步長/學習率
alpha = 0.00001
# 初始化誤差，每個維度的theta都應該有一個誤差，所以誤差是一個4維。
error = np.zeros((n, 1))  # 列向量
​
# 初始化偏導數
diff = np.zeros((input_data.shape[1], 1))
​
# 初始化循環次數
count = 0
​
while count < loop_max:
    count += 1
    sum_m = np.zeros((n, 1))
​
    for i in range(m):
        for j in range(input_data.shape[1]):
            # 計算每個維度的theta
            diff[j] = (input_data[i].dot(theta) - y_train[i]) * input_data[i, j]
        # 求每個維度的梯度的累加和
        sum_m = sum_m + diff
    # 利用這個累加和更新梯度
    theta = theta - alpha * sum_m
    # else中將前一個theta賦值給error,theta - error便表示前后兩個梯度的變化，當梯度
    # 變化很小（在接收的范圍內）時，便停止迭代。
    if np.linalg.norm(theta - error) < epsilon:
        break
    else:
        error = theta
​
print(theta)

 

隨機梯度線性回歸:

# -*- coding: utf-8 -*-
"""
Created on Wed Nov 28 14:02:17 2018
​
@author: Administrator
"""
​
import numpy as np
​
# 構造訓練數據集
x_train = np.array([[2, 0., 3], [3, 1., 3], [0, 2., 3], [4, 3., 2], [1, 4., 4]])
​
# 構建一個權重作為數據集的真正的權重，theta1主要是用來構建y_train，然后通過模型計算
# 擬合的theta，這樣可以比較兩者之間的差異，驗證模型。
theta1 = np.array([[2 ,3, 4]]).T
​
# 構建標簽數據集,y=t1*x1+t2*x2+t3*x3+b即y=向量x_train乘向量theta+b, 這里b=2
y_train = (x_train.dot(theta1) + np.array([[2],[2],[2],[2],[2]])).ravel()
# 構建一個5行1列的單位矩陣x0，然它和x_train組合，形成[x0, x1, x2, x3]，x0=1的數據形式，
# 這樣可以將y=t1*x1+t2*x2+t3*x3+b寫為y=b*x0+t1*x1+t2*x2+t3*x3即y=向量x_train乘向
# 量theta其中theta應該為[b, *, * , *]，則要擬合的theta應該是[2,2,3,4]，這個值可以
# 和算出來的theta相比較，看模型的是否達到預期
x0 = np.ones((5, 1))
input_data = np.hstack([x0, x_train])
m, n = input_data.shape
​
# 設置兩個終止條件
loop_max = 10000000
epsilon = 1e-6
​
# 初始化theta（權重）
np.random.seed(0)
theta = np.random.rand(n).T # 隨機生成10以內的，n維1列的矩陣
​
# 初始化步長/學習率
alpha = 0.000001
# 初始化迭代誤差（用於計算梯度兩次迭代的差）
error = np.zeros(n)
​
# 初始化偏導數矩陣
diff = np.zeros(n)
​
# 初始化循環次數
count = 0
​
while count < loop_max:
    count += 1  # 沒運行一次count加1，以此來總共記錄運行的次數
    # 計算梯度
    for i in range(m):
        # 計算每個維度theta的梯度，並運用一個梯度去更新它
        diff = input_data[i].dot(theta)-y_train[i]
        theta = theta - alpha * diff*(input_data[i])
    
    # else中將前一個theta賦值給error,theta - error便表示前后兩個梯度的變化，當梯度
    #變化很小（在接收的范圍內）時，便停止迭代。
    if np.linalg.norm(theta - error) < epsilon: # 判斷theta與零向量的距離是否在誤差內
        break
    else:
        error = theta  
print(theta)

 

​

小批量梯度下降

import numpy as np
​
# 構造訓練數據集
x_train = np.array([[2, 0., 3], [3, 1., 3], [0, 2., 3], [4, 3., 2], [1, 4., 4]])
m = len(x_train)
​
x0 = np.full((m, 1), 1)
​
# 構造一個每個數據第一維特征都是1的矩陣
input_data = np.hstack([x0, x_train])
m, n = input_data.shape
theta1 = np.array([[2, 3, 4]]).T
# 構建標簽數據集，后面的np.random.randn是將數據加一點噪聲，以便模擬數據集。
# y_train = (input_data.dot(np.array([1, 2, 3, 4]).T)).T
y_train = x_train.dot(theta1) + np.array([[2], [2], [2], [2], [2]])
​
# 設置兩個終止條件
loop_max = 1000000
epsilon = 1e-5
​
# 初始theta
np.random.seed(0)  # 設置隨機種子
theta = np.random.randn(n, 1)  # 隨機取一個1維列向量初始化theta
​
# 初始化步長/學習率
alpha = 0.00001
# 初始化誤差，每個維度的theta都應該有一個誤差，所以誤差是一個4維。
error = np.zeros((n, 1))  # 列向量
​
# 初始化偏導數
diff = np.zeros((input_data.shape[1], 1))
​
# 初始化循環次數
count = 0
​
# 設置小批量的樣本數
minibatch_size = 2
​
while count < loop_max:
    count += 1
    sum_m = np.zeros((n, 1))
    for i in range(1, m, minibatch_size):
        # for j in range(i - 1, i + minibatch_size - 1, 1):
        #     # 計算每個維度的theta
        #     diff[j] = (input_data[i].dot(theta) - y_train[i]) * input_data[i, j]
        diff = (1/minibatch_size)*input_data[i:i+minibatch_size,:].T.dot(input_data[i:i+minibatch_size,:].dot(theta) - y_train[i:i+minibatch_size,:])
        # 求每個維度的梯度的累加和
        sum_m = sum_m + diff
    # 利用這個累加和更新梯度
    theta = theta - alpha * (1.0 / minibatch_size) * sum_m
    # else中將前一個theta賦值給error,theta - error便表示前后兩個梯度的變化，當梯度
    # 變化很小（在接收的范圍內）時，便停止迭代。
    if np.linalg.norm(theta - error) < epsilon:
        break
    else:
        error = theta
# 和sklearn庫相對比
# from sklearn.linear_model import LinearRegression
# model = LinearRegression(fit_intercept=False)
# model.fit(input_data,y_train)
# print(model.coef_)
print(theta)

 

線性回歸和正則化

向量范數

1-范數:即向量元素絕對值之和，調用函數np.linalg.norm(x, 1)

 

 

2-范數:Euclid范數（歐幾里得范數，常用計算向量長度），即向量元素絕對值的平方 和再開方，調用函數np.linalg.norm(x) (默認為2)

 

 

無窮大范數:即所有向量元素絕對值中的最大值，調用函數np.linalg.norm (x, np.inf)

 

 

無窮小范數:即所有向量元素絕對值中的最小值，調用函數np.linalg.norm (x, -np.inf)

 

 

p-范數即向量元素絕對值的p次方和的1/p次冪，調用函數np.linalg.norm (x, p)

 

 

下圖表示了p從無窮到0變化時，三維空間中到原點的距離（范數）為1的點構成的圖形的變化情況。

 

經驗風險和結構風險(可跳過)

期望風險表示的是全局的概念，表示的是決策函數對所有的樣本(X,Y)預測能力的大小，而經驗風險則是局部的概念，僅僅表示決策函數對訓練數據集里樣本的預測能力。理想的模型（決策）函數應該是讓所有的樣本的損失函數最小的（也即期望風險最小化），但是期望風險函數往往是不可得到的，即上式中，X與Y的聯合分布函數不容易得到。現在我們已經清楚了期望風險是全局的，理想情況下應該是讓期望風險最小化，但是呢，期望風險函數又不是那么容易得到的。怎么辦呢？那就用局部最優的代替全局最優這個思想吧。這就是經驗風險最小化的理論基礎。

通過上面的分析可以知道，經驗風險與期望風險之間的聯系與區別。現在再總結一下：

經驗風險是局部的，基於訓練集所有樣本點損失函數最小化的。

期望風險是全局的，是基於所有樣本點的損失函數最小化的。

經驗風險函數是現實的，可求的；

期望風險函數是理想化的，不可求的；

只考慮經驗風險的話，會出現過擬合的現象，過擬合的極端情況便是模型f(x)對訓練集中所有的樣本點都有最好的預測能力，但是對於非訓練集中的樣本數據，模型的預測能力非常不好。怎么辦呢？這個時候就引出了結構風險。結構風險是對經驗風險和期望風險的折中。在經驗風險函數后面加一個正則化項（懲罰項）便是結構風險了。如下式：

 

相比於經驗風險，結構風險多了一個懲罰項，其中是一個l是一個大於0的系數。J(f)表示的是模型f的復雜度。結構風險可以這么理解：

經驗風險越小，模型決策函數越復雜，其包含的參數越多，當經驗風險函數小到一定程度就出現了過擬合現象。也可以理解為模型決策函數的復雜程度是過擬合的必要條件，那么我們要想防止過擬合現象的發生，就要破壞這個必要條件，即降低決策函數的復雜度。也即，讓懲罰項J(f)最小化，現在出現兩個需要最小化的函數了。我們需要同時保證經驗風險函數和模型決策函數的復雜度都達到最小化，一個簡單的辦法把兩個式子融合成一個式子得到結構風險函數然后對這個結構風險函數進行最小化。

 

簡而言之,經驗風險最小化可以理解為：最小化代價函數結構風險最小化可以理解為：最小化目標函數（代價函數+正則化項）

為什么要引入正則項

 

為了防止數據過擬合，也就是的θ值在樣本空間中不能過大/過小，可以在目標函數之上增加一個平方和損失：

 

Ridge回歸

L2正則化：權值向量w中各個元素的平方和

上面我們一直使用θ的平方作為正則項也就是L2正則,它的學名叫做Ridge回歸或者嶺回歸,

Lasso回歸

L1正則化：權值向量w中各個元素的絕對值之和

L1正則化 VS L2正則化

L1正則化可以產生稀疏權值矩陣，即產生一個稀疏模型，可以用於特征選擇。

L2正則化可以防止模型過擬合（overfitting）

為什么 L1 正則可以產生稀疏模型（很多參數=0），而L2 正則不會出現很多參數為0的情況？

 

由圖可知左面式L1正則，右面是L2正則，這里左面更容易和超平面在y軸上形成最優解，此時x=0也就是產生了稀疏解，這正是我們大多說時候想要的。

Elasitc Net

同時使用L1正則和L2正則的線性回歸模型就稱為Elasitc Net算法(彈性網絡算法)

 

到此線性回歸以及優化已經全部總結，訓練結束后可以使用MSE，RMSE，R^2這些指標進行評估。