ML(4)——逻辑回归

  Logistic Regression虽然名字里带“回归”,但是它实际上是一种分类方法,“逻辑”是Logistic的音译,和真正的逻辑没有任何关系。

模型

线性模型

  由于逻辑回归是一种分类方法,所以我们仍然以最简的二分类为例。与感知机不同,对于逻辑回归的分类结果,y ∈ {0, 1},我们需要找到最佳的hθ(x)拟合数据。

  这里容易联想到线性回归。线性回归也可以用于分类,但是很多时候,尤其是二分类的时候,线性回归并不能很好地工作,因为分类不是连续的函数,其结果只能是固定的离散值。设想一下有线性回归得到的拟合曲线hθ(x),当x→∞时,有可能y→∞,这就无法对y ∈ {0, 1}进行有效解释。

  对于二分类,逻辑回归的目的是找到一个函数,使得无论x取何值,都有:

  满足这个式子的典型函数是sigmoid函数,也称为logistic函数:

 

  在sigmoid函数g(z)中:

  现在,将hθ(x)赋予sigmoid函数g(z)的特性:

 

  其中:

  最终,逻辑回归的模型函数:

 

  假设给定一些输入,现在需要根据逻辑回归模型预测肿瘤是否是良性,最终得到hθ(x) = 0.8,可以用概率表述:

 

  上式表示在当前输入下,y=1的概率是0.8,y=0的概率是0.2,因为是分类,所以判断y = 1。

  需要注意的是,sigmoid函数不是样本点的分隔曲线,它表示的是逻辑回归的测结果;θTx才是分隔曲线,它将样本点分为θTx ≥ 0和θTx < 0两部分:

分隔曲线

 

Sigmoid函数,最终模型

  由此看来,逻辑回归的线性模型同样是找到最佳的θ,使两类样本点分离,这就在很大程度上和感知机相似。

多项式模型

  直观地看,线性模型的决策边界就是将两类样本点分离开的分隔曲线,我们之前已经多次接触过,只是没有给它起一个专业的名字。假设在一个模型中,hθ(x) = g(θ0 + θ1x1 + θ2x2) = g(-3 + x1 + x2),那么决策边界就是 -3 + x1 + x2 = 0:

  很多时候,直线并不能很好地作为决策边界,如下图所示:

  此时需要使用多项式模型添加更多的特征:

  这相当于添加了两个新的特征:

  

  深入了解不同函数的特征有助于选择正确的模型。添加的特征越多,曲线越复杂,对训练样本的拟合度越高,同时也更容易导致过拟合而丧失泛化性:

 

  关于过拟合及其处理,可参考《ML(附录3)——数据拟合和正则化》。

多分类

  我们已经知道了二分类的模型,然而实际问题中往往不只是二分类,那么逻辑回归如何处理多分类呢?

  一种可行的方法就是化繁为简,将多分类转换为二分类:

  如上图所示,三分类转换为三个二分类,其中上标表示分类的类别:

 

  对于输入x,其预测结果是所有hθ(x)中值最大的一个。对于最后的预测结论,以上面的三分类为例,如果输入一个标签为2的特征集,对于hθ(0)(x) 来说,hθ(0)(x) < 0.5:

 

  对于hθ(1)(x) 来说,hθ(1)(x) < 0.5:

 

  对于hθ(3)(x) 来说,hθ(3)(x) ≥ 0.5:

  因此,对于输入x,其预测结果是所有hθ(x)中值最大的一个。至于每个样本标签值是多少,无所谓了,在训练每个hθ(i)(x)前,都需要把y(i) 转换为1,其余转换为0。

  在上面的图组中也可以看出,对于有k个标签的多分类,需要训练k个不同的逻辑回归模型。

