因变量为无序多分类资料，或者因变量虽为有序多分类但不满足比例优势假定(平行性检验P>0.05)，可采用无序多分类的logistic回归进行分析。当然当结局变量无序，自变量只有一个且为分类变量时，可以直接采用卡方检验；结局变量有序，自变量只有一个且为分类变量时，可以直接采用非参数检验。
无序多分类的logistic回归模型与有序多分类的logistic回归模型不同。有序多分类的logistic回归采用的是累积logit模型，进行logit变换的是因变量有序取值水平的累积概率；而无序多分类的logistic回归采用的是广义logit模型，是用因变量的各个水平(除参照水平外)与参照水平比值的自然对数来建立模型方程，当水平数为2时，该模型等价于二分类资料的logistic回归，因此该模型可以看做是二分类logistic回归模型的扩展。因变量y是有n个水平的无序多分类变量，进行无序多分类的logistic回归时可以产生n-1个广义logit模型。参照水平R的阳性概率记为πR，第k个水平(k=1,2,…n)的阳性概率分别为πk，则有π1+π2+…πn=1。自变量x有m个，第k个水平第i个自变量(i=1,2，…m)Xi系数为βki。
以因变量为4水平(第4水平为参照水平)，自变量有m个为例，模型可以表示为：
显然π1+π2+π3+π4=1，如希望比较1和2，可以将相应的两个公式相减即可得到相应的函数，同理可以比较1和3，或者2和3。当然我们也可以直接对参照水平进行修改。
示例：某研究人员欲了解不同社区和性别之间成年居民获取健康知识途径是否不同，对2个社区的314名成人进行了调查，结果见下表。变量赋值为：社区(社区A=0，社区B=1)、性别(男=0，女=1)、获取健康知识途径(传统大众媒介=1，网络=2，社区宣传=3)。请拟合社区和性别对居民获取健康知识途径的多分类logistic回归模型。
1、数据录入
2、数据加权：Data>>Weight Cases…，将[频数]进行加权
3、多变量回归分析：Analyze>>Regression>>Multinomial Logistic…
l 因变量：获得途径
l 因素：社区、性别
因变量和因素必须是分类变量，协变量是试验设计中不为研究者关注但对结果会产生影响的独立解释变量，可以分类变量也可以是连续型变量。在[因变量]下方的[参考类别(Reference Category…)]中可以对参考类别和类别顺序进行设置，默认的参考类别是最后一个类别，默认的类别顺序为升序。升序时因变量取值最小的为第一类别，而降序则将取值最小的为最后一类。
【模型】：可以指定分析的模型，默认是只分析主效应，也可以进行全因素分析(主效应+交互作用)，当然也可以进行自定义分析。选中自定义/逐步(Custom/Stepwise)后，除了可以自定义模型外，还可以实现变量的筛选，类似于二分类logistics回归中的Block和Method。本例取默认的主效应分析。
【统计量】:除默认选项，同时选中信息准则(输出AIC和BIC)、单元格概率、分类表及拟合优度检验。定义子总体默认选项为对所有的自变量和协变量计算单元格概率并进行拟合优度检验。
【收敛准则】：主要对迭代进行设置。
【选项】：可对进入和剔除标准及其检验方法进行设置。【保存】：可以保存新变量[估计反应概率]、[预测分类]、[预测分类概率]、[实际分类概率]。4、结果
【案例处理概要】：分析示例的基本情况。
【模型拟合信息】：与只含常数项的初始模型相比，最终模型的AIC(Akaike信息准则)、BIC(贝叶斯信息准则)、负2倍的对数似然值(-2LL)均有下降。-2LL值从80.877下降至36.821，下降了44.056(卡方值)，似然比卡方检验具有统计学意义(P<0.001)，说明模型纳入性别和社区两个变量中至少有一个偏回归系数不为0。
【拟合优度检验】：显示Pearson拟合优度检验和Deviance拟合优度检验结果。此两种方法实际是检验当前模型预测值与样本实测值的比较，两者结果P值均大于0.05，表示拟合良好。但要注意这两种方法对自变量的样本量有一定要求，当自变量较多或者含有连续变量时一般不采用这两种方法的检验结果。
【伪R2】：输出三种伪决定系数。对分类数据的统计分析，不需要过于在意这三种伪决定系数过低的问题。
【似然比检验】：表格显示最终模型的AIC、BIC、-2LL值(与【模型拟合信息】表的结果一致)，以及简约模型(去掉某个自变量效应后的模型)的AIC、BIC、-2LL值，卡方检验统计量为简约模型与最终模型的-2LL差值。结果显示社区和性别对模型的贡献均有统计学意义。
【参数估计】SPSS中因变量默认以取值高水平为参照水平(本例为社区宣传)，如希望将其他取值水平作为参照水平，可在数据中修改因变量各水平的赋值，或者通过[参考类别(Reference Category…)]来指定。自变量也默认取值水平高的为参照水平，也可以修改自变量各水平的赋值来改变参照水平，如果将变量作为协变量纳入分析则低水平会被默认为参照水平。因此本例中社区B(社区=1)和女(性别=1)为参照水平，其参数值为0，一般是研究者不感兴趣的参数，即冗余参数。
从结果来看，社区A(社区=0)的回归系数为负值，P=0.001<0.05，OR=0.370。具有统计学意义表明社区A回归系数不为0(社区B的回归系数为0)。回归系数为负，表明与社区宣传相比，社区A(比社区B)更不愿意通过传统大众传媒获得健康知识，或者说社区A更愿意通过社区宣传获得健康知识；OR=0.370，即相比社区宣传，社区A通过传统大众传媒获得健康知识是社区B的0.37倍，或者更符合表达逻辑的说法是社区A通过社区宣传获得健康知识是社区B的2.70倍(1/0.370)，社区B通过传统大众传媒获得健康知识是社区A的2.70倍。当然严格来说OR的表达应该是：社区B选择传统大众传媒与选择社区宣传的比值 是 社区A相应比值的2.70倍。
同理可知与社区宣传相比，男性(比女性)更愿意通过传统大众传媒获得健康知识，OR=3.410。与网络宣传相比，社区A(与社区B)在通过传统大众媒体获得健康知识方面没有统计过学差异(Wald χ2=1.7,P=0.192>0.05)，但男性更倾向于选择网络获得健康知识(Wald χ2=8.126,P=0.004<0.05，OR=2.213)。
如想比较传统大众媒介与网络，可直接将相应的模型方程相减，
大体可以判断与网络途径相比，社区A更不倾向传统大众传媒(也就是说更倾向于网络)，男性更倾向传统大众传媒，但是否有统计学意义尚需进一步检验。可在多变量回归对话框中，通过[参考类别(Reference Category…)]将参考类别自定义为网络(Custom Value=2)，可获得如下结果，同上述计算结果一致，解读略。
另外，自变量为多分类时也要遵循同进同出的原则。
【分类表】:观测频数与预测频数的差异。对角线为正确判断的频数，而非对角线为判断错误的个数。预测正确率一般，有待改善。
【观测频数和预测频数】:较为接近，拟合良好。
—— END ——
	

