機器學習之局部加權線性回歸、邏輯回歸演算法（附演算法程序）

05-22

來自專欄小山的知識庫

一、緒言

在看這篇文章之前，如果你對標題的兩個演算法沒有任何的接觸，我建議你可以先到網易公開課觀看對應的視頻（網址：lecture 3），看完後歡迎你回來批判性地閱讀該文章，而且文章中會附上課堂上每種演算法的相關項目程序，相信會讓你對課堂學到的知識有更深的理解。

接下來要講述的演算法分別是局部加權線性回歸演算法(locally weighted linear regression)，該演算法克服了梯度下降法（相關文章：機器學習之梯度下降法）不能擬合非線性曲線的缺點。而邏輯回歸(logistics regression)則不同於前兩種演算法，下面會講述如何用該演算法來解決分類的問題。

二、局部加權線性回歸演算法(LWLR)

在上一篇文章的線性回歸梯度下降法中，它的最終目的是求出使得代價函數

式2-1

達到最小值的線性擬合函數：

不過梯度下降法存在一個問題，它只能進行線性擬合，而不能進行非線性擬合，就像下圖的訓練庫：

圖2-1

不難發現，圖中的數據更適合用曲線去擬合，如果強行用梯度下降法對其進行擬合，效果就會顯得不好！好吧，可能現在說不定又忘了用梯度下降實現線性回歸的方法，在這先複習一下，有三種方法可以實現：

批量梯度下降法（batch gradient descent）：

隨機梯度下降法（stochastic gradient descent）：

矩陣實現形式：

如果想了解這三種演算法的推導過程，請參考相關文章：機器學習之梯度下降法。為了方便起見，我用梯度下降法的矩陣實現形式來對圖2-1進行線性擬合，相關程序如下：

#代碼1clc;clear;close all; %初始化load(datax.dat); %導入訓練庫的輸入變數dataxload(datay.dat); %導入訓練庫的輸入變數datayhold on;scatter(datax,datay); %畫出出訓練庫所有樣本對應的點X = [ones(size(datax,1),1) datax];Y = datay;theta = [];theta = inv(X*X)*X*Y; %對應梯度下降法的矩陣實現形式plot_x = -5:0.01:15;plot_y = theta(1) + theta(2) * plot_x;plot(plot_x,plot_y); %畫出擬合直線hold off

用梯度下降法進行線性擬合後，得到的結果如下：

圖2-2

一看就知道，這個擬合顯得太過勉強，是不合格的，我們得要對這個樣本庫進行非線性擬合才行。不過學過微積分的人都知道只要是可導函數，只要在合適的範圍下都能進行以導數為斜率的線性近似。同樣的思想，對於圖2-2中每個橫坐標x，我們只對其一定範圍內的樣本進行線性擬合，即非線性曲線的局部線性化，就可以實現非線性曲線的擬合，舉個例子，如下圖（舉了兩個點x1和x2，原諒畫得不好看=_=）：

圖2-3

那如何實現呢？不難發現線性回歸是在考慮全部樣本值後妥協得出的線性擬合，而局部線性回歸則只重視某一局部中的樣本值妥協得出的線性擬合，並且橫坐標不同，對應的線性擬合也會不同。

而線性回歸又是如何對全部樣本值進行考慮的呢？是通過對代價函數 $J( heta)$ （式2-1）進行偏導，來求取使 $J( heta)$ 達到最小值的一組線性參數 $heta$ 。因此根據線性回歸的原理，只需把梯度下降法的代價函數（式2-1）做點小動作，改成：

式2-2

其中 $J_{x_{j}}( heta)$ 代表是點 $x_{j}$ 的代價函數，每一個不同的橫坐標都有自己的代價函數，就如同每個點都有不同的斜率一樣，而 $W_{x_{j}}^{(i)}$ 代表是點 $x_{j}$ 的加權函數，該函數不固定，但必須要有的特徵為樣本橫坐標 $x^{(i)}$ 離 $x_{j}$ 越遠，其值越小，否則反之，通俗點說，就是離 $x_{j}$ 越近的樣本就越重視！然後我們就像梯度下降法一樣，對該點代價函數 $J_{x_{j}}( heta)$ 求偏導，對參數進行更新迭代，最後得到某一個點 $x_{j}$ 的局部線性回歸直線，就如圖2-3一樣，因此該方法稱為局部加權線性回歸。在很多情況下，我們會使加權函數 $W_{x_{j}}^{(i)}$ 取：

其中 $sigma$ 可以對局部線性擬合的範圍進行調節。類似第一篇文章梯度下降法的推導過程，根據式2-2，對於任意一點 $x_{j}$ ，我們可以得到實現LWLR的參數更新迭代方程：

式2-3

以及它的矩陣實現形式：

式2-4

其中

矩陣推理過程與上一篇文章梯度下降法的矩陣推理一樣，只不過在局部加權線性回歸中，某一點的代價函數矩陣形式發生了變化，變為 $J( heta)=frac{1}{2m}(X heta-Y)^{T}W(X heta-Y)$ 。