学习策略

  在感知机中,模型函数使用了sign,由于sign是阶跃函数,不易优化,所以为了求得损失函数,我们使用函数间隔进行了一系列转换。对于逻辑回归,由于函数本身是连续曲线,所以不会存在这样的问题,其J(θ)使用对数损失函数。对数损失函数的原型:

 

  将其用在逻辑回归上:

  这里的对数是以e为底的对数,即:

  注意,0 < hθ(x) < 1。上图是y = 1时的costfunction,可以看到,当hθ(x)→1时,Cost(hθ(x), y)→0;当hθ(x)→0时,Cost(hθ(x),y)→∞。也就是说当分类是1时,sigmoid的值越接近于1,损失值越小;sigmoid的值越接近于0,损失值越大。损失值越大,分类点越接近决策面,其分类越模糊。与此类似,下图是 y = 0时的cost function:

  Cost function可以把y = 1和y=1两种情况合并到一起:

 

  理解这种形式需要再次将y分开:

  最后得到J(θ)的最终形式:

 

  注意,这里hθ(x)是sigmoid函数:

 

  如果写成矩阵的形式:

  上式中,hθ(x)的操作对应矩阵中的每个元素,1-Y,log hθ(x)也一样,可参照后文的代码实现来理解。

算法

  与之前的算法一样,我们的目的是找到最佳的θ使得J(θ)最小化,将求解θ转换为最优化问题,即:

梯度下降

  梯度下降是一个适用范围很广的方法,这里同样可以使用梯度下降求解θ:

 

  关于梯度下降的更多内容可参考《ML(附录1)——梯度下降》。

  在对J(θ)求偏导时先做一些准备工作,计算sigmoid函数的导数(关于偏导和一元函数的导数,可参考《多变量微积分》和《单变量微积分》的相关章节):

 

  现在计算m=1时J(θ)的偏导,此时可以删除上标:

  推广到m个样本:

  在《机器学习实战》中提到对极大似然数使用梯度上升求最大值,最后得到:

 

  这和对损失函数采用梯度下降求最小值是一样的,因为损失函数使用了似然数的负数形式,Cost(X, Y) = -logP(Y|X),所以对-logP(Y|X)梯度下降和对+logP(Y|X)梯度上升将得到同样的结果。

  对于多项式模型,需要预先添加特征,使得每个θj都有唯一的xj对应:

 

  如果用矩阵表示:

  

  在此基础上使用L2正则化(关于正则化,可参考《ML(附录3)——过拟合与欠拟合》):

 

代码实现

