這個例子是《機器學習實戰》（<machine learning in action>）邏輯迴歸的一個例項：從疝氣病症預測病馬的死亡率。

疝病是描述馬胃腸痛的術語。該資料集中包含了醫院檢查馬疝病的一些指標，我們的目標是通過這些指標（特徵），來預測馬是否會死亡。

資料集包括229個訓練樣本和67各測試樣本，特徵數量為22。資料集中包含缺失值，採取的措施是用0替換所有缺失值，這樣做的原因是，我們在更新theta值時使用下式（具體參考上一篇），當某個特徵值缺失時，我們如果用0來替代，則對應的theta值不會被更新（公式右邊第二項=0），即缺失值不會對引數造成影響。

以下為python程式碼，由於訓練資料比較少，這邊使用了批處理梯度下降法，沒有使用增量梯度下降法。

##author:lijiayan
##data:2016/10/27
##name:logReg.py
from numpy import *
import matplotlib.pyplot as plt


def loadData(filename):
    data = loadtxt(filename)
    m,n = data.shape
    print 'the number of  examples:',m
    print 'the number of features:',n-1
x = data[:,0:n-1]
    y = data[:,n-1:n]
    return 
 
x,y


#the sigmoid function
def sigmoid(z):
    return 1.0 / (1 + exp(-z))


#the cost function
def costfunction(y,h):
    y = array(y)
    h = array(h)
    J = sum(y*log(h))+sum((1-y)*log(1-h))
    return J


# the batch gradient descent algrithm
def gradescent(x,y):
    m,n = shape(x)     #m: number of training example; n: number of features
 
x = c_[ones(m),x]     #add x0
x = mat(x)      # to matrix
y = mat(y)
    a = 0.0000025       # learning rate
maxcycle = 4000
theta = zeros((n+1,1))  #initial theta
J = []
    for i in range(maxcycle):
        h = sigmoid(x*theta)
        theta = theta + a * (x.T)*(y-h)
        cost = costfunction(y,h)
        J.append(cost)

    plt.plot(J)
    plt.show()
    return theta,cost


#the stochastic gradient descent (m should be large,if you want the result is good)
def stocGraddescent(x,y):
    m,n = shape(x)     #m: number of training example; n: number of features
x = c_[ones(m),x]     #add x0
x = mat(x)      # to matrix
y = mat(y)
    a = 0.01       # learning rate
theta = ones((n+1,1))    #initial theta
J = []
    for i in range(m):
        h = sigmoid(x[i]*theta)
        theta = theta + a * x[i].transpose()*(y[i]-h)
        cost = costfunction(y,h)
        J.append(cost)
    plt.plot(J)
    plt.show()
    return theta,cost


#plot the decision boundary
def plotbestfit(x,y,theta):
    plt.plot(x[:,0:1][where(y==1)],x[:,1:2][where(y==1)],'ro')
    plt.plot(x[:,0:1][where(y!=1)],x[:,1:2][where(y!=1)],'bx')
    x1= arange(-4,4,0.1)
    x2 =(-float(theta[0])-float(theta[1])*x1) /float(theta[2])

    plt.plot(x1,x2)
    plt.xlabel('x1')
    plt.ylabel(('x2'))
    plt.show()


def classifyVector(inX,theta):
    prob = sigmoid((inX*theta).sum(1))
    return where(prob >= 0.5, 1, 0)


def accuracy(x, y, theta):
    m = shape(y)[0]
    x = c_[ones(m),x]
    y_p = classifyVector(x,theta)
    accuracy = sum(y_p==y)/float(m)
    return accuracy

呼叫上面程式碼：

from logReg import *
x,y = loadData("horseColicTraining.txt")
theta,cost = gradescent(x,y)
print 'J:',cost

ac_train = accuracy(x, y, theta)
print 'accuracy of the training examples:', ac_train

x_test,y_test = loadData('horseColicTest.txt')
ac_test = accuracy(x_test, y_test, theta)
print 'accuracy of the test examples:', ac_test

學習速率=0.0000025，迭代次數=4000時的結果：

似然函式走勢（J = sum(y*log(h))+sum((1-y)*log(1-h))），似然函式是求最大值，一般是要穩定了才算最好。

下圖為計算結果，可以看到訓練集的準確率為73%，測試集的準確率為78%。

這個時候，我去看了一下資料集，發現沒個特徵的數量級不一致，於是我想到要進行歸一化處理：

歸一化處理句修改列loadData(filename)函式：

def loadData(filename):
    data = loadtxt(filename)
    m,n = data.shape
    print 'the number of  examples:',m
    print 'the number of features:',n-1
x = data[:,0:n-1]
    max = x.max(0)
    min = x.min(0)
    x = (x - min)/((max-min)*1.0)     #scaling
y = data[:,n-1:n]
    return x,y

在沒有歸一化的時候，我的學習速率取了0.0000025（加大就會震盪，因為有些特徵的值很大，學習速率取的稍大，波動就很大），由於學習速率小，迭代了4000次也沒有完全穩定。現在當把特徵歸一化後（所有特徵的值都在0~1之間），這樣學習速率可以加大，迭代次數就可以大大減少，以下是學習速率=0.005，迭代次數=500的結果：

此時的訓練集的準確率為72%，測試集的準確率為73%

從上面這個例子，我們可以看到對特徵進行歸一化操作的重要性。