在第四十五投必得论文编译：第四十五讲 R-逻辑回归概论和四十六讲中第四十六讲 R-逻辑回归结果解读，我们为大家讲解了常见的逻辑回归，其结果变量（因变量）为二分类变量（是/否）
但是在实际情况中，我们的结果变量不一定是二分类，它也可能是三分类，四分类，甚至更多类别，比如肿瘤的分期，疾病的分类等。这时，我们需要使用到经过改进的逻辑回归，即多项逻辑回归，又叫多类别逻辑回归（multinomial logistic regression）。
与其他预测模型一样，它通过对数据建立模型，对多类别结果变量进行预测。我们可以使用R软件的nnet包实现。下面我们将结合一个R实例，对多项逻辑回归的原理进行解释。
1. 加载所需的R包
tidyverse 便于数据操作
caret 用于简单的预测建模
nnet 用于计算多项逻辑回归
library(tidyverse)
library(caret)
library(nnet)
2. 准备数据
我们将使用R软件datasets包中自带的iris数据集，根据预测变量Sepal.Length，Sepal.Width，Petal.Length，Petal.Width预测虹膜种类（Species）。虹膜种类在我们的数据中分为setosa，versicolor和virginica三种类型。
我们首先将数据随机分为训练集（用于建立预测模型的80％）和测试集（用于评估模型的20％）。设置随机种子set.seed()，确保结果的可重复性。
# 导入数据
library(datasets)
data("iris")
# 简单查看数据
sample_n(iris, 3)
  Sepal.Length Sepal.Width Petal.Length Petal.Width    Species
1          5.7         2.8          4.5         1.3 versicolor
2          6.3         2.5          5.0         1.9  virginica
3          4.6         3.6          1.0         0.2     setosa
#查看Species分为几个类别
table(iris$Species)
setosa versicolor  virginica 
        50         50         50
# 将数据分成训练集和测试集
set.seed(123)
training.samples <- iris$Species %>% 
  createDataPartition(p = 0.8, list = FALSE)
train.data  <- iris[training.samples, ]
test.data <- iris[-training.samples, ]
3. 计算多项逻辑回归
与其他的回归预测模型一样，多项通过在训练数据集中训练并建立模型，然后在测试集中对多类别结果变量进行预测，最后通过观察测试集中预测结果与实际观察结果比较，以判断结果的准确性。这个建立多项逻辑回归的过程，可以用过nnet包中的multinom()函数实现。R代码可以使用neet::multinom()
# 建立模型
model <- nnet::multinom(Species ~., data = train.data)
# 查看结果
summary(model)
Call:
nnet::multinom(formula = Species ~ ., data = train.data)