ex2data1.txt:

  1 34.62365962451697,78.0246928153624,0
  2 30.28671076822607,43.89499752400101,0
  3 35.84740876993872,72.90219802708364,0
  4 60.18259938620976,86.30855209546826,1
  5 79.0327360507101,75.3443764369103,1
  6 45.08327747668339,56.3163717815305,0
  7 61.10666453684766,96.51142588489624,1
  8 75.02474556738889,46.55401354116538,1
  9 76.09878670226257,87.42056971926803,1
 10 84.43281996120035,43.53339331072109,1
 11 95.86155507093572,38.22527805795094,0
 12 75.01365838958247,30.60326323428011,0
 13 82.30705337399482,76.48196330235604,1
 14 69.36458875970939,97.71869196188608,1
 15 39.53833914367223,76.03681085115882,0
 16 53.9710521485623,89.20735013750205,1
 17 69.07014406283025,52.74046973016765,1
 18 67.94685547711617,46.67857410673128,0
 19 70.66150955499435,92.92713789364831,1
 20 76.97878372747498,47.57596364975532,1
 21 67.37202754570876,42.83843832029179,0
 22 89.67677575072079,65.79936592745237,1
 23 50.534788289883,48.85581152764205,0
 24 34.21206097786789,44.20952859866288,0
 25 77.9240914545704,68.9723599933059,1
 26 62.27101367004632,69.95445795447587,1
 27 80.1901807509566,44.82162893218353,1
 28 93.114388797442,38.80067033713209,0
 29 61.83020602312595,50.25610789244621,0
 30 38.78580379679423,64.99568095539578,0
 31 61.379289447425,72.80788731317097,1
 32 85.40451939411645,57.05198397627122,1
 33 52.10797973193984,63.12762376881715,0
 34 52.04540476831827,69.43286012045222,1
 35 40.23689373545111,71.16774802184875,0
 36 54.63510555424817,52.21388588061123,0
 37 33.91550010906887,98.86943574220611,0
 38 64.17698887494485,80.90806058670817,1
 39 74.78925295941542,41.57341522824434,0
 40 34.1836400264419,75.2377203360134,0
 41 83.90239366249155,56.30804621605327,1
 42 51.54772026906181,46.85629026349976,0
 43 94.44336776917852,65.56892160559052,1
 44 82.36875375713919,40.61825515970618,0
 45 51.04775177128865,45.82270145776001,0
 46 62.22267576120188,52.06099194836679,0
 47 77.19303492601364,70.45820000180959,1
 48 97.77159928000232,86.7278223300282,1
 49 62.07306379667647,96.76882412413983,1
 50 91.56497449807442,88.69629254546599,1
 51 79.94481794066932,74.16311935043758,1
 52 99.2725269292572,60.99903099844988,1
 53 90.54671411399852,43.39060180650027,1
 54 34.52451385320009,60.39634245837173,0
 55 50.2864961189907,49.80453881323059,0
 56 49.58667721632031,59.80895099453265,0
 57 97.64563396007767,68.86157272420604,1
 58 32.57720016809309,95.59854761387875,0
 59 74.24869136721598,69.82457122657193,1
 60 71.79646205863379,78.45356224515052,1
 61 75.3956114656803,85.75993667331619,1
 62 35.28611281526193,47.02051394723416,0
 63 56.25381749711624,39.26147251058019,0
 64 30.05882244669796,49.59297386723685,0
 65 44.66826172480893,66.45008614558913,0
 66 66.56089447242954,41.09209807936973,0
 67 40.45755098375164,97.53518548909936,1
 68 49.07256321908844,51.88321182073966,0
 69 80.27957401466998,92.11606081344084,1
 70 66.74671856944039,60.99139402740988,1
 71 32.72283304060323,43.30717306430063,0
 72 64.0393204150601,78.03168802018232,1
 73 72.34649422579923,96.22759296761404,1
 74 60.45788573918959,73.09499809758037,1
 75 58.84095621726802,75.85844831279042,1
 76 99.82785779692128,72.36925193383885,1
 77 47.26426910848174,88.47586499559782,1
 78 50.45815980285988,75.80985952982456,1
 79 60.45555629271532,42.50840943572217,0
 80 82.22666157785568,42.71987853716458,0
 81 88.9138964166533,69.80378889835472,1
 82 94.83450672430196,45.69430680250754,1
 83 67.31925746917527,66.58935317747915,1
 84 57.23870631569862,59.51428198012956,1
 85 80.36675600171273,90.96014789746954,1
 86 68.46852178591112,85.59430710452014,1
 87 42.0754545384731,78.84478600148043,0
 88 75.47770200533905,90.42453899753964,1
 89 78.63542434898018,96.64742716885644,1
 90 52.34800398794107,60.76950525602592,0
 91 94.09433112516793,77.15910509073893,1
 92 90.44855097096364,87.50879176484702,1
 93 55.48216114069585,35.57070347228866,0
 94 74.49269241843041,84.84513684930135,1
 95 89.84580670720979,45.35828361091658,1
 96 83.48916274498238,48.38028579728175,1
 97 42.2617008099817,87.10385094025457,1
 98 99.31500880510394,68.77540947206617,1
 99 55.34001756003703,64.9319380069486,1
100 74.77589300092767,89.52981289513276,1
View Code

 

  

Octave

