一：matlab實現
1.資料的Excel處理

西瓜資料集3.0

2.程式碼

# -*- coding: utf-8 -*-

old_l = 0;
n = 0;
b = [0;0;1]; %對應書中（3.25）下的B=(w;b)，因為x有兩個屬性：密度，含糖率，所以有b三行，還有一個是w*x+b中的b。

x = xlsread('E:\Program Files\octave\西瓜3.0.xlsx', 'Sheet1', 'A1:Q3')
y = xlsread('E:\Program Files\octave\西瓜3.0.xlsx', 'Sheet1', 'A4:Q4')

while(1
 
)
  cur_l = 0;
  bx = zeros(17,1);

  for i=1:17
    bx(i) = b.'*x(:,i);
    cur_l = cur_l + (-y(i)*bx(i)) + log(1+exp(bx(i))); %對應式（3.27）
  end

  if abs(cur_l-old_l) < 0.0001
    break;
  end

  n += 1;
  old_l = cur_l;
  dl = 0;
  d2l = 0;

  for i=1:17 %牛頓法求解（3.29）
    pl(i) = 1 - 1 / (1+exp(bx(i)));
    dl = dl - x(:,i) * (y(i) - pl(i)); %（3.30
 
）
    d2l += x(:,i) * x(:,i).'* pl(i)*(1-pl(i)); %（3.31）
  end
  b -= d2l \ dl;
end

%繪圖
for i=1:17
  if y(i) == 1
    plot(x(1,i),x(2,i),'+r');
    hold on;
  else if y(i) == 0
    plot(x(1,i),x(2,i),'og');
    hold on;
  end
end
end  


ply = -(0.1*b(1)+b(3))/b(2);
pry = -(0.9*b(1)+b(3))/b(2);     
line([0.1
 
 0.9],[ply pry]);        

xlabel('density');
ylabel('Sugar content');
title('Rate regression');

3.繪圖結果

二：Python實現
1.Excel處理資料

2.程式碼

# -*- coding: utf-8 -*-
#對率迴歸分類
import numpy as np
from numpy import linalg
import pandas as pd
#讀取資料集
inputfile = 'E:/Program Files/octave/xigua.xlsx'
data_original = pd.read_excel(inputfile, 'Sheet1')
#資料的初步轉化與操作--屬性x變數2行17列陣列，並新增一組1作為吸入的偏置x^=（x;1）
x=np.array([list(data_original[u'密度']),list(data_original[u'含糖率']),[1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1]])
y=np.array([1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0])
#定義初始引數
beta = np.array([[0],[0],[1]])       #β列向量
old_l = 0 #3.27式l值的記錄，這是上一次迭代的l值
n=0

while 1:
    beta_T_x = np.dot(beta.T[0], x)  # 對β進行轉置取第一行（因為β轉置後是array([[0, 0, 1]]，取第一行得到array([0, 0, 1])
                                      # ，再與x相乘（dot）,beta_T_x表示β轉置乘以x)
    cur_l = 0   #當前的l值
    for i in range(17):
        cur_l = cur_l + ( -y[i]*beta_T_x[i]+np.log(1+np.exp(beta_T_x[i])) )#計算當前3.27式的l值，這是目標函式，希望他越小越好
    #迭代終止條件
    if np.abs(cur_l - old_l) <= 0.000001:   #精度，二者差在0.000001以內就認為可以了，說明l已經很收斂了
        break               #滿足條件直接跳出迴圈

    #牛頓迭代法更新β
    #求關於β的一階導數和二階導數
    n=n+1
    old_l = cur_l
    dbeta = 0
    d2beta = 0
    for i in range(17):
        dbeta = dbeta - np.dot(np.array([x[:,i]]).T,( y[i]-(  np.exp(beta_T_x[i])/(1+np.exp(beta_T_x[i])) ) )) #一階導數
        d2beta =d2beta + np.dot(np.array([x[:,i]]).T,np.array([x[:,i]]).T.T) * (  np.exp(beta_T_x[i])/(1+np.exp(beta_T_x[i])) ) * (1-(  np.exp(beta_T_x[i])/(1+np.exp(beta_T_x[i])) ))
    beta = beta - np.dot(linalg.inv(d2beta),dbeta)
print('模型引數是：',beta)
print('迭代次數：',n)

3.執行結果