下面首先對LWLR的參數更新迭代方程（式2-3）進行程序的模擬：

註：該位置用前面的#代碼1代替#代碼2：LWLR參數更新迭代m = length(Y);LWLR_Xp = -5:0.01:13; %生成要進行局部線性回歸的橫坐標sigma = 0.8;alpha = 0.1;for i = 1:length(LWLR_Xp) theta = zeros(2,1); for j = 1:m W(j,j) = exp((-1/(2*(sigma^2))) * (LWLR_Xp(i)-X(j,2))^2); %求取加權矩陣 end for iter = 1:200 theta(1) = theta(1) - (alpha/m)*(X*theta-Y)*W*X(:,1); theta(2) = theta(2) - (alpha/m)*(X*theta-Y)*W*X(:,2); %參數的迭代更新 end LWLR_Yp(i) = theta(1) + theta(2)*LWLR_Xp(i); %根據生成對應的縱坐標endplot(LWLR_Xp,LWLR_Yp);

得到下面的結果：

感覺結果還行，不過用迭代的方法求取參數，首先在運行速度上會比較慢，其次還要選取合適的係數 $alpha$ 才能達到比較好的擬合效果，而矩陣實現形式（式2-4）可以很好地解決以上問題，其相關代碼如下：

註：該位置用前面的#代碼1代替#代碼3：LWLR矩陣實現形式sigma = 0.8;m = length(Y);LWLR_Xp = -5:0.01:13;for i = 1 : length(LWLR_Xp) for j = 1:m W(j,j) = exp((-1/(2*(sigma^2))) * (LWLR_Xp(i)-X(j,2))^2); end theta = inv(X*W*X)*X*W*Y; %對應文章式2-4 LWLR_Yp(i) = theta(1) + theta(2)*LWLR_Xp(i);endplot(LWLR_Xp,LWLR_Yp);

具體結果如下：

三、探究代價函數 $J( heta)$

到這裡，我覺得很有必要搞清楚為什麼線性回歸的代價函數為式2-1：

式2-1

因為在接下來講邏輯回歸時，兩者代價函數的求取思路具有一定的相似性，所以我們最好把線性回歸的代價函數求取思路給過一遍。假設樣本庫中的樣本相互獨立，且樣本為：

它們輸入與輸出的關係為：

式2-6

其中 $heta$ 為線性回歸參數集， $epsilon^{(i)}$ 為樣本輸出與線性回歸輸出的誤差，假設 $epsilon^{(i)}$ 服從高斯分布，即 $epsilon^{(i)}sim(0,sigma^{2})$ ，所以可得：

式2-7

根據式2-6，可以將式2-7變為：

因為樣本庫中的樣本相互獨立，因此它們聯合概率為：

式2-8

此時，為方便運算，對式2-8進行取對數，得到：

式2-9

觀察式2-9，當使式子：

式2-10

達到最小值時， $L( heta)$ 達到最大值，此時樣本庫輸出與線性回歸函數接近程度最大！接著再對式2-10取平均，我們就可以得到代價函數 $J( heta)$ 。通俗點來講，代價函數衡量採用擬合函數 $h_{ heta}(x)$ 估計時所要付出的代價程度，代價函數值越小，代表擬合之後要付出的代價就越小，擬合程度也就越好！

可能到這裡，大家可能會對其中的論證邏輯表示懷疑，說實話，我自己也有不小的疑惑，不過這個論證的套路在線性回歸和邏輯回歸中都適用，裡面肯定存在一定的合理性。對於為什麼可以假定誤差 $varepsilon^{(i)}$ 服從高斯分布？主要有兩個原因，第一，高斯分布便於數學化簡運算；第二，現實中大部分的誤差概率分布都服從高斯分布！

四、邏輯回歸演算法

在邏輯回歸中，輸出y是離散的，而輸入x可以是連續的。我們先集中討論輸出只有0,1兩種情況。對於輸出只有0,1，我們可以用下面函數對樣本庫進行擬合：

式4-1

該函數圖像為：

對該函數求導，得到：

式4-2

現在我們開始假設輸出概率分布函數為：

式4-3

結合式4-3，我們可以得到一般輸入與輸出概率函數：

式4-4

與上面第三小節推理思路相似，根據式4-4，我們來求取存有m個樣本的訓練庫的聯合概率對數形式：

因此，類似線性回歸代價函數的推導過程，我們便可以得到此次邏輯回歸的代價函數為：

式4-5

於是，我們採用梯度下降法的思想，求取代價函數關於參數集的偏導，以此求出最優參數：

式4-6

根據式4-2，可以把式4-6化簡為：

式4-7

根據式4-7，我們可以得到參數迭代更新方程：

式4-8