%% Machine Learning Online Class - Exercise 2: Logistic Regression
%
%  Instructions
%  ------------
% 
%  This file contains code that helps you get started on the logistic
%  regression exercise. You will need to complete the following functions 
%  in this exericse:
%
%     sigmoid.m
%     costFunction.m
%     predict.m
%     costFunctionReg.m
%
%  For this exercise, you will not need to change any code in this file,
%  or any other files other than those mentioned above.
%

%% Initialization
clear ; close all; clc

%% Load Data
%  The first two columns contains the exam scores and the third column
%  contains the label.

data = load('ex2data1.txt');
X = data(:, [1, 2]); y = data(:, 3);

%% ==================== Part 1: Plotting ====================
%  We start the exercise by first plotting the data to understand the 
%  the problem we are working with.

fprintf(['Plotting data with + indicating (y = 1) examples and o ' ...
         'indicating (y = 0) examples.\n']);

plotData(X, y);

% Put some labels 
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;

fprintf('\nProgram paused. Press enter to continue.\n');
pause;


%% ============ Part 2: Compute Cost and Gradient ============
%  In this part of the exercise, you will implement the cost and gradient
%  for logistic regression. You neeed to complete the code in 
%  costFunction.m

%  Setup the data matrix appropriately, and add ones for the intercept term
[m, n] = size(X);

% Add intercept term to x and X_test
X = [ones(m, 1) X];

% Initialize fitting parameters
initial_theta = zeros(n + 1, 1);

% Compute and display initial cost and gradient
[cost, grad] = costFunction(initial_theta, X, y);

fprintf('Cost at initial theta (zeros): %f\n', cost);
fprintf('Expected cost (approx): 0.693\n');
fprintf('Gradient at initial theta (zeros): \n');
fprintf(' %f \n', grad);
fprintf('Expected gradients (approx):\n -0.1000\n -12.0092\n -11.2628\n');

% Compute and display cost and gradient with non-zero theta
test_theta = [-24; 0.2; 0.2];
[cost, grad] = costFunction(test_theta, X, y);

fprintf('\nCost at test theta: %f\n', cost);
fprintf('Expected cost (approx): 0.218\n');
fprintf('Gradient at test theta: \n');
fprintf(' %f \n', grad);
fprintf('Expected gradients (approx):\n 0.043\n 2.566\n 2.647\n');

fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============= Part 3: Optimizing using fminunc  =============
%  In this exercise, you will use a built-in function (fminunc) to find the
%  optimal parameters theta.

%  Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 400);
%  Run fminunc to obtain the optimal theta
%  This function will return theta and the cost 
[theta, cost] = fminunc(@(t)(costFunction(t, X, y)), initial_theta, options);
% Print theta to screen
fprintf('Cost at theta found by fminunc: %f\n', cost);
fprintf('Expected cost (approx): 0.203\n');
fprintf('theta: \n');
fprintf(' %f \n', theta);
fprintf('Expected theta (approx):\n');
fprintf(' -25.161\n 0.206\n 0.201\n');

% Plot Boundary
plotDecisionBoundary(theta, X, y);

% Put some labels 
hold on;
% Labels and Legend
xlabel('Exam 1 score')
ylabel('Exam 2 score')

% Specified in plot order
legend('Admitted', 'Not admitted')
hold off;


fprintf('\nProgram paused. Press enter to continue.\n');
pause;

%% ============== Part 4: Predict and Accuracies ==============
%  After learning the parameters, you'll like to use it to predict the outcomes
%  on unseen data. In this part, you will use the logistic regression model
%  to predict the probability that a student with score 45 on exam 1 and 
%  score 85 on exam 2 will be admitted.
%
%  Furthermore, you will compute the training and test set accuracies of 
%  our model.
%
%  Your task is to complete the code in predict.m

%  Predict probability for a student with score 45 on exam 1 
%  and score 85 on exam 2 

prob = sigmoid([1 45 85] * theta);
fprintf(['For a student with scores 45 and 85, we predict an admission ' ...
         'probability of %f\n'], prob);
