邏輯回歸模型（Logistic Regression）及Python實現

邏輯回歸模型（Logistic Regression）及Python實現

http://www.cnblogs.com/sumai

1.模型

　　在分類問題中，比如判斷郵件是否為垃圾郵件，判斷腫瘤是否為陽性，目標變量是離散的，只有兩種取值，通常會編碼為0和1。假設我們有一個特征X，畫出散點圖，結果如下所示。這時候如果我們用線性回歸去擬合一條直線：h_θ(X) = θ₀+θ₁X，若Y≥0.5則判斷為1，否則為0。這樣我們也可以構建出一個模型去進行分類，但是會存在很多的缺點，比如穩健性差、准確率低。而邏輯回歸對於這樣的問題會更加合適。

     

　　邏輯回歸假設函數如下，它對θ^TX作了一個函數g變換，映射至0到1的范圍之內，而函數g稱為sigmoid function或者logistic function，函數圖像如下圖所示。當我們輸入特征，得到的h_θ(x)其實是這個樣本屬於1這個分類的概率值。也就是說，邏輯回歸是用來得到樣本屬於某個分類的概率。

       

      

 

 

2.評價

   回想起之前線性回歸中所用到的損失函數：

   

　如果在邏輯回歸中也運用這種損失函數，得到的函數J是一個非凸函數，存在多個局部最小值，很難進行求解，因此需要換一個cost函數。重新定義個cost函數如下：

   

   當實際樣本屬於1類別時，如果預測概率也為1，那么損失為0，預測正確。相反，如果預測為0，那么損失將是無窮大。這樣構造的損失函數是合理的，並且它還是一個凸函數，十分方便求得參數θ，使得損失函數J達到最小。

   

 

3.優化

    我們已經定義好了損失函數J(θ)，接下來的任務就是求出參數θ。我們的目標很明確，就是找到一組θ，使得我們的損失函數J(θ)最小。最常用的求解方法有兩種：批量梯度下降法(batch gradient descent)， 牛頓迭代方法((Newton's method)。兩種方法都是通過迭代求得的數值解，但是牛頓迭代方法的收斂速度更加快。

   批量梯度下降法： 

   牛頓迭代方法：  (H為海瑟矩陣)

 

4.python代碼實現

# -*- coding: utf-8 -*-
"""
Created on Wed Feb 24 11:04:11 2016

@author: SumaiWong
"""

import numpy as np
import pandas as pd
from numpy import dot
from numpy.linalg import inv

iris = pd.read_csv('D:\iris.csv')
dummy = pd.get_dummies(iris['Species']) # 對Species生成啞變量
iris = pd.concat([iris, dummy], axis =1 )
iris = iris.iloc[0:100, :] # 截取前一百行樣本

# 構建Logistic Regression , 對Species是否為setosa進行分類 setosa ~ Sepal.Length
# Y = g(BX) = 1/(1+exp(-BX))
def logit(x):
    return 1./(1+np.exp(-x))

temp = pd.DataFrame(iris.iloc[:, 0])
temp['x0'] = 1.
X = temp.iloc[:,[1,0]]
Y = iris['setosa'].reshape(len(iris), 1) #整理出X矩陣 和 Y矩陣

# 批量梯度下降法
m,n = X.shape #矩陣大小
alpha = 0.0065 #設定學習速率
theta_g = np.zeros((n,1)) #初始化參數
maxCycles = 3000 #迭代次數
J = pd.Series(np.arange(maxCycles, dtype = float)) #損失函數

for i in range(maxCycles):
    h = logit(dot(X, theta_g)) #估計值  
    J[i] = -(1/100.)*np.sum(Y*np.log(h)+(1-Y)*np.log(1-h)) #計算損失函數值      
    error = h - Y #誤差
    grad = dot(X.T, error) #梯度
    theta_g -= alpha * grad
print theta_g
print J.plot()   

# 牛頓方法
theta_n = np.zeros((n,1)) #初始化參數
maxCycles = 10 #迭代次數
C = pd.Series(np.arange(maxCycles, dtype = float)) #損失函數
for i in range(maxCycles):
    h = logit(dot(X, theta_n)) #估計值 
    C[i] = -(1/100.)*np.sum(Y*np.log(h)+(1-Y)*np.log(1-h)) #計算損失函數值      
    error = h - Y #誤差
    grad = dot(X.T, error) #梯度
    A =  h*(1-h)* np.eye(len(X)) 
    H = np.mat(X.T)* A * np.mat(X) #海瑟矩陣, H = X`AX
    theta_n -= inv(H)*grad
print theta_n
print C.plot()   

代碼所用的數據下載地址：http://files.cnblogs.com/files/sumai/iris.rar