如何优化逻辑回归(logistic regression)？ - 知乎352被浏览21618分享邀请回答653 条评论分享收藏感谢收起2添加评论分享收藏感谢收起查看更多回答机器学习笔记(四)之Logistic回归 - gzj_1101的专栏 - CSDN博客
机器学习笔记(四)之Logistic回归
machine learning
回归的概念
假设有一些数据点，我们利用一条直线对这些点进行拟合(该直线为最佳拟合直线)，这个拟合过程称为回归。
logistic回归思想
根据根据现有数据集对分类边界线建立回归公式，以此进行分类。
logisitc回归一般过程
收集数据：采用任意方法收集数据
准备数据：由于需要计算距离，所以要求数据类型是数值型。另外，结构化数据格式最佳
分析数据：采用任意方法对数据进行分析
训练算法：大部分时间用于训练，训练的目的用于找到最佳分类回归系数
测试算法：一旦训练步骤完成，分类将会很快
使用算法：首先，我们需要输入一些数据，然后转化成对应的结构化数值；接着，用训练好的数据对这些数值进行简单的回归计算判断他们输入那个类别；
逻辑斯蒂回归模型
逻辑斯蒂分布
设X是连续随机变量，X服从逻辑斯谛分布是指X具有下面的分布函数和密度函数：
式中，u为位置参数，y&0为形状参数。这里可以类比
逻辑斯谛分布的密度函数f(x)和分布函数F(X)的图形如下：分布函数属于logistic函数，其图形属于一条S型曲线(sigmoid curve,神经网络里常用的一种激活函数)。该曲线以点(u,1/2)为中心对称。
曲线在中心附近增长速度较快，在两端的增长速度较慢。形状参数y的值越小，曲线在中心附近增长的越快。
二项逻辑斯蒂回归模型
二项逻辑斯蒂回归模型(binary logistic regression model)是一种分类模型，由条件概率P(Y|X)表示，形式为条件化的逻辑斯蒂分。随机变量X取值为实数，随机变量Y取值为1或者0
逻辑斯第回归模型
w称为权值参数，b称为偏置，w*x是w和x的內积。
对于给定的输入实例x，按照式子(6.3)和(6.4)可以求得P(Y=1|x)和P(Y=0|x)。逻辑斯蒂回归比较两个条件概率的大小，将实例x划分到概率大的分类。
即我们将每一个特征乘以一个回归系数，将所有结果值相加，得到的总和带入到sigmoid函数中，进而得到一个0-1之间的数据，如果数值大于0.5，则被划分到1类中，否则将被划分到0类中。所以logistic回归也可以看做是一种概率估计。
逻辑斯蒂回归模型的特点
一个时间的几率(odds)是指该事件发生的概率与不发生的概率的比值。如果该事件发生的概率为p，那么该事件发生的几率为p/(1-p),该事件的对数几率或者logit函数是
对于逻辑斯蒂回归而言，由式(6.5)与式(6.6)得
这就是说，在逻辑斯蒂回归模型当中，输出Y=1的对对书几率时输入x的线性函数或者说，输出Y=1的对数几率是由输入x的线性函数表示的模型。
模型参数估计
逻辑斯蒂回归模型学习时，对于给定的训练集T={(x1,y1),(x2,y3),(x3,y3)…..(xn,yn)}，其中xi属于R，yi属于{0,1}，可以应用计参数模型，从而得到逻辑斯蒂回归模型：
设似然函数为
对数似然函数为:
对L(w)求极大值，得到w的估计值
Logistic回归的优缺点
优点：计算代价不高易于理解和实现。
缺点：容易欠拟合，分类精度不高。
适用数据类型：数值型和标称型
最佳回归系数的确定
梯度上升算法
梯度上升法的基本思想：要找到函数的最大值，最好的方法是沿着该函数的梯度方向探寻，如果梯度标记为，则函数f(x,y)的梯度由下面式子表示
这个梯度意味着沿着x方向的梯度为?f(x,y)?x
沿着y方向的梯度为?f(x,y)?y
函数f(x,y)在(x,y)处可微。
梯度算子总是指向函数值增长最快的方向。我们规定一个步长a表示的是移动量的大小。所以梯度上升算法的迭代公式为w=w+α?f(w)
该公式将一直被迭代执行，直到到达某个停止条件为止比如迭代次数下降到某个指定值或者误差在允许范围内。
梯度下降算法就是沿着反方向执行迭代，实际原理是一样的。
梯度上升算法的伪代码
每个回归系数初始化为1
计算整个数据集的梯度
使用alpha*gradient跟新回归系数的向量
返回回归系数
Logistic回归梯度上升优化算法
def loadDataSet():
dataMat = []
labelMat = []
fr = open('testSet.txt')
for line in fr.readlines():
lineArr = line.strip().split()
dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
labelMat.append(int(lineArr[2]))
return dataMat,labelMat
首先loadDataset函数从test.txt里面讲数据读取出来，testSet.txt每行包括三个数据，每行的前两个值分别为X1和X2，第三个值表示的对应的标签类别。
下面是梯度上升的代码
def sigmoid(intX):
return 1.0/(1+np.exp(-intX))
def gradAscent(dataMatIn,classLabels):
dataMatrix = np.mat(dataMatIn)
labelMat = np.mat(classLabels).transpose()
m,n = np.shape(dataMatrix)
alpha = 0.001
maxCycles = 500
weights = np.ones((n,1))
for k in range(maxCycles):
t = dataMatrix * weights
h = sigmoid(t)
error = labelMat-h
weights = weights + alpha *dataMatrix.transpose()*error
return weights
最后weights = weights + alpha *dataMatrix.transpose()*error的理解
w=w+αf(x)△x
其中α对应alpha，error对应△x
分析数据，画出决策树
def plotBestFit(weights):
dataMat,labelMat = loadDataSet()
dataArr = np.array(dataMat)
n = np.shape(dataArr)[0]
xcord1 = [];ycord1 = []
xcord2 = [];ycord2 = []
for i in range(n):
if int(labelMat[i]) == 1:
xcord1.append(dataArr[i,1])
ycord1.append(dataArr[i,2])
xcord2.append(dataArr[i,1])
ycord2.append(dataArr[i,2])
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(xcord1,ycord1,s=30,c='red',marker='s')
ax.scatter(xcord2,ycord2,s=30,c='green')
x=np.arange(-3.0,3.0,0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
需要注意的是我们，我们设置了sigmoid函数为0，0是两个分类的分界线。因此我们设定0=w0x0+w1x1+w2x2 (x0=1),所以可以解出X2与X1的关系
最后拟合的直线如图所示
随机梯度上升算法
梯度上升算法每次在更新回归系数的时候需要遍历整个数据集，如果数据量过大，计算复杂度会过高。一种改进方法时每一次只用一个样本更新回归系数，该方法成为随机梯度上升算法。由于新的样本可以不断的更新回归系数因此随机梯度上升算法是一种
随机梯度上升算法的伪代码
所有回归系数初始化为1
对数据集中每一个样本：
计算该样本的梯度
使用alpha*gradient更新回归系数
返回回归系数
下面是随机梯度上升算法的代码
def stocGradAscent0(dataMatrix,classLabels):
m,n = np.shape(dataMatrix)
alpha = 0.01
weights = np.ones(n)
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weights))
error = classLabels[i] - h
weights = weights + alpha*error *dataMatrix[i]
return weights
和梯度上升算法的区别就是变量以及误差都是数值，其他的地方差异性不大。
最后迭代100次后的效果(有100组数据)，由于迭代次数的原因，导致有一部分被分错了。
对于随机梯度算法进行改进
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
m,n = np.shape(dataMatrix)
weights=np.ones(n)
for j in range(numIter):
dataIndex = list(range(m))
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001
randIndex = int(np.random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights+alpha*error*dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
最后得到的效果相对于随机梯度上升就提升了很多。
示例：从疝气病症预测病马的死亡概率
我们有一个数据集包含368个样本和28个特征。但是有接近30%的数据是缺失的。下面介绍如何处理数据集的缺失。
准备数据：处理数据中的缺失值
数据有时候某些特征会缺失或者损坏，但是因此丢掉整个数据集代价十分昂贵，十分不可取。所以必须采用一些方法解决这个问题：
下面给出一些可选的做法：
使用可用特征均值来填补缺失值
使用特殊值来填补缺失值，例如1
忽略有缺失值的样本
使用相似样本的均值来填补缺失值
使用机器学习方法来预测缺失值
这里我们选择实数0来代替所有的缺失值，，这样恰好能够适用于logistic回归。
weights = weights + alpha * error * dataMatrix[randIndex]
如果dataMatrix的特征值为0，那么就不更新权值。
测试数据：用Logistic回归进行分类
def classifyVector(inX,weights):
prob = sigmoid(sum(inX*weights))
if prob & 0.5:
return 1.0
return 0.0
def colicTest():
frTrain = open('horseColicTraining.txt')
frTest = open('horseColicTest.txt')
trainingSet = []
trainingLabels = []
for line in frTrain.readlines():
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[21]))
trainWeights = stocGradAscent1(np.array(trainingSet),trainingLabels,500)
errorCount = 0
numTestVec = 0.0
for line in frTest.readlines():
numTestVec += 1.0
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(float(currLine[i]))
if int(classifyVector(np.array(lineArr),trainWeights)) != int(currLine[21]):
errorCount += 1
errorRate = (float(errorCount)/numTestVec)
print('the error rate of this test is %f' % errorRate)
return errorRate
def multiTest():
numTests = 10
errorSum = 0.0
for k in range(numTests):
errorSum += colicTest()
print('after %d iterations the average error rate is: %f' % (numTests,errorSum/(float(numTests))))
最后运行结果
参考资料：
机器学习实战
统计学习方法
另外看到了两篇很好的博客，便于更好的理解其中的数学原理，传送门：
我的热门文章本作业使用逻辑回归(logistic regression)和神经网络(neural networks)识别手写的阿拉伯数字(0-9)
关于逻辑回归的一个编程练习，可参考：/hapjin/p/6078530.html
下面使用逻辑回归实现多分类问题：识别手写的阿拉伯数字(0-9)，使用神经网络实现：识别手写的阿拉伯数字(0-9)，请参考：
数据加载到Matlab中的格式如下：
一共有5000个训练样本，每个训练样本是400维的列向量（20X20像素的 grayscale image），用矩阵 X 保存。样本的结果(label of training set)保存在向量 y 中，y 是一个5000行1列的列向量。
比如 y = (1,2,3,4,5,6,7,8,9,10......)T，注意，由于Matlab下标是从1开始的，故用 10 表示数字 0
①样本数据的可视化
随机选择100个样本数据，使用Matlab可视化的结果如下：
&②使用逻辑回归来实现多分类问题(one-vs-all)
所谓多分类问题，是指分类的结果为三类以上。比如，预测明天的天气结果为三类：晴(用y==1表示)、阴(用y==2表示)、雨(用y==3表示)
分类的思想，其实与逻辑回归分类(默认是指二分类，binary classification)很相似，对“晴天”进行分类时，将另外两类(阴天和下雨)视为一类：(非晴天)，这样，就把一个多分类问题转化成了二分类问题。示意图如下：（图中的圆圈 表示：不属于某一类的 所有其他类）
对于N分类问题(N&=3)，就需要N个假设函数(预测模型)，也即需要N组模型参数θ（θ一般是一个向量）
然后，对于每个样本实例，依次使用每个模型预测输出，选取输出值最大的那组模型所对应的预测结果作为最终结果。
因为模型的输出值，在sigmoid函数作用下，其实是一个概率值。，注意：hθ(1)(x)，hθ(2)(x)，hθ(3)(x)三组 模型参数θ 一般是不同的。比如：
hθ(1)(x)，输出 预测为晴天(y==1)的概率
hθ(2)(x)，输出 预测为阴天(y==2)的概率
hθ(3)(x)，输出 预测为雨天(y==3)的概率
③Matlab代码实现
对于上面的识别阿拉伯数字的问题，一共需要训练出10个逻辑回归模型，每个逻辑回归模型对应着识别其中一个数字。
我们一共有5000个样本，样本的预测结果值就是：y=(1,2,3,4,5,6,7,8,9,10)，其中 10 代表 数字0
我们使用Matlab fmincg库函数 来求解使得代价函数取最小值的 模型参数θ
function [all_theta] = oneVsAll(X, y, num_labels, lambda)
%ONEVSALL trains multiple logistic regression classifiers and returns all
%the classifiers in a matrix all_theta, where the i-th row of all_theta
%corresponds to the classifier for label i
[all_theta] = ONEVSALL(X, y, num_labels, lambda) trains num_labels
logisitc regression classifiers and returns each of these classifiers
in a matrix all_theta, where the i-th row of all_theta corresponds
to the classifier for label i
% Some useful variables
m = size(X, 1);% num of samples
n = size(X, 2);% num of features
% You need to return the following variables correctly
all_theta = zeros(num_labels, n + 1);
% Add ones to the X data matrix
X = [ones(m, 1) X];
% ====================== YOUR CODE HERE ======================
% Instructions: You should complete the following code to train num_labels
logistic regression classifiers with regularization
parameter lambda.
% Hint: theta(:) will return a column vector.
% Hint: You can use y == c to obtain a vector of 1's and 0's that tell use
whether the ground truth is true/false for this class.
% Note: For this assignment, we recommend using fmincg to optimize the cost
function. It is okay to use a for-loop (for c = 1:num_labels) to
loop over the different classes.
fmincg works similarly to fminunc, but is more efficient when we
are dealing with large number of parameters.
% Example Code for fmincg:
% Set Initial theta
initial_theta = zeros(n + 1, 1);
% Set options for fminunc
options = optimset('GradObj', 'on', 'MaxIter', 50);
% Run fmincg to obtain the optimal theta
% This function will return theta and the cost
[theta] = ...
fmincg (@(t)(lrCostFunction(t, X, (y == c), lambda)), ...
initial_theta, options);
initial_theta = zeros(n + 1, 1);
options = optimset('GradObj','on','MaxIter',50);
for c = 1:num_labels %num_labels 为逻辑回归训练器的个数，num of logistic regression classifiers
all_theta(c, :) = fmincg(@(t)(lrCostFunction(t, X, (y == c),lambda)), initial_theta,options );
% =========================================================================
lrCostFunction，完全可参考：/hapjin/p/6078530.html 里面的 正则化的逻辑回归模型实现costFunctionReg.m文件
下面来解释一下 for循环：
num_labels 为分类器个数，共10个，每个分类器(模型)用来识别10个数字中的某一个。
我们一共有5000个样本，每个样本有400中特征变量，因此：模型参数θ 向量有401个元素。
initial_theta = zeros(n + 1, 1); % 模型参数θ的初始值(n == 400)
all_theta是一个10*401的矩阵，每一行存储着一个分类器(模型)的模型参数θ 向量，执行上面for循环，就调用fmincg库函数 求出了 所有模型的参数θ 向量了。
求出了每个模型的参数向量θ，就可以用 训练好的模型来识别数字了。对于一个给定的数字输入(400个 feature variables) input instance，每个模型的假设函数hθ(i)(x) 输出一个值(i = 1,2,...10)。取这10个值中最大值那个值，作为最终的识别结果。比如g(hθ(8)(x))==0.96 比其它所有的 g(hθ(i)(x)) (i = 1,2,...10,但 i 不等于8) 都大，则识别的结果为 数字 8
function p = predictOneVsAll(all_theta, X)
%PREDICT Predict the label for a trained one-vs-all classifier. The labels
%are in the range 1..K, where K = size(all_theta, 1).
p = PREDICTONEVSALL(all_theta, X) will return a vector of predictions
for each example in the matrix X. Note that X contains the examples in
rows. all_theta is a matrix where the i-th row is a trained logistic
regression theta vector for the i-th class. You should set p to a vector
of values from 1..K (e.g., p = [1; 3; 1; 2] predicts classes 1, 3, 1, 2
for 4 examples)
m = size(X, 1);
num_labels = size(all_theta, 1);
% You need to return the following variables correctly
p = zeros(size(X, 1), 1);
% Add ones to the X data matrix
X = [ones(m, 1) X];
% ====================== YOUR CODE HERE ======================
% Instructions: Complete the following code to make predictions using
your learned logistic regression parameters (one-vs-all).
You should set p to a vector of
(from 1 to
num_labels).
% Hint: This code can be done all vectorized using the max function.
In particular, the max function can also return the index of the
max element, for more information see 'help max'. If your examples
are in rows, then, you can use max(A, [], 2) to obtain the max
for each row.
[~,p] = max( X * all_theta',[],2); % 求矩阵(X*all_theta')每行的最大值，p 记录矩阵每行的最大值的索引
% =========================================================================
阅读(...) 评论()Web Developer写在前面的废话
2014，又到了新的一年，首先祝大家新年快乐，也感谢那些关注我的博客的人。
现在想想数据挖掘课程都是去年的事了，一直预告着，盘算着年内完工的分类算法也拖了一年了。
本来打算去年就完成分类算法，如果有人看的话也顺带提提关联分析，聚类神马的，
借着新年新气象的借口来补完这一系列的文章，
可是，这明明就是在发。
尽管这个是预告里的最后一篇，但是我也没打算把这个分类算法就这么完结。尽管每一篇都很浅显，每个算法都是浅尝辄止的，在deep learning那么火的今天，掌握这些东西算起来屌丝得不能再屌丝了。考虑到一致性与完备性，最后补上两篇一样naive的：组合方法提高分类效率、几种分类方法的绩效讨论。希望读到的人喜欢。
算法六：logistic回归
& & & & 由于我们在前面已经讨论过了神经网络的分类问题（参见《》），如今再从最优化的角度来讨论logistic回归就显得有些不合适了。Logistic回归问题的最优化问题可以表述为：寻找一个非线性函数sigmoid的最佳拟合参数，求解过程可使用最优化算法完成。它可以看做是用sigmoid函数作为二阈值分类器的感知器问题。
& & & &今天我们将从统计的角度来重新考虑logistic回归问题。
一、logistic回归及其MLE
& & & &当我们考虑解释变量为分类变量如考虑一个企业是否会被并购，一个企业是否会上市，你的能否考上研究生这些问题时，考虑线性概率模型P（yi =1）= β0 + β1xi 显然是不合适的，它至少有两个致命的缺陷：1、概率估计值可能超过1，使得模型失去了意义；（要解决这个问题并不麻烦，我们将预测超过1的部分记为1，低于0的部分记为0，就可以解决。这个解决办法就是计量里有一定历史的tobit模型）2、边际效应假定为不变，通常来说不合经济学常识。考虑一个边际效应递减的模型（假定真实值为蓝线），可以看到线性模型表现很差。
& & & & 但是sigmoid函数去拟合蓝线确实十分合适的。于是我们可以考虑logistic回归模型：
& & & & &假定有N个观测样本Y1,Y2,…,YN，设P(Yi=1|Xi)=π(Xi)为给定条件Xi下得到结果Yi=1的条件概率；而在同样条件下得到结果Yi=0的条件概率为P(Yi=0|Xi)=1-π(Xi)，于是得到一个观测值的概率P（Yi）=π(Xi)Yi[1-π(Xi)]1-Yi假设各观测独立，对logistic回归模型来说，其对数似然函数为：
于是便可求解出logistic模型的MLE。
二、logit还是probit？
& & & &虽说sigmoid函数对边际递减的模型拟合良好，但是我们也要知道S型函数并非仅sigmoid函数一个，绝大多数的累积分布函数都是S型的。于是考虑F-1（P）（F为标准正态分布的累积分布函数）也不失为一个很好的选择。像这样的，对概率P做一点变换，让变换后的取值范围变得合理，且变换后我们能够有办法进行参数估计的，就涉及到广义线性模型理论中的连接函数。在广义线性模型中我们把log(P/(1-P))称为logit，F-1（P）（F为标准正态分布的累积分布函数）称为probit。那么这里就涉及到一个选择的问题：连接函数选logit还是probit？logistic回归认为二分类变量服从伯努利分布，应当选择logit，而且从解释的角度说，p/（1-p）就是我们常说的odds
ratio，也就是软件报告中出现的OR值。
& & & & & & &但是probit也有它合理的一面，首先，中心极限定理告诉我们，伯努利分布在样本够多的时候就是近似正态分布的；其次，从不确定性的角度考虑，probit认为我们的线性概率模型服从正态分布，这也是更为合理的。
& & & & &我们来看一下经过变换后，自变量和P的关系是什么样子的：
如果你确实想知道到底你的数据用哪一个方法好，也不是没有办法，你可以看一下你的残差到底是符合logit函数呢还是符合probit函数，当然，凭肉眼肯定是看不出来的，因为这两个函数本来就很接近，你可以通过函数的假定，用拟合优度检验一下。但通常，估计不会有人非要这么较真，因为没有必要。但是有一点是要注意的，logit模型较probit模型而言具有厚尾的特征，这也是为什么经济学论文爱用logit的原因。
& & & & &我们以鸢尾花数据中的virginica,versicolor两类数据分类为例，看看两种办法分类有无差别。
& & & & & & & & & & & & probit.predictions
& & & & & & & & & & & & versicolor&&& virginica
& versicolor & & & & & & & 47&&&&&&&& 3
& virginica & & & & & & & & & 3&&&&&&& 47
& & & & & & & & & & & & & & &logit.predictions
& & & & & & & & & & & & & versicolor&&& virginica
& versicolor & & & & & & & &47&&&&&&&& 3
& virginica & & & & & & & & & &3&&&&&&& 47
从上图与比较表格均可以看出，两者差别不大。
三、多项式logit与order logit
& & & &对于二项分类模型的一个自然而然的推广便是多项分类模型。
& & & &我们借鉴神经网络里提到的，有：
& & & & &按照上面这种方法，给定一个输入向量x，获得最大hθ(x)的类就是x所分到的类。
& & & & &选择最大的&hθ(x)十分好理解：在类别选择问题中，不论要选的类别是什么，每一个类别对做选择的经济个体来说都有或多或少的效用(没有效用的类别当然不会被考虑) ，一个类别的脱颖而出必然是因为该类别能产生出最高的效用。
& & & & & & 我们将多项logit模型的数学表述叙述如下：
& & & & & &假定对于第i个观测，因变量Yi有M个取值：1,2，…，M，自变量为Xi，则多项logit模型为：
& & & & & & 与logistic回归的似然估计类似，我们可以很容易写出多项logit的对数似然函数：
多项 Logit模型虽然好用，但从上面的叙述可以看出，多项 Logit 模型最大的限制在于各个类别必须是对等的，因此在可供选择的类别中，不可有主要类别和次要类别混杂在一起的情形。例如在研究旅游交通工具的选择时，可将交通工具的类别粗分为航空、火车、公用汽车、自用汽车四大类，但若将航空类别再依三家航空公司细分出三类而得到总共六个类别，则多项 Logit 模型就不适用，因为航空、火车、公用汽车、自用汽车均属同一等级的主要类别，而航空公司的区别则很明显的是较次要的类别，不应该混杂在一起。在这个例子中，主要类别和次要类别很容易分辨，但在其他的研究中可能就不是那么容易，若不慎将不同层级的类别混在一起，则由多项
Logit 模型所得到的实证结果就会有误差。
& & & & & &对于分类模型，我们还会遇到被解释变量中有分类变量的情形。对于连续变量解释离散变量，且被解释的离散变量是有顺序的（这个是和多项logit最大的区别）的情形，我们就需要考虑到order logit模型。
其数学模型叙述如下：
&其中，F(.)表示累积分布函数，当F表示正态分布的分布函数时，对应的是order probit；F表示logistic分布时，变对应order logit。
& & & &与logistic分布类似，我们可以很容易写出其对数似然函数：
四、dummy variable
& & & & & & 在logistic回归中，经常会遇到解释变量为分类变量的情形，比如收入：高、中、低；地域：北京、上海、广州等。这里对分类变量而言就涉及一个问题：要不要将分类变量设置dummy variable？
& & & & & 这个问题的答案在线性模型中很显然，必须要这么做！！！如果我们不设置哑变量，而是单纯地赋值：北京=1，上海=2，广州=3，即我们将自变量视作连续性的数值变量，但这仅仅是一个代码而己，并不意味着地域间存在大小次序的关系，即并非代表被解释变量（响应变量）会按此顺序线性增加或减少。即使是有序多分类变量，如家庭收入分为高、中、低三档，各类别间的差距也是无法准确衡量的，按编码数值来分析实际上就是强行规定为等距，这显然可能引起更大的误差。
& & &但是在logistic回归中，由于logit（p)变化的特殊性，在解释定序变量时，为了减少自由度（即解释变量个数），我们常常将定序变量（如家庭收入分为高、中、低）视为连续的数值变量，而且经济解释可以是XX每提高一个档次，相应的概率会提高expression（delta(XX))（expression的表达式还是很复杂的，不打了）。当然减少变量个数是以牺牲预测精度为代价的。毕竟数据处理是一门艺术而非一门技术，如何取舍还得具体问题具体分析。当然，非定序的分类变量是万万不可将其视为数值变量的。
五、广义线性模型的R实现
& & & & & R语言提供了广义线性模型的拟合函数glm()，其调用格式如下：
glm(formula, family = gaussian, data,weights, subset,
&&&na.action, start = NULL, etastart, mustart, offset,
&&& control= list(...), model = TRUE, method = &glm.fit&,
&&& x =FALSE, y = TRUE, contrasts = NULL, ...)
参数说明：
Formula：回归形式，与lm()函数的formula参数用法一致
Family：设置广义线性模型连接函数的典则分布族，glm()提供正态、指数、gamma、逆高斯、Poisson、二项分布。我们的logistic回归使用的是二项分布族binomial。Binomial族默认连接函数为logit，可设置为probit。
Data：数据集
& & & & 鸢尾花例子使用的R代码：
logit.fit &- glm(Species~Petal.Width+Petal.Length,
family = binomial(link = 'logit'),
data = iris[51:150,])
logit.predictions &- ifelse(predict(logit.fit) & 0,'virginica', 'versicolor')
table(iris[51:150,5],logit.predictions)
probit.fit &- glm(Species~Petal.Width+Petal.Length,
family = quasibinomial(link = 'probit'),
data = iris[51:150,])
probit.predictions &- ifelse(predict(probit.fit) &0,'virginica', 'versicolor')
table(iris[51:150,5],probit.predictions)
& & & & 程序包mlogit提供了多项logit的模型拟合函数：
mlogit(formula, data, subset, weights,na.action, start = NULL,
&&&&&&alt.subset = NULL, reflevel = NULL,
&&&&&&nests = NULL, un.nest.el = FALSE, unscaled = FALSE,
&&&&&&heterosc = FALSE, rpar = NULL, probit = FALSE,
&&&&&&R = 40, correlation = FALSE, halton = NULL,
&&&&&&random.nb = NULL, panel = FALSE, estimate = TRUE,
&&&&&&seed = 10, ...)
mlogit.data(data, choice, shape = c(&wide&,&long&), varying = NULL,
&&&&&&&&&&& sep=&.&,alt.var = NULL, chid.var = NULL, alt.levels = NULL,id.var = NULL, opposite = NULL, drop.index = FALSE,
&&&&&&&&&&& ranked = FALSE, ...)
参数说明：
formula：mlogit提供了条件logit，多项logit，混合logit多种模型，对于多项logit的估计模型应写为：因变量~0|自变量,如：mode ~ 0 | income
data：使用mlogit.data函数使得数据结构符合mlogit函数要求。
Choice：确定分类变量是什么
Shape：如果每一行是一个观测，我们选择wide，如果每一行是表示一个选择，那么就应该选择long。
alt.var：对于shape为long的数据，需要标明所有的选择名称
&&&&&&& 选择wide的数据示例：
&&&&&&&&&& 选择long的数据示例：
以fishing数据为例，来说明如何使用mlogit。
library(mlogit)
data(&Fishing&, package = &mlogit&)
Fish &- mlogit.data(Fishing,shape = &wide&,choice = &mode&)
summary(mlogit(mode ~ 0 | income, data = Fish))
这个输出的结果与nnet包中的multinom()函数一致。由于mlogit包可以做的logit模型更多，所以这里就不在对nnet包的multinom作介绍了，可以参见《根据Econometrics in R一书，将回归方法总结一下》一文。
程序包MASS提供polr()函数可以进行ordered logit或probit回归。用法如下：
polr(formula, data, weights, start, ..., subset, na.action,
&&&& contrasts = NULL, Hess = FALSE, model = TRUE,
&&&& method = c(&logistic&, &probit&, &cloglog&, &cauchit&))
参数说明：
Formula：回归形式，与lm()函数的formula参数用法一致
Data：数据集
Method：默认为order logit，选择probit时变为order probit模型。
以housing数据为例说明函数用法：
house.plr &- polr(Sat ~ Infl + Type + Cont, weights = Freq, data = housing)
summary(house.plr, digits = 3)
这些结果十分直观，易于解读，所以我们在这里省略所有的输出结果。
再看手写数字案例：
& & & & & &最后，我们回到最开始的那个手写数字的案例，我们试着利用多项logit重做这个案例。（这个案例的描述与数据参见《》一章）
& & & & & &特征的选择可参见《》一章。
由于手写数字的特征选取很容易导致回归系数矩阵是降秩的，所以我们使用nnet包的multinom()函数代替mlogit()。
运行下列代码：
setwd(&D:/R/data/digits/trainingDigits&)
names&-list.files(&D:/R/data/digits/trainingDigits&)
data&-paste(&train&,1:1934,sep=&&)
for(i in 1:length(names))
assign(data[i],as.matrix(read.fwf(names[i],widths=rep(1,32))))
label&-factor(rep(0:9,c(189,198,195,199,186,187,195,201,180,204)))
feature&-matrix(rep(0,length(names)*25),length(names),25)
for(i in 1:length(names)){
feature[i,1]&-sum(get(data[i])[,16])
feature[i,2]&-sum(get(data[i])[,8])
feature[i,3]&-sum(get(data[i])[,24])
feature[i,4]&-sum(get(data[i])[16,])
feature[i,5]&-sum(get(data[i])[11,])
feature[i,6]&-sum(get(data[i])[21,])
feature[i,7]&-sum(diag(get(data[i])))
feature[i,8]&-sum(diag(get(data[i])[,32:1]))
feature[i,9]&-sum((get(data[i])[17:32,17:32]))
feature[i,10]&-sum((get(data[i])[1:8,1:8]))
feature[i,11]&-sum((get(data[i])[9:16,1:8]))
feature[i,12]&-sum((get(data[i])[17:24,1:8]))
feature[i,13]&-sum((get(data[i])[25:32,1:8]))
feature[i,14]&-sum((get(data[i])[1:8,9:16]))
feature[i,15]&-sum((get(data[i])[9:16,9:16]))
feature[i,16]&-sum((get(data[i])[17:24,9:16]))
feature[i,17]&-sum((get(data[i])[25:32,9:16]))
feature[i,18]&-sum((get(data[i])[1:8,17:24]))
feature[i,19]&-sum((get(data[i])[9:16,17:24]))
feature[i,20]&-sum((get(data[i])[17:24,17:24]))
feature[i,21]&-sum((get(data[i])[25:32,17:24]))
feature[i,22]&-sum((get(data[i])[1:8,25:32]))
feature[i,23]&-sum((get(data[i])[9:16,25:32]))
feature[i,24]&-sum((get(data[i])[17:24,25:32]))
feature[i,25]&-sum((get(data[i])[25:32,25:32]))
data1 &- data.frame(feature,label)
#降秩时mlogit不可用
#data10&- mlogit.data(data1,shape = &wide&,choice = &label&)
#m1&-mlogit(label~0|X1+X2+X3+X4+X5+X6+X7+X8+X9+X10+X11+X12+X13+X14+X15+X16+X17+X18+X19+X20+X21+X22+X23+X24+X25,data=data10)
library(nnet)
m1&-multinom(label ~ ., data = data1)
pred&-predict(m1,data1)
table(pred,label)
sum(diag(table(pred,label)))/length(names)
setwd(&D:/R/data/digits/testDigits&)
name&-list.files(&D:/R/data/digits/testDigits&)
data1&-paste(&train&,1:1934,sep=&&)
for(i in 1:length(name))
assign(data1[i],as.matrix(read.fwf(name[i],widths=rep(1,32))))
feature&-matrix(rep(0,length(name)*25),length(name),25)
for(i in 1:length(name)){
feature[i,1]&-sum(get(data1[i])[,16])
feature[i,2]&-sum(get(data1[i])[,8])
feature[i,3]&-sum(get(data1[i])[,24])
feature[i,4]&-sum(get(data1[i])[16,])
feature[i,5]&-sum(get(data1[i])[11,])
feature[i,6]&-sum(get(data1[i])[21,])
feature[i,7]&-sum(diag(get(data1[i])))
feature[i,8]&-sum(diag(get(data1[i])[,32:1]))
feature[i,9]&-sum((get(data1[i])[17:32,17:32]))
feature[i,10]&-sum((get(data1[i])[1:8,1:8]))
feature[i,11]&-sum((get(data1[i])[9:16,1:8]))
feature[i,12]&-sum((get(data1[i])[17:24,1:8]))
feature[i,13]&-sum((get(data1[i])[25:32,1:8]))
feature[i,14]&-sum((get(data1[i])[1:8,9:16]))
feature[i,15]&-sum((get(data1[i])[9:16,9:16]))
feature[i,16]&-sum((get(data1[i])[17:24,9:16]))
feature[i,17]&-sum((get(data1[i])[25:32,9:16]))
feature[i,18]&-sum((get(data1[i])[1:8,17:24]))
feature[i,19]&-sum((get(data1[i])[9:16,17:24]))
feature[i,20]&-sum((get(data1[i])[17:24,17:24]))
feature[i,21]&-sum((get(data1[i])[25:32,17:24]))
feature[i,22]&-sum((get(data1[i])[1:8,25:32]))
feature[i,23]&-sum((get(data1[i])[9:16,25:32]))
feature[i,24]&-sum((get(data1[i])[17:24,25:32]))
feature[i,25]&-sum((get(data1[i])[25:32,25:32]))
labeltest&-factor(rep(0:9,c(87,97,92,85,114,108,87,96,91,89)))
data2&-data.frame(feature,labeltest)
pred1&-predict(m1,data2)
table(pred1,labeltest)
sum(diag(table(pred1,labeltest)))/length(name)
经整理，输出结果如下：（左边为训练集，右边为测试集）
Tips: oddsratio=p/1-p 相对风险指数贝努力模型中 P是发生A事件的概率，1-p是不发生A事件的概率所以p/1-p是 发生与不发生的相对风险。OR值等于1，表示该因素对A事件发生不起作用；OR值大于1，表示A事件发生的可能性大于不发生的可能性；OR值小于1，表示A事件不发生的可能性大于发生的可能性。
Further reading：
Yves Croissant：
本文之后待写的文章有：
算法提升：组合方法提高分类效率
算法评价：几种分类方法的绩效讨论
（to be continue)
本文已收录于以下专栏：
相关文章推荐
error theta
#虽然是python 标签 但缺少
sigmoid &- function(x){
return (1.0/(1+exp(-x)))
#用什么样的优化算...
虽说线性回归无法直接用于分类预测，但可以对其加层映射：将连续无穷输出映射到指定的有限输出。逻辑回归(Logistic Regression, LR)便是基于此思想在线性回归的结果上加上一个逻辑函数，将...
本节将一下逻辑回归和R语言实现，逻辑回归(LR，LogisticRegression)其实属于广义回归模型，根据因变量的类型和服从的分布可以分为，普通多元线性回归模型，和逻辑回归，逻辑回归是指因变量是...
在感知器模型、神经网络模型、深度学习模型中均会看见激活函数的声影。激活函数又被称为转移函数、激励函数、传输函数或限幅函数，其作用就是将可能的无限域变换到一指定的有限范围内输出，这类似于生物神经元具有的...
直接上代码，如下：
data_sample &- iris[51:150,];
m &- dim(data_sample)[1]
#获取数据集记录条数
val &- sample(m, size ...
线性判别分析(Linear Discriminant Analysis, LDA)，也叫 做Fisher线性判别(Fisher Linear Discriminant ,FLD)，是模式识别的经典算法...
模拟退火其实也是一种贪心算法，但是它的搜索过程引入了随机因素。在迭代更新可行解时，以一定的概率来接受一个比当前解要差的解，因此有可能会跳出这个局部的最优解，达到全局的最优解。
kNN算法是k近邻分类(k-nearest neighbor classification)算法的简称。基本流程是从训练集中找到和新数据最接近的k条记录，然后根据他们的主要分类来决定新数据的类别。该算...
逻辑回归是啥？Logistic 回归是一个二分类算法，用来预测给定独立变量集的二分类输出。我们使用哑变量代替二分类输出。也可以把逻辑回归看成输出为类别变量的特殊的线性回归(使用对数几率作为依赖变量)。...
他的最新文章
讲师：王哲涵
讲师：韦玮
您举报文章：
举报原因：
原文地址：
原因补充：
(最多只允许输入30个字)}