1，引言
logistic回归是机器学习中最常用最经典的分类方法之一，有人称之为逻辑回归或者逻辑斯蒂回归。虽然他称为回归模型，但是却处理的是分类问题，这主要是因为它的本质是一个线性模型加上一个映射函数Sigmoid，将线性模型得到的连续结果映射到离散型上。它常用于二分类问题，在多分类问题的推广叫softmax。

　　本文首先阐述Logistic回归的定义，然后介绍一些最优化算法，其中包括基本的梯度上升法和一个改进的随机梯度上升法，这些最优化算法将用于分类器的训练，最好本文将给出一个Logistic回归的实例，预测一匹病马是否能被治愈。

　　在我们的日常生活中遇到过很多最优化问题，比如如何在最短时间内从A点到达B点？如何投入最少工作量却获得最大的效益？如何设计发动机使得油耗最少而功率最大？可见，最优化的作用十分强大，所以此处我们介绍几个最优化算法，并利用它们训练出一个非线性函数用于分类。

　　现在假设有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称作回归。利用logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类，这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法，下面我们首先介绍一下这个二值型输出分类器的数学原理。

那么逻辑回归与线性回归是什么关系呢？

逻辑回归（Logistic Regression）与线性回归（Linear Regression）都是一种广义线性模型（generalized linear model）。逻辑回归假设因变量 y 服从伯努利分布，而线性回归假设因变量 y 服从高斯分布。 因此与线性回归有很多相同之处，去除Sigmoid映射函数的话，逻辑回归算法就是一个线性回归。可以说，逻辑回归是以线性回归为理论支持的，但是逻辑回归通过Sigmoid函数引入了非线性因素，因此可以轻松处理0/1分类问题。
2，Logistic回归的一般过程