fprintf('Expected value: 0.775 +/- 0.002\n\n');

% Compute accuracy on our training set
p = predict(theta, X);

fprintf('Train Accuracy: %f\n', mean(double(p == y)) * 100);
fprintf('Expected accuracy (approx): 89.0\n');
fprintf('\n');

  plotData.m

function plotData(X, y)
%PLOTDATA Plots the data points X and y into a new figure 
%   PLOTDATA(x,y) plots the data points with + for the positive examples
%   and o for the negative examples. X is assumed to be a Mx2 matrix.

% Create New Figure
  figure; hold on;

% Instructions: Plot the positive and negative examples on a
%               2D plot, using the option 'k+' for the positive
%               examples and 'ko' for the negative examples.
  pos = find(y==1); neg = find(y == 0);

  plot(X(pos, 1), X(pos, 2), 'k+','LineWidth', 2, 'MarkerSize', 7);
  plot(X(neg, 1), X(neg, 2), 'ko', 'MarkerFaceColor', 'y', 'MarkerSize', 7);

  hold off;

end

  sigmoid.m

function g = sigmoid(z)
%SIGMOID Compute sigmoid function
%   g = SIGMOID(z) computes the sigmoid of z.

% You need to return the following variables correctly 
  g = ones(size(z)) ./ (1 + exp(-1 * z));

end

  costFunction.m

function [J, grad] = costFunction(theta, X, y)
  %COSTFUNCTION Compute cost and gradient for logistic regression
  %   J = COSTFUNCTION(theta, X, y) computes the cost of using theta as the
  %   parameter for logistic regression and the gradient of the cost
  %   w.r.t. to the parameters.

  % Initialize some useful values
  m = length(y); % number of training examples

  % You need to return the following variables correctly 
  J = 0;
  grad = zeros(size(theta));

  % Instructions: Compute the cost of a particular choice of theta.
  %               You should set J to the cost.
  %               Compute the partial derivatives and set grad to the partial
  %               derivatives of the cost w.r.t. each parameter in theta
  %
  % Note: grad should have the same dimensions as theta
  %

  %  % use intrator to compue J
  %  for i = 1:m
  %    theta_X = X(i,:) * theta;
  %    h = 1 / (1 + exp(-1 * theta_X));
  %    J += y(i) * log(h) + (1 - y(i)) * log(1 - h);
  %  end
  %  J /= -1 * m;
  
  % use matrix to compute gradient
  h = sigmoid(X * theta);
  J = (y' * log(h) + (1 - y)' * log(1 - h)) / (-1 * m);
  #J /= -1 * m;
  
  grad = X' * (h - y) / m;
end