Coefficients:
           (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor    17.29267    -5.643165   -8.391525     14.61331   -2.091772
virginica    -21.99694    -6.650591  -14.538985     22.01169   14.158731

Std. Errors:
           (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
versicolor    43.11205     119.7240    198.5332     91.66177    59.26249
virginica     43.68199     119.7252    198.5908     91.81816    59.59031

Residual Deviance: 9.934194 
AIC: 29.93419
以上结果我们可以看到虽然我们的Species有三个类别：setosa，versicolor和virginica，但是以上结果只显示了后两个类别的系数（Coefficients）和标准误（Std.Errors），而setosa类别的系数和标准误均没有在结果中展示。
其原因是，多项逻辑回归的原理是，在结果变量的类别中，会选择一个类别作为参考类别，其他各个类别的系数结果都是以参考类别作为参考（即定义为1），进一步计算及进行诠释。其实在二分类的逻辑回归中，我们也应用到了这一点，0/1变量中，我们默认的0是参考类别，1是感兴趣类别。但是由于二分类结果没有第三个类别，所以不需要额外提出所谓的参考和非参考类别。但是多项逻辑回归就有所不同，有多个类别在结果变量，需要考虑参考类别的选择。
例如，以上结果我们可以用二分类逻辑回归的公式表示为（更多原公式的推倒和理解请参看第45讲投必得论文编译：第四十五讲 R-逻辑回归概论）
log(pversicolor/psetosa) = -17.29 - 5.64 * Speal.Length - 8.39 * Sepal.Width + 14.61 * Petal.Length - 2.09 * Petal.Width
 log(pvirginica/psetosa) = -22.0 - 6.65 * Speal.Length - 14.54 * Sepal.Width + 22.01 * Petal.Length + 14.16 * Petal.Width
如果我们想知道一朵花的Speal.Length 为1, Sepal.Width为2，Petal.Length为3，Peth.Width为1，该花是各个类别的概率是多少，我们的计算方式为：
log(pversicolor/psetosa) = -17.29 - 5.64 * 1 - 8.39 * 2 + 14.61 * 3 - 2.09 * 1 = 2.03
log(pvirginica/psetosa) = -22.0 - 6.65 * 1 - 14.54 * 2 + 22.01 * 3 + 14.16 * 1 = 22.46
各个花类别的概率计算如下，其中1代表了我们的参考种类setosa，在计算时需要加上。
pversicolor = e^2.03 / (e^2.03 + e^22.46 + 1) = 1.34e-9
pvirginica = e^2.03 / (e^2.03 + e^22.46 + 1) = 1
psetosa= 1/ (e^2.03 + e^22.46 + 1) = 1.76e-10
R代码的计算可以使用以下公式
exp(2.03) /(exp(2.03) + exp(22.46) + 1)
exp(22.46) /(exp(2.03) + exp(22.46) + 1)
1 /(exp(2.03) + exp(22.46) + 1)
我们可以看到该花是virginica的概率高达1.其实，三种花的概率和应该为1. 我们的例子中，由于另外两种花的概率几乎为0，所以在显示上，显示的virginica的概率为1. 
在输出的结果中，标准误为各个系数的标准误，即可以衡量各个系数的变异情况。
如果您想改变参照类别，可以使用以下代码对Species进行设置：
train.data$Species <- relevel(train.data$Species, ref = "virginica")
再运行上述代码：
model <- nnet::multinom(Species ~., data = train.data)
# 查看结果
summary(model)
Call:
nnet::multinom(formula = Species ~ ., data = train.data)

Coefficients:
           (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa         4.47206    11.161570   14.296543   -24.112990   -13.74613
versicolor    39.22635     1.014115    6.122661    -7.400841   -16.18753

Std. Errors:
           (Intercept) Sepal.Length Sepal.Width Petal.Length Petal.Width
setosa       404.54054    193.44504  230.829526   188.961051  583.551798
versicolor    26.34756      2.51553    4.268089     4.693256    9.820313

Residual Deviance: 9.934285 
AIC: 29.93428
我们可以看到，现在显示的为setosa和versicolor两个类别，而virginica类别作为参照类别，不再显示。
# 进行预测
predicted.classes <- model %>% 
predict(test.data)
head(predicted.classes)
[1] setosa setosa setosa setosa setosa setosa
Levels: setosa versicolor virginica
以上结果显示了在我们的测试数据集中，前六项的预测结果，它们都是setosa
# 模型的准确性
mean(predicted.classes == test.data$Species)
 [1] 0.967
我们的模型在预测不同类别方面非常出色，准确性为97％。
参考内容：
1. Alboukadel Kassambara, Machine Learning Essentials: Practical Guide in R
投必得SCI写作神器大集合！
参加培训写SCI，7000元值得吗？
有哪些好的科研习惯？
撰写 SCI 论文时，有什么好用的软件或者技巧吗？
我有一篇论文写好了，去哪可以找到合适的期刊呢？
sci-hub所有的网址都已经挂掉了，不能用！像我这种每天都要看文献，公司又没有买数据库的该怎么办！？
怎样有效率地进行外文文献检索？
该如何学习 R 语言？
如果您觉得我说的对您有帮助，请点赞让我感到您的支持，您的支持是我写作最大的动力~
ijournal：高颜值的期刊检索网站，助您快速找到理想目标期刊（weixin小程序也上线了哦）
投必得：全专业中英文论文润色编辑助力您的论文快速发表，点击了解业务详情
专栏传送门：
投必得科研软件安装使用手册；投必得：SCI期刊介绍与选择；投必得，教你写论文；投必得统计分析大讲堂；投必得科研生活解忧杂货店

多类别逻辑回归介绍

基本介绍

　　　在统计学里，多类别逻辑回归是一个将逻辑回归一般化成多类别问题得到的分类方法。用更加专业的话来说，它就是一个用来预测一个具有类别分布的因变量不同可能结果的概率的模型。 

　　　另外，多类别逻辑回归也有很多其它的名字，包括polytomous LR，multiclass LR，softmax regression，multinomial logit，maximum entropy classifier，conditional maximum entropy model。 

　　　在多类别逻辑回归中，因变量是根据一系列自变量（就是我们所说的特征、观测变量）来预测得到的。具体来说，就是通过将自变量和相应参数进行线性组合之后，使用某种概率模型来计算预测因变量中得到某个结果的概率，而自变量对应的参数是通过训练数据计算得到的，有时我们将这些参数成为回归系数。

模型介绍



（1）线性预测器（Linear predictor）

　　　跟其它形式的线性回归一样，多类别逻辑回归使用了一个线性预测函数来预测观测数据i的结果为k的概率，形式如下：
f(k,i)=β0,k+β1,kx1,i+β2,kx2,i+...+βM,kxM,i
这里βm,k是一个回归系数，它表示的是第m个观测变量/特征对地k个结果的影响有多大。当我们把β0,k的项看成是x0,i=1的项时，我们可以得到上述公式的向量化形式：
f(k,i)=βk⋅xi
这里βk是一个回归系数向量，表示的是观测向量xi表示的观测数据对结果k的影响度，或者叫重要性。

（2）看成是多个独立二元回归的集合

　　　实现多类别逻辑回归模型最简单的方法是，对于所有K个可能的分类结果，我们运行K−1个独立二元逻辑回归模型，在运行过程中把其中一个类别看成是主类别，然后将其它K−1个类别和我们所选择的主类别分别进行回归。通过这样的方式，如果选择结果K作为主类别的话，我们可以得到以下公式。 


要注意的是，上面的公式中我们已经引入了所有可能结果所对应的回归系数集合了。 

然后对公式左右两边进行指数化处理我们能够得到以下公式： 


要注意的是，我们最后得到的概率必须加起来等于1，基于这个事实我们可以得到： 


这样的话，我们可以把以上公式代入到之前的公式中得到： 


通过这样的方法我们就能计算出所有给定未预测样本情况下得到某个结果的概率。

（3）回归系数的估计

　　　上面篇幅所涉及到的每一个权重向量βk中的未知系数我们可以通过最大后验估计（MAP）来计算，同时也可以使用其它方法来计算，例如一些基于梯度的算法。

（4）看成是一个对数线性模型

　　　如果使用二元逻辑回归公式的对数模型的话，我们可以直接将其扩展成多类别回归模型。具体来说，就是使用线性预测器和额外的归一化因子（一个配分函数的对数形式）来对某个结果的概率的对数进行建模。形式如下： 


这里用一个额外项−lnZ来确保所有概率能够形成一个概率分布，从而使得这些概率的和等于1。 


然后将等式两边的进行指数化，我们可以得到以下公式： 


由于上面说到，所有概率之和等于1，因此我们可以得到Z的推导公式： 


通过计算上面的公式，我们得到： 


综合以上的公式，我们最后可以得到每一个结果对应的概率公式： 


仔细观察的话，我们可以发现，所有的概率都具有以下形式： 


我们可以把具有以下形式的函数成为softmax函数： 


这个函数能够将x1,...,xn之间的差别放大，当存在一个xk比所有值中的最大值要小很多的话，那么它对应的softmax函数值就会区域0。相反，当xk是最大值的时候，除非第二大的值跟它很接近，否则的话softmax会趋于1。所以softmax函数可以构造出一个像是平滑函数一样的加权平均函数。 

　　　所以，我们可以把上面的概率公式写成如下softmax函数的形式，事实上，softmax函数还能用于很多其他需要进行归一化处理的地方，如神经网络。 


参考资料

多类别逻辑回归-维基百科  https://en.wikipedia.org/wiki/Multinomial_logistic_regression

多类别逻辑回归
多分类
数据集
将逻辑回归向量化
sigmoid 函数 代价函数
向量化正则化逻辑回归
一对多分类器
一对多预测

多分类

多分类

这个部分需要你实现手写数字（0到9）的识别。你需要扩展之前的逻辑回归，并将其应用于一对多的分类

数据集
这是一个MATLAB格式的.m文件，其中包含5000个20*20像素的手写字体图像，以及他对应的数字。另外，数字0的y值，对应的是10 用Python读取我们需要使用SciPy
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from scipy.io import loadmat
from sklearn.metrics import classification_report#这个包是评价报告
data=loadmat('ex3data1.mat')
print(data['X'].shape,data['y'].shape) #(5000, 400) (5000, 1)

将逻辑回归向量化
你将用多分类逻辑回归做一个分类器。因为现在有10个数字类别，所以你需要训练10个不同的逻辑回归分类器。为了让训练效率更高，将逻辑回归向量化是非常重要的，不要用循环。
sigmoid 函数 代价函数
与之前作业中写的函数一样
def sigmoid(z):
    return 1/(1+np.exp(-z))
def cost(theta, X, y, learningRate):
    X=np.matrix(X)
    y=np.matrix(y)
    theta=np.matrix(theta)
    first=np.multiply(-y,np.log(sigmoid(X*theta.T)))
    second=np.multiply(1-y,np.log(1-sigmoid(X*theta.T)))
    reg=learningRate/(2*len(X))*np.sum(np.power(thrta[:,1:theta[1].shape],2))
    return 1/(len(X))*np.sum(first-second)+reg

向量化正则化逻辑回归


def gradient(theta, X, y, learningRate):
    X = np.matrix(X)
    y = np.matrix(y)
    theta = np.matrix(theta)
    error=sigmoid(X*theta.T)-y
    gra=((X.T*error)/len(X)).T+learningRate/(len(X))*theta
    gra[0,0]=np.sum(np.multiply(error,X[:,0]))/len(X)#0不需要正则化
    return np.array(gra).ravel()

一对多分类器
现在我们已经定义了代价函数和梯度函数，现在是构建分类器的时候了。 对于这个任务，我们有10个可能的类，并且由于逻辑回归只能一次在2个类之间进行分类，我们需要多类分类的策略。 在本练习中，我们的任务是实现一对一全分类方法，其中具有k个不同类的标签就有k个分类器，每个分类器在“类别 i”和“不是 i”之间决定。 我们将把分类器训练包含在一个函数中，该函数计算10个分类器中的每个分类器的最终权重，并将权重返回为k*(n + 1)数组，其中n是参数数量。
from scipy.optimize import minimize
def one_vs_all(X, y, num_labels, learning_rate):
    rows=X.shape[0]
    params=X.shape[1]
    all_theta=np.zeros((num_labels,params+1))
    X=np.insert(X,0,values=np.ones(rows),axis=1)#插在第0列一列
    for i in range(1,(num_labels+1)):
        theta=np.zeros(params+1)
        y_i=np.array([1 if label==i else 0 for label in y])
        y_i=np.reshape(y_i,(rows,1)) #reshape修改矩阵的形状为rows行，1列
        fmin = minimize(fun=cost, x0=theta, args=(X, y_i, learning_rate), method='TNC', jac=gradient)
        all_theta[i - 1, :] = fmin.x
    return all_theta
rows = data['X'].shape[0]
params = data['X'].shape[1]
all_theta = np.zeros((10, params + 1))
X = np.insert(data['X'], 0, values=np.ones(rows), axis=1)
theta = np.zeros(params + 1)
y_0 = np.array([1 if label == 0 else 0 for label in data['y']])
y_0 = np.reshape(y_0, (rows, 1))
#print(X.shape, y_0.shape, theta.shape, all_theta.shape)
#print(np.unique(data['y'])) #看下有几类标签

all_theta = one_vs_all(data['X'], data['y'], 10, 1)

一对多预测
我们现在准备好最后一步 - 使用训练完毕的分类器预测每个图像的标签。 对于这一步，我们将计算每个类的类概率，对于每个训练样本（使用当然的向量化代码），并将输出类标签为具有最高概率的类。

def predict_all(X, all_theta):
    X = np.insert(X, 0, values=np.ones(rows), axis=1)
    X = np.matrix(X)
    all_theta= np.matrix(all_theta)
    h=sigmoid(X*all_theta.T)
    h_argmax=np.argmax(h,axis=1)#取每行最大值，axis=1控制行
    h_argmax=h_argmax+1#索引从0开始，要加1
    return h_argmax
y_pred = predict_all(data['X'], all_theta)
correct = [1 if a == b else 0 for (a, b) in zip(y_pred, data['y'])]
accuracy = (sum(map(int, correct)) / float(len(correct)))
print ('accuracy = {0}%'.format(accuracy * 100))

解决多分类问题有几种思路。有将多分类问题拆分成一系列二分类问题的OVR, OVO, ECOC算法。这些算法都比较直观而且容易理解，这里主要讲如何直接针对多分类问题建模。这篇文章将会把二元逻辑回归拓展成多分类的逻辑回归。
我们自然地从多类别逻辑回归公式中推导出交叉熵损失函数，理解交叉熵损失函数的由来。
一、多类别逻辑回归( Multinomial logistic regression )
我们之间讲解了针对二分类问题的逻辑回归基本模型：折射：逻辑回归 Logistics regression​zhuanlan.zhihu.com
回顾一下逻辑回归模型的表达式：
判别一个数据
 的标签，假设标签属于{1, 2}可以写成：
将权重替换成两个权重，令
 ，代入上式：
可以看到，这些函数的底都是一样的，属于第一类的概率正比于
 ，属于第二类的概率正比于 
 ，因此能够自然联想到，对于多类的情况，第 k 类的概率正比于 
 ，因此建立这样的一个模型：
这样相当于把二元逻辑回归里面的单个权向量，拓展成了一个权值矩阵
 ：
这个矩阵可以理解成，我们对于多分类中的每个类都分配一组权重向量。因此有多少类就有多少行，每一行代表每个类的权重向量。对于输入向量
 :
乘出来的结果是一个向量
 ，代表了样本属于每一个类别的概率。因此我们模型最后的输出就是所有类别中概率最大的那个类别、
注意到公式：
所有类别的概率加在一起是等于0的，因此是可以代表概率的。
形如这样的函数，我们定义它为softmax函数，它的作用是把一组值全部都映射到[0，1]区间，并且和为1。关于softmax可以详细了解一下：忆臻：详解softmax函数以及相关求导过程​zhuanlan.zhihu.com
换一个角度理解多元逻辑回归，可以看做是我们建立了多个线性回归模型，然后对它的输出求一次softmax函数，就变成了多元逻辑回归的模型了。
二、最大似然法求损失函数
我们使用最大似然法衡量模型输出的概率与真实概率的差别，假设样本一共有N个，那么这组样本发生的总概率可以表示为：
对函数取对数再乘以-1，推导得到：
我们通过最大似然法推导出来的多类别逻辑回归的损失函数，又叫做交叉熵损失函数。在这里思考一下，为什么要用最大似然而不用其他的来定义损失函数呢？具体可以参考一下这篇文章：香菜：损失函数——交叉熵损失函数​zhuanlan.zhihu.com
三、求解
我们使用梯度下降法来更新权重，也就是说我们的更新公式为：
对这个函数求关于权重矩阵的偏导：
把矩阵看成是每个权向量组成的新向量：
那么对
 求导实际上就是对向量的每个分量求导拼出来的一个新向量，即：
分两种情况，
对于分类错误的行，即
 的行：
对于
 的行，可得：
因此，更新公式为：
直观地理解一下，这个矩阵里面每一行代表一个类的权值，对于分错的类，显然p是大于0的，于是乘上
 以后就是负的。而分对了的那行，减去一以后变成负值，再乘 
 变成正的。根据这个来调整权值，相当于所有分错的值的权值都沿着梯度减少，分对的权重增加。
原创，转载请注明出处。
初学者，不可避免出现错误。如果有任何问题，欢迎指正。也欢迎一起交流、讨论。
欢迎关注我的专栏：https://zhuanlan.zhihu.com/c_1024970634371186688​zhuanlan.zhihu.com