（1）收集数据：采用任意方法收集数据


（2）准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳


（3）分析数据：采用任意方法对数据进行分析


（4）训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数


（5）使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定他们属于哪个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。


3，Logistic回归的优缺点
**优点：**计算代码不多，易于理解和实现，计算代价不高，速度快，存储资源低

　　**缺点：**容易欠拟合，分类精度可能不高

　　**适用数据类型：**数值型和标称型数据
4，基于Logistic回归和Sigmoid函数的分类
**我们想要的函数应该是：能接受所有的输入，然后预测出类型。**例如，在两个类的情况下，上述函数输出0或1。该函数称为海维赛德阶跃函数（Heaviside step function），或者直接称为单位阶跃函数。然而，海维赛德阶跃函数的问题在于：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质（可以输出0或者1），且数学上更易处理，这就是Sigmoid函数。Sigmoid函数具体的计算公式如下：

　　自变量取值为任意实数，值域[0, 1]

　　图5-1给出了Sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时，Sigmoid函数值为0.5。随着x的增大，对应的Sigmoid值将逼近于1；而随着x的减少，Sigmoid值将逼近于0.如果横坐标刻度足够大，Sigmoid函数看起来很像一个阶跃函数。


**解释Sigmoid函数：**将任意的输入映射到了 [0, 1]区间，我们在线性回归中可以得到一个预测值，再将该值映射到 Sigmoid函数中这样就完成了由值到概率的转换，也就是分类任务。

　　因此，为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和带入Sigmoid函数中，进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分入1类，小于0.5即被归入0类，所以，Logistic回归也可以被看成是一种概率估计。

　　确定了分类器的函数形式之后，现在的问题变成了：最佳回归系数是多少？如何确定其大小。
5，基于最优化方法的最佳回归系数确定
Sigmoid函数的输入记为z，由下面公式得到：


如果采用向量的写法，上述公式可以写成  z = wTx  ，它表示将这两个数值向量对应元素相乘，然后全部加起来即得到z值。



　　其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数（系数），从而使得分类器尽可能的准确，为了寻找该最佳参数，需要用到最优化理论的一些知识。

　　然后再看看我们的Logistic回归模型的公式：



　　这里假设 W>0，Y与X各维度叠加的图形关系，如下图所示（x为了方便取1维）：



　　下面首先学习梯度上升的最优化方法，我们将学习到如何使用该方法求得数据集的最佳参数，接下来，展示如何绘制梯度上升法产生的决策边界图，该图将梯度上升法的分类效果可视化的呈现出来，最后我们将学习随机梯度上升算法，以及如何对其进行修改以获得很好地结果。

可能我们最常听到的是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法，梯度上升算法用来求函数的最大值，而梯度下降算法是用来求函数的最小值

6，梯度上升法
梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻，如果梯度记为，则函数 f(x,y) 的梯度由下面式子表示：


这个梯度意味着要沿着x的方向移动

，沿着y方向移动

，其中函数f(x,y)必须要在待计算的点上有定义并且可微，一个具体的函数例子见图5-2：