Python

  1 from __future__ import division
  2 import numpy as np
  3 import random
  4 import matplotlib.pyplot as plt
  5 
  6 def train(X, Y, iterateNum=10000000, alpha=0.003):
  7     '''
  8     :param X: 训练样本的特征集
  9     :param Y: 训练样本的标签
 10     :param iterateNum: 梯度下降的迭代次数
 11     :param alpha: 学习率
 12     :return:theta
 13     '''
 14     m, n = np.shape(X)
 15     theta = np.zeros((n + 1, 1))
 16     # 在第一列添加x0
 17     X_new = np.c_[np.ones(m), X]
 18 
 19     for i in range(iterateNum):
 20         m = np.shape(X_new)[0]
 21         h = h_function(X_new, theta)
 22         theta -= alpha * (np.dot(X_new.T, h - Y) / m)
 23 
 24         if i % 100000 == 0:
 25             print('\t---------iter=' + str(i) + ', J(θ)=' + str(J_function(X_new, Y, theta)))
 26 
 27     print( str(J_function(X_new, Y, theta)))
 28     return theta
 29 
 30 def h_function(X, theta):
 31     return sigmoid(np.dot(X, theta))
 32 
 33 def sigmoid(X):
 34     return 1 / (1 + np.exp(-X ))
 35 
 36 # 计算J(θ)
 37 def J_function(X, Y, theta):
 38     h = h_function(X, theta)
 39     J_1 = np.dot(Y.T, np.log(h))
 40     J_2 = np.dot(1 - Y.T, np.log(1 - h))
 41     m = np.shape(X)[0]
 42     J = (-1 / m) * (J_1 + J_2)
 43 
 44     return J
 45 
 46 def predict(x, theta):
 47     if h_function(x, theta) >= 0.5:
 48         return 1
 49     else:
 50         return 0
 51 
 52 # 归一化处理
 53 def normalization(X):
 54     m, n = np.shape(X)
 55     X_new = np.zeros((m, n))
 56 
 57     for j in range(n):
 58         max = np.max(X[:,j])
 59         min = np.min(X[:,j])
 60         d_value = max - min
 61         for i in range(m):
 62             X_new[i, j] = (X[i, j] - min) / d_value
 63 
 64     return X_new
 65 
 66 def plot_datas(X, Y, theta):
 67     plt.figure()
 68 
 69     # 绘制分隔直线 g = 0
 70     x1 = [0, 1]
 71     x2 = [(-1 / theta[2]) * (theta[0] + theta[1] * x1[0]),
 72           (-1 / theta[2]) * (theta[0] + theta[1] * x1[1])]
 73     plt.xlabel('x1')
 74     plt.ylabel('x2')
 75 
 76     plt.plot(x1, x2, color='b')
 77 
 78     # 绘制数据点
 79     admit_x1, admit_x2 = [],[]
 80     not_admit_x1, not_admit_x2 = [],[]
 81     for i in range(len(X)):
 82         if (Y[i] == 1):
 83             admit_x1.append(X[i][0])
 84             admit_x2.append(X[i][1])
 85         else:
 86             not_admit_x1.append(X[i][0])
 87             not_admit_x2.append(X[i][1])
 88 
 89     plt.scatter(admit_x1, admit_x2, color='g')
 90     plt.scatter(not_admit_x1, not_admit_x2, marker='x', color='r')
 91 
 92     plt.legend(['logistic line', 'Admitted', 'Not admitted'])
 93     plt.show()
 94 
 95 if __name__ == '__main__':
 96     train_datas = np.loadtxt('ex2data1.txt', delimiter=',')
 97     X = train_datas[:,[0, 1]]
 98     X = normalization(X)
 99     Y = train_datas[:,[2]]
100     theta = train(X, Y)
101 
102     print(theta)
103     plot_datas(X, Y, theta)

 

 

  对于样本数据,上面的梯度下降并不是非常有效,无论是否预处理数据(归一化或其他方法),都必须反复调整学习率和迭代次数。如果学习率过大,算法将不会收敛;如果过小,算法收敛的十分缓慢,需要增加迭代次数。代码中的参数最终将使算法收敛于0.203。

Sklearn

 1 from sklearn.linear_model import LogisticRegression
 2 import numpy as np
 3 
 4 if __name__ == '__main__':
 5     train_datas = np.loadtxt("ex2data1.txt", delimiter=',')
 6     X_train = train_datas[:,[0, 1]]
 7     Y_train = train_datas[:,[2]]
 8 
 9     logistic = LogisticRegression()
10     logistic.fit(X_train, Y_train)
11 
12     theta = [logistic.intercept_[0], logistic.coef_[0]]
13     print(theta)

 

  参考:

  Ng视频《Logistic Regression》

  周志华《机器学习》

  《机器学习导论》

  Peter Flach《机器学习》

 


   作者:我是8位的

  出处:http://www.cnblogs.com/bigmonkey

  本文以学习、研究和分享为主,如需转载,请联系本人,标明作者和出处,非商业用途! 

  扫描二维码关注公众号“我是8位的”

  

posted on 2018-06-29 17:31  我是8位的  阅读(1142)  评论(0编辑  收藏  举报

导航