不過在這裡要告訴大家一個很遺憾的消息，我對迭代方程式4-8在matlab上進行模擬，很遺憾未能找到對應 $J( heta)$ 最小值的那一組參數 $heta$ ，同時也在網上找了好多資料，未發現用迭代方程求取最佳參數的程序，對此，我給出迭代方程程序模擬失敗的兩個原因：1、有可能我的程序問題，但我沒發現。2、可能由於代價函數 $J( heta)$ 非線性的原因，存在多個極值點，普通的參數迭代演算法找不到通往最小值的道路！如果有用式4-8找出最佳參數的網友，望請賜教。

雖然用在程序上式4-8找不到使得代價函數 $J( heta)$ 最小的參數集，但是得益於Matlab的強大運算功能，fminunc函數可以幫你找到，至於它內部演算法就沒有去了解過，反正它可以幫你找到使目標函數達到最小值的點！具體函數的使用請參考Find minimum of unconstrained multivariable function。

下面要進行邏輯回歸的背景是某一個老師要通過考察學生兩次考試的成績來決定是否通過考核，於是她拿出以前學生的一些考核數據作為參考，如下圖：

相關顯示代碼如下：

#代碼4：數據顯示代碼clc;clear; close all;data = load(ex2data1.txt);X = data(:,1:2);Y = data(:,3);hold on;pos = find(Y == 1); %y為1代表通過審核neg = find(Y == 0); %y為0代表通過審核plot(X(pos,1),X(pos,2),b+); plot(X(neg,1),X(neg,2),ro);xlabel(第一次考試成績);ylabel(第一次考試成績);legend(通過, 未通過);hold off;

下面要使用函數fminunc找出使得代價函數最小的參數集 $heta$ ，不過在這之前，你必須為代價函數建立一個函數文件，要注意的是這個函數文件要有兩個輸出，第一個輸出是代價函數值，第二是代價函數的梯度，否則等下調用fminunc函數會出錯，如下代碼：

#代碼5：cost_functions.m文件function [J, grad] = cost_functions(f_theta,Xp,Yp) %兩個輸出，代價函數以及梯度m = length(Yp);sigmoid_vector = 1./(1+exp(-(Xp*f_theta))); %對應式4-3J = (-1/m) * (Yp*(log(sigmoid_vector)) + (1-Yp)*log(1-sigmoid_vector));%求取代價函數grad = (1/m) * (sigmoid_vector-Yp)*Xp; %求取代價函數的梯度end

建立好代價函數 $J( heta)$ 的程序函數後，開始調用fminunc函數，相關代碼如下：

#代碼6：調用fminunc函數X = [ones(size(X,1),1) X];theta = zeros(size(X,2),1);options = optimset(GradObj, on, MaxIter, 5000); %模式設置[theta,cost] = fminunc(@(t)(cost_functions(t,X,Y)),theta,options); %調用fminunc函數fprintf(theta(1) is %f theta(2) is %f theta(3) is %f the cost is %f ,theta(1),theta(2),theta(3),cost);

分析一下代碼：GradObj表示告訴fminunc，cost_functions函數將會返回代價函數以及它的梯度，然後fminunc就會使用得到的梯度值來尋找代價函數最小值的點。MaxIter則代表fminunc運算過程中迭代次數上限值。而@(t)(cost_functions(t,X,Y))則代表cost_functions輸入參數中的第一個參數為變數，其餘為固定常數。而參考#代碼5，cost_functions函數的第一個參數便是 $heta$ ，也就是找出使得代價函數為最小值的那一組 $heta$ 。運行結果為：

最優參數求出來後，後面就到如何在數據點根據參數來畫分界線。在一開始的推理中，我們就假設 $h_{ heta}(x)$ 作為概率分布函數，如式4-3。因此，當 $h_{ heta}( heta^{T}x)>frac{1}{2}$ 時，我們可以判定該學生通過考核，否則反之。因此我們要找出使得 $h_{ heta}( heta^{T}x)=frac{1}{2}$ 的邊界值，根據式4-1可得當x處於邊界值時，有：

根據該關係式，便可以求出第一成績x1與第二次成績x2的函數關係，進而畫出分界線。下面是相關代碼：

#代碼7：畫出分界線hold on;plot_x = [min(X(:,2)),max(X(:,2))];plot_y = (-1/theta(3)) * (theta(1) + theta(2)*plot_x);plot(plot_x,plot_y);

相關運行結果如下：

從圖中可以看到，邏輯回歸演算法對數據進行比較正確的一個分類。

五、結語

本次文章在梯度下降法（相關文章：機器學習演算法之梯度下降（Gradient descent）附相關項目）的基礎講述了局部加權線性回歸演算法，並對其中的參數迭代更新方程以及矩陣實現形式都進行了matlab程序的模擬，並且按照類似線性回歸的推理思路，對邏輯回歸演算法進行了一番推導，並在matlab程序上使用fminunc函數實現了邏輯回歸演算法，如有錯誤，望請賜教，謝謝！

關於本次文章的項目程序以及訓練數據已經上傳至網盤：文章二相關程序，不過由於在文章中已對程序進行分析，所以上傳的matlab程序不會再有注釋！