　　上图中的梯度上升算法沿梯度方向移动了一步，可以看出，梯度算子总是指向函数值增长最快的方向。这里所说的移动方向，而未提到移动量的大小。该量值称为步长，记为。用向量来表示的话，梯度算法的迭代公式如下：



　　该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

　　基于上面的内容，我们来看一个Logistic回归分类器的应用例子，从图5-3可以看到我们采用的数据集。



梯度上升法的公式推导（LR 损失函数）

　　在LR中，应用极大似然估计法估计模型参数，由于Sigmoid函数的特性，我们可以做如下的假设：


上式即为在已知样本X和参数θ的情况下。样本X属性正类（y=1）和负类（y=0）的条件概率，将两个公式合并成一个，如下：



　　假定样本与样本之间相互独立，那么整个样本集生成的概率即为所有样本生成概率的乘积（也就是n个独立样本出现的似然函数如下）：



　　为了简化问题，我们对整个表达式求对数（即为LR 损失函数）：



　　满足似然函数（θ）的最大的θ值即时我们需要求解的模型。

　　那么梯度上升法就像爬坡一样，一点一点逼近极值，而上升这个动作用数学公式表达即为：



　　其中，α 为步长。
回到Logistic回归问题，我们同样对函数求偏导。

对这个公式进行分解，先看：



　　我们可以看到，对函数求偏导，分解为三部分，然后我们对这三部分分布求导。
其中：



再由：


可得：


接下来：

最后：



综合三部分即得到：



　　如果上面链式分解不好理解的话，可以看下面直接求导（结果是一样的）：



　　注意上面是将梯度上升求最大值，转换为梯度下降了，本质没变。

因此梯度迭代公式为：



如果为梯度下降，我们注意符号的变化，如下：


7，训练算法：使用梯度上升找到最佳参数
上图有100个样本点，每个点包含两个数值型特征：X1和X2，在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。
所以我们的目标：建立分类器，求解出theta参数

　　设定阈值，根据阈值判断结果
8 Python中的sklearn.linear_model.LogisticRegression
sklearn.linear_model.LogisticRegression官方API：http://scikitlearn.org/stable/modules/generated/sklearn.linear_model.LogisticRegression.html
class sklearn.linear_model.LogisticRegression(penalty='l2', dual=False, tol=0.0001, C=1.0,fit_intercept=True, intercept_scaling=1, class_weight=None, random_state=None,solver='liblinear', max_iter=100, multi_class='ovr', verbose=0,warm_start=False, n_jobs=1)

ogistic回归（aka logit，MaxEnt）分类器。
在多类情况下，如果将“ multi_class”选项设置为“ ovr”，则训练算法将使用“一对多休息”（OvR）方案；如果将“ multi_class”选项设置为“多项式”，则使用交叉熵损失’。（当前，只有“ lbfgs”，“ sag”，“ saga”和“ newton-cg”求解器支持“多项式”选项。）
该类使用“ liblinear”库，“ newton-cg”，“ sag”，“ saga”和“ lbfgs”求解器实现正则逻辑回归。请注意，默认情况下将应用正则化。它可以处理密集和稀疏输入。使用C排序的数组或包含64位浮点数的CSR矩阵可获得最佳性能；其他任何输入格式将被转换（并复制）。
“ newton-cg”，“ sag”和“ lbfgs”求解器仅支持带有原始公式的L2正则化，不支持正则化。'liblinear’求解器支持L1和L2正则化，仅针对L2罚分采用对偶公式。仅“ saga”求解器支持Elastic-Net正则化。





笔记
底层的C实现在拟合模型时使用随机数生成器选择特征。因此，对于相同的输入数据具有略微不同的结果并不罕见。如果发生这种情况，请尝试使用较小的tol参数。
在某些情况下，预测输出可能与独立liblinear的输出不匹配。请参见 叙述文档中与liblinear的区别。
参考文献
L-BFGS-B –大规模约束优化软件

朱次有，理查德·伯德，豪尔赫·诺德达尔和何塞·路易斯·莫拉莱斯。 http://users.iems.northwestern.edu/~nocedal/lbfgsb.html
LIBLINEAR –大型线性分类的库

https://www.csie.ntu.edu.tw/~cjlin/liblinear/
凹陷-马克·施密特（Mark Schmidt），尼古拉斯·勒·鲁（Nicolas Le Roux）和弗朗西斯·巴赫（Francis Bach）

用随机平均梯度最小化有限求和 https://hal.inria.fr/hal-00860051/document
SAGA – Defazio，A.，Bach F.和Lacoste-Julien S.（2014）。

SAGA：一种支持非强凸复合物镜的快速增量梯度方法 https://arxiv.org/abs/1407.0202
于祥富，黄芳兰，林志仁（2011）。双坐标下降

逻辑回归和最大熵模型的方法。机器学习85（1-2）：41-75。 https://www.csie.ntu.edu.tw/~cjlin/papers/maxent_dual.pdf
例子
>>> from sklearn.datasets import load_iris
>>> from sklearn.linear_model import LogisticRegression
>>> X, y = load_iris(return_X_y=True)
>>> clf = LogisticRegression(random_state=0).fit(X, y)
>>> clf.predict(X[:2, :])
array([0, 0])
>>> clf.predict_proba(X[:2, :])
array([[9.8...e-01, 1.8...e-02, 1.4...e-08],
       [9.7...e-01, 2.8...e-02, ...e-08]])
>>> clf.score(X, y)
0.97...



使用实例sklearn.linear_model.LogisticRegression

正则化逻辑回归
1、什么是正则化？？？


尽可能保留所有特征值
减小除零次项以外所有参数，假设参数索引从0开始，那么即减少$\theta _j(j \geq 1)$的大小
通过增加损失函数中一个惩罚，正则化损失函数公式如下：


$J\left( \theta  \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}$
其中$\lambda$称为正则化参数，$\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}$ 则为其代价函数所增加的惩罚值，以期通过此值来减小$\theta$，从而避免过拟合。






2、如何进行正则化逻辑回归？？？
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
plt.style.use('fivethirtyeight') #样式美化
from sklearn.metrics import classification_report#这个包是评价报告

'''读入数据'''
df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
df.head()
#这里的test1和test2已经做好归一化
sns.set(context="notebook", style="ticks", font_scale=1.5)
sns.lmplot('test1', 'test2', hue='accepted', data=df, 
           size=6, 
           fit_reg=False, 
           scatter_kws={"s": 50}
          )

plt.title('Regularized Logistic Regression')
plt.show()




根据上图可以发现，这样的分类边界不是一个线性方程能解释的，所以需要进行特征映射，将低纬度特征映射成高维的。




2.1 特征映射
polynomial expansion（多项式扩展）
for i in 0…i

for p in 0…i:

output x^(i-p) * y^p

假设扩展的维度为6，即将一个2维特征扩展成28维的多项式特征


def feature_mapping(x, y, power, as_ndarray=False):
    '''
    将特征映射生成一个矩阵或者DataFrame
    前面相当于命名列名，后面相当于将一列的值填充进去
    用的是列表生成式
    '''
    data = {"f{}{}".format(i - p, p): np.power(x, i - p) * np.power(y, p)
                for i in np.arange(power + 1)
                for p in np.arange(i + 1)
            }
    if as_ndarray:
        return pd.DataFrame(data).as_matrix()
    else:
        return pd.DataFrame(data)

'''分别取出两个特征进行多项式映射'''
x1 = np.array(df.test1)
x2 = np.array(df.test2)
data = feature_mapping(x1, x2, power=6)
print(data.shape)
data.head()





2.2正则化损失函数
$J\left( \theta  \right)=\frac{1}{m}\sum\limits_{i=1}^{m}{[-{{y}^{(i)}}\log \left( {{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)-\left( 1-{{y}^{(i)}} \right)\log \left( 1-{{h}_{\theta }}\left( {{x}^{(i)}} \right) \right)]}+\frac{\lambda }{2m}\sum\limits_{j=1}^{n}{\theta _{j}^{2}}$
theta=np.zeros(X.shape[1])
def regular_cost(theta, X, y,l=1):
    ''' cost + lambda/2m*sum(theta)'''
    #这里需要补充矩阵乘法和dataframe中apply的用法，可以写在csdn里面
    costfir = np.mean(-y*np.log(sigmoid(X@theta))-(1-y)*np.log(1-sigmoid(X@theta)))
    theta_2_to_n = theta[1:]
    regular_cost = (l/(2*len(X)))*(np.sum(np.power(theta_2_to_n,2))) 
    return costfir + regular_cost



2.3 正则化梯度下降函数
$\frac{\partial J\left( \theta  \right)}{\partial {{\theta }_{j}}}=\left( \frac{1}{m}\sum\limits_{i=1}^{m}{\left( {{h}_{\theta }}\left( {{x}^{\left( i \right)}} \right)-{{y}^{\left( i \right)}} \right)} \right)x_{_{j}}^{(i)}+\frac{\lambda }{m}{{\theta }_{j}}\text{ }\text{             for  j}\ge \text{1}$
'''构造正则化梯度下降函数，lambda为默认参数'''
def regularized_gradient(theta, X, y, l=0.01):
    #theta0不进行正则化处理，因此从第二项开始
    theta_2_to_n = theta[1:]
    regularized_term = (l/len(X))*theta_2_to_n
    regularized_term_all = np.insert(regularized_term,0,0)
    #regularized_term_all = np.concatenate([np.array([0]),regularized_term]) #这段代码也行，是数组之间的粘贴
    return gradient(theta, X, y) + regularized_term_all
gradient(theta, X, y)
regularized_gradient(theta, X, y)
#上面两个应该一样，因为现在的theta是0


'''拟合参数'''
import scipy.optimize as opt
print('init cost = {}'.format(regularized_cost(theta, X, y)))
res = opt.minimize(fun=regularized_cost, x0=theta, args=(X, y), method='Newton-CG', jac=regularized_gradient)
res
#res里面x表示的就是所有theta的分量
#批量梯度下降
'''可以看下cost的下降'''
def lr_reg_gradientDescent(X,y,theta,alpha,iters):
    theta1 = theta.copy()
    theta_cn = theta1.shape[0]  #不用len()是因为theta2是矩阵，返回的是1
    N = len(X)  #样本容量
    cost2 = np.zeros(iters)  #用来存每次迭代后计算的代价函数
    temp = np.mat(np.zeros(theta1.shape))  #用来传临时theta2，后面做递归用

    for i in range(iters):
        lr_re_gr = regularized_gradient(theta1, X, y)
        # 把temp的两个值其实也就是theta0、theta1传给theta，再计算第i次迭代出来的代价函数的值    
        theta1 = theta1 - alpha*lr_re_gr
        cost2[i] = regularized_cost(theta1, X, y)
    return theta1, cost2

alpha = 0.01
iters = 10000
# 输出theta和cost
lr_theta, lr_cost = lr_reg_gradientDescent(X, y,theta, alpha, iters)
#可视化一下
fig, ax = plt.subplots(figsize=(10,8))
ax.plot(np.arange(iters), lr_cost, 'r')#第一个参数是x，第二个参数是y，最后的参数是线条的颜色
ax.set_xlabel('Iterations')
ax.set_ylabel('Cost')
ax.set_title('Error vs. Training Epoch')
plt.show()


'''预测结果报告'''
final_theta = res.x
y_pred = predict(X, final_theta)

print(classification_report(y, y_pred))

       precision    recall  f1-score   support

       0       0.88      0.75      0.81        60
       1       0.78      0.90      0.83        58

accuracy                           0.82       118

'''根据不同lambda画出决策边界'''
def draw_boundary(power, l):
#     """
#     power: polynomial power for mapped feature
#     l: lambda constant
#     """
    density = 1000
    threshhold = 2 * 10**-3

    final_theta = feature_mapped_logistic_regression(power, l)
    x, y = find_decision_boundary(density, power, final_theta, threshhold)

    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
    sns.lmplot('test1', 'test2', hue='accepted', data=df, size=6, fit_reg=False, scatter_kws={"s": 100})

    plt.scatter(x, y, c='R', s=10)
    plt.title('Decision boundary')
    plt.show()
    
'''特征映射'''
def feature_mapped_logistic_regression(power, l):
     """画图需要，并不是广义逻辑回归
     power: int
         raise x1, x2 to polynomial power
     l: int
        lambda constant for regularization term
    """
    df = pd.read_csv('ex2data2.txt', names=['test1', 'test2', 'accepted'])
    x1 = np.array(df.test1)
    x2 = np.array(df.test2)
    y = get_y(df)

    X = feature_mapping(x1, x2, power, as_ndarray=True)
    theta = np.zeros(X.shape[1])

    res = opt.minimize(fun=regularized_cost,
                       x0=theta,
                       args=(X, y, l),
                       method='TNC',
                       jac=regularized_gradient)
    final_theta = res.x
    return final_theta

'''寻找决策边界，密集打点'''
def find_decision_boundary(density, power, theta, threshhold):
    t1 = np.linspace(-1, 1.5, density)
    t2 = np.linspace(-1, 1.5, density)

    cordinates = [(x, y) for x in t1 for y in t2]
    x_cord, y_cord = zip(*cordinates)
    mapped_cord = feature_mapping(x_cord, y_cord, power)  # this is a dataframe

    inner_product = mapped_cord.as_matrix() @ theta
    decision = mapped_cord[np.abs(inner_product) < threshhold]
    return decision.f10, decision.f01
 
draw_boundary(power=6, l=0.7)     #set lambda = 1

1. 引言
此前的博客中，我们已经介绍了几个分类算法。

k 近邻算法

决策树的构建算法 – ID3 与 C4.5 算法

朴素贝叶斯算法的推导与实践
本文介绍的是另一个分类算法 – 逻辑斯蒂回归。

他凭借易于实现与优异的性能，拥有着十分广泛的使用，他不仅可以进行二分类问题的解决，也可以解决多分类问题，简单起见，本文我们只讨论二分类问题。

它的基本思想是利用一条直线将平面上的点分为两个部分，即两个类别，要解决的问题就是如何拟合出这条直线，这个拟合的过程就称之为“回归”。


2. 逻辑斯蒂算法推导
2.1. 几率函数与 logit 函数
假设一个事件发生的概率是 p，不发生的概率就是 1-p，那么 p/(1-p) 被称为事件发生的“几率”。

f§ = p/(1-p) 就是几率函数，举个简单的例子，足球赛上，A队对抗B队，胜率是 90%，那么通过几率函数可以求得 f(0.9) = 9，也就是说，在10场比赛中，A队可以平均获胜9场，在实际的生活中，几率函数比概率更为实用。

对几率函数取对数就是 logit 函数：


2.2. logistic 函数
我们将上面的例子看作在条件 X = 10 场比赛下，A队获胜概率为 90% 时得到的 A 队获胜 9 场的特征结果。

但我们需要解决的是符合一定特征的样本其属于某个类别的概率，这就需要队上面的函数求反，从而得到 logistic 函数，也称为 Sigmoid 函数：

2.3. 绘制 Sigmoid 函数图
通过下面的 python 代码我们可以绘制出 Sigmoid 函数的函数图：
import matplotlib.pyplot as plt
import numpy as np
 
 
def sigmoid(z):
    return 1.0/(1.0 + np.exp(-z))
 
 
if __name__ == '__main__':
    z=np.arange(-6, 6, 0.05)
    plt.plot(z, sigmoid(z))
    plt.axvline(0.0, color='k')
    plt.axhline(y=0.0, ls='dotted', color='k')
    plt.axhline(y=1.0, ls='dotted', color='k')
    plt.axhline(y=0.5, ls='dotted', color='k')
    plt.yticks([0.0, 0.5, 1.0])
    plt.ylim(-0.1, 1.1)
    plt.xlabel('z')
    plt.ylabel('$\phi (z)$')
    plt.show()


3. 代价函数
假设我们的样本 x 具有 x1、x2、x3 。。。xn 等多个特征，通过添加系数 θ 可以得到下面的公式：

从而得到：

因为我们要解决二分类问题，所以我们只关心发生和不发生的概率，也就是 y=0 与 y=1 的解：

因为 y 只有 0 和 1 两个取值，我们可以通过下面函数将上面两个函数整合在一起：

这个函数就是代价函数。
4. 极大似然估计
有了代价函数，我们只要找到所有的 θ 使得对于我们的所有样本都能成立就可以了，这个找 θ 的过程就是极大似然估计。

极大似然估计的目的就是利用已知的样本结果，反推最有可能（最大概率）导致这样结果的参数值，从而将概率密度估计问题转化为参数估计问题。

但他有一个明显的前提需要保证：训练样本的分布能代表样本的真实分布。

p(x1, x2, x3, …, xn|θ) 被称为 {x1, x2, x3, …, xn} 的 θ 的似然函数。

在假设全部特征 x1, x2, x3 …, xn 均为独立特征的前提下，我们可以得到：

5. 梯度上升法求解极大似然函数
5.1. 准备求解
让 l(θ) 函数获得最大值的 θ 值就是我们想要的值。

正如上一篇日志中所写的，连续积运算通过对数转换为累加是数学中的常用方法，因此我们对上面函数两边求 log。

假设函数连续且任一点可导，那么，要计算一个函数结果取到极值时的参数值，我们通常只需要让函数的导数为 0 即可计算出函数取极值时的参数值。
我们定义梯度因子为：

那么我们希望求解：

这个函数的直接求解太过复杂，因此不能直接求解，有几种方法可以通过迭代的方式得到近似的估计值，例如牛顿-拉菲森迭代方法。

本文我们介绍另一种迭代方法 – 梯度上升法，梯度上升算法和牛顿迭代相比更加易于理解，但收敛速度慢，因为梯度上升算法是一阶收敛，而牛顿迭代属于二阶收敛。
5.2. 梯度上升法
梯度上升法的原理是像爬坡一样，一点一点的向上爬，等到爬到山顶，那就已经非常接近我们的极值点了。

也就是说，每次比较上一次的函数值和本次的函数值，如果本次的函数值大于上一次的函数值，说明正处于上升阶段，否则说明已经进入了下降阶段，那么上次的点就是极值点。

但是这里有一个上升步长的问题，如果我们步子迈得太大，那么很容易就越过山顶了，最后只把半山腰上的某个点认成了极值点，而如果步子迈得太小，那么就会耗费更加大量的时间去不断地爬坡。

根据我们的经验，函数的导数是随着我们接近极值点而逐渐减小的，所以我们只要让步长正比于函数的导数，就可以让我们的步子随着极值点的接近越来越小，因此我们定义上升函数：

α 是步长系数。
5.3. 求解 l(θ) 的导数
现在，我们只要求出 l(θ) 导数带入上面的公式就可以利用梯度上升算法，求解 l(θ) 的极大值了。

接下来就是数学推导过程了：


通过添加项的方法，我们可以将求导问题归类为三部分的求解：

第一部分



第二部分

根据：


可得：


第三部分


综上：

从而得到梯度上升迭代公式：

6. 参考资料
https://en.wikipedia.org/wiki/Maximum_likelihood_estimation。

https://en.wikipedia.org/wiki/Gradient_descent。

https://zlatankr.github.io/posts/2017/03/06/mle-gradient-descent。

https://blog.csdn.net/c406495762/article/details/77723333。

https://blog.csdn.net/zengxiantao1994/article/details/72787849。

https://blog.csdn.net/sinat_29957455/article/details/78944939。

https://blog.csdn.net/amds123/article/details/70243497。

https://blog.csdn.net/gwplovekimi/article/details/80288964。

https://www.cnblogs.com/HolyShine/p/6395965.html。