三、【机器学习作业】逻辑（Logistic）回归算法 & 正则化（python版ex2）

学习完了机器学习的逻辑回归课程，现在来将所学记录下来。

（1）Logistic回归算法模型

更新
简而言之，逻辑回归是解决二分类问题的，简单来说就是 “线性回归+sigmoid函数”。

---

概念：
逻辑回归（Logistic Regression）是用于处理因变量为分类变量的回归问题，常见的是二分类或二项分布问题，也可以处理多分类问题，它实际上是属于一种分类方法。比如常见的分类是：把邮件分成垃圾邮件和非垃圾邮件；把肿瘤分成良性和恶性，等等。对于类别我们通常称为正类（positive class） 和负类（negative class），垃圾邮件的例子中，正类就是正常邮件，负类就是垃圾邮件。

逻辑回归算法产生的原因：
我们知道，线性回归的模型是求出输出特征向量Y和输入样本矩阵X之间的线性关系系数θ，满足Y=Xθ。此时我们的Y是连续的，所以是回归模型。如果我们想要Y是离散的话，怎么办呢？一个可以想到的办法是，我们对于这个Y再做一次函数转换，变为g(Y)。如果我们令g(Y)的值在某个实数区间的时候是类别A，在另一个实数区间的时候是类别B，以此类推，就得到了一个分类模型。如果结果的类别只有两种，那么就是一个二元分类模型了。逻辑回归的出发点就是从这来的。

“分类问题”和“回归问题”的差别：
其实分类和回归的本质是一样的，都是对输入做出预测，其区别在于输出的类型。
分类问题：分类问题的输出是离散型变量(如: +1、-1)，是一种定性输出。(预测明天天气是阴、晴还是雨)
回归问题：回归问题的输出是连续型变量，是一种定量输出。(预测明天的温度是多少度)。

“分类问题”和 “回归问题”之间的转化：
在某些情况下，可以将回归问题转换为分类问题。例如，要预测的数量可以转换为离散桶。
例如，$ 0到$ 100之间连续范围内的金额可以转换为2个桶：

0级：0美元到49美元
1级：50美元到100美元

这通常称为离散化，结果输出变量是标签具有有序关系（称为序数）的分类。
在某些情况下，分类问题可以转换为回归问题。例如，标签可以转换为连续范围。
可以对类值进行排序并映射到连续范围：

1级 $ 0到$ 49
2级 $ 50至$ 100

特点：

• 逻辑回归的数学模型和求解都相对比较简洁，实现相对简单。
• 逻辑回归适用于二分类问题，或者多分类问题。
• 逻辑回归可以处理线性问题，也可以处理非线性问题。

优点：

• 可以适用于连续性和类别性自变量。
• 预测结果是介于0和1之间的概率。
• 容易使用和解释。

缺点：

• 对模型中自变量多重共线性较为敏感。
• 对数据和场景的适应能力有局限性。
• 容易欠拟合，分类精度不高。

逻辑回归原理总结可以参考 https://www.cnblogs.com/pinard/p/6029432.html

（2）公式 原理推导

Logistic回归算法也分成三个步骤：

1. 需要确定一个预测函数，即预测出一个值来判断归属哪一类，可定义预测值大于某个阈值判断为一类，反之为另一类；
2. 为了计算参数，我们需要定义一个代价函数，代价函数用来衡量真实值和预测值之间的差异，并且该差值越小，则说明预测效果越好；
3. 使用梯度下降法计算参数，使代价函数不断减小，计算出回归参数。

【1】 预测函数公式：

由于逻辑回归算法和线性回归算法的不同，逻辑回归算法的输出值是0或1，那么需要产生0或1的值就需要导入sigmoid函数，该函数的公式为：$$g(z)=\frac{1}{1+e^{-z}}$$ 该函数的图像为：

它有一个非常好的性质，即当z趋于正无穷时，g(z)趋于1，而当z趋于负无穷时，g(z)趋于0，这非常适合于我们的分类概率模型。

这样就符合逻辑回归算法所需要的输出值为$[0,1]$。如果令g(z)中的z为：$z={\theta }^{T}x$，那么就可以将逻辑回归算法的假设函数改写为：$$h_{\theta }(x)=g({\theta }^{T}x)=\frac{1}{1+e^{-{\theta }^{T}x}}$$
其中$x$为样本输入，$h_{\theta }(x)$为模型输出，可以理解为某一分类的概率大小。而$\theta$为分类模型的要求出的模型参数。

$h_{\theta }(x)$的值越小，而分类为0的的概率越高，反之，值越大的话分类为1的的概率越高。如果靠近临界点，则分类准确率会下降。

由上面的函数图像可以看出， sigmoid 函数可以很好地将 $(-\infty ,+\infty )$ 内的数映射到 $[0,1]$ 上。可以将 $g(z)\ge 0.5$ 时分为"正"类， g(z)<0.5 时分为"负"类，这就可以判断边界了。需要注意的是，这里选择0.5作为阈值，只是一般的做法，在实际应用中，需要根据实际情况选择不同的阈值。

【2】 代价函数公式：

回顾下线性回归的代价函数，由于线性回归是连续的，所以可以使用模型误差的的平方和来定义代价函数。但是逻辑回归不是连续的，线性回归代价函数定义的经验就用不上了。不过我们可以用最大似然法来推导出代价函数。

根据二元逻辑回归的定义，假设样本输出是0或者1两类。那么使用公式来表达就是：$$P(y=1|x;\theta )=h_{\theta }(x)——正类$$ $$P(y=0|x;\theta )=1-h_{\theta }(x)——负类$$ 合并起来：$$P(y|x;\theta )=(h_{\theta }(x){)}^y(1-h_{\theta }(x){)}^{1-y}$$ 对于二分类的任务，$y$的取值只能是0或者1。当y=0时，只留下$(1-h_{\theta }(x){)}^{1-y}$；当y=1时，只留下$(h_{\theta }(x){)}^y$。

得到了$y$的概率分布函数表达式，就可以使用似然函数最大化来求解所需要的模型系数$\theta$了。

由上面的式子可以得到二分类的最大似然函数：$$L(\theta )=\prod _{i=1}^{m}P(y^{(i)}|x^{(i)};\theta )=\prod _{i=1}^m(h_{\theta }(x^{(i)}){)}^{y^{(i)}}(1-h_{\theta }(x^{(i)}){)}^{1-y^{(i)}}$$其中，$m$为样本的个数。

对数的最大似然函数为：$$logL(\theta )=\sum _{i=1}^m[y^{(i)}\cdot log(h_{\theta }(x^{(i)}))+(1-y^{(i)})\cdot log(1-h_{\theta }(x^{(i)}))]$$

根据上面的推导，便可以得到了代价函数，对 $logL(\theta )$ 求它的最大值来求得参数 $\theta$ 的值。由于上面的公式是求似然函数的最大值，即求梯度上升的值，那么求梯度下降的最小值就需要添加“-”号。

对似然函数对数化取反的表达式，即代价函数表达式（未正则化） 为：

$$\begin{gathered} J(\theta )=\frac{1}{m}\sum _{i=1}^{m}Cost(h_{\theta }(x),y) \\ =-\frac{1}{m}[\sum _{i=1}^m{y}^{(i)}\cdot log(h_{\theta }(x^{(i)}))+(1-y^{(i)})\cdot log(1-h_{\theta }(x^{(i)}))] \end{gathered}$$
其中，$$\begin{gathered} Cost(h_{\theta }(x),y)=-[y^{(i)}\cdot log(h_{\theta }(x^{(i)}))+(1-y^{(i)})\cdot log(1-h_{\theta }(x^{(i)}))] \\ =\{\begin{array}{ccc}-log(h_{\theta }(x^{(i)}))& if& y=1\\ -log(1-h_{\theta }(x^{(i)}))& if& y=0\end{array} \end{gathered}$$
根据上面的公式，对$J(\theta )$求最小值，即可得到参数$\theta$。

向量正则化的逻辑回归代价函数：
$$J(\theta )=-\frac{1}{m}[\sum _{i=1}^m{y}^{(i)}\cdot log(h_{\theta }(x^{(i)}))+(1-y^{(i)})\cdot log(1-h_{\theta }(x^{(i)}))]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$

【3】 梯度下降参数更新公式：

对于二元逻辑回归的代价函数最小化，有比较多的方法，最常见的有梯度下降法，坐标轴下降法，等牛顿法等。这里推导出梯度下降法中$\theta$每次迭代的公式。

想要$J(\theta )$求得最小值，并求得参数$\theta$，通常使用的是梯度下降法。梯度下降公式与线性回归算法中运用的公式类似，不同的是线性回归算法和逻辑回归算法中的 假设函数$h_{\theta }(x)$是不一样的，也可以说这是两种形式一样，而数值计算方法不一样的梯度下降算法。特别说明一点，特征缩放也可以应用在逻辑回归算法中，帮助提高梯度下降的收敛速度。

梯度下降参数更新公式为：$${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}$$ 其中，$\alpha$是学习率（也称为“步长”）， $i=1,2,...,m$，表示样本数；$j=1,2,..,n$，表示特征数。$(h_{\theta }(x^{(i)})-y^{(i)})$表示误差，即预测值减去真实值。需要注意的是，同时更新${\theta }_j$的值。

梯度下降算法的矩阵表达式：
$$\begin{array}{r}\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right]:=\left[\begin{array}{c}{\theta }_0\\ {\theta }_1\\ ⋮\\ {\theta }_n\end{array}\right]-\alpha \frac{1}{m}\left[\begin{array}{ccc}{x}_0^{(1)}& \cdots & x_n^{(1)}\\ ⋮& \ddots & ⋮\\ x_0^{(m)}& \cdots & x_n^{(m)}\end{array}\right]E\\ =\theta -\alpha \frac{1}{m}{X}^{T}E\\ =\theta -\alpha \frac{1}{m}{X}^T(h_{\theta }(X)-y)\\ =\theta -\alpha \frac{1}{m}{X}^T(g(X\theta )-y)\end{array}$$

逻辑回归的梯度函数：
当$j=0$时，
$$\frac{\partial J(\theta )}{\partial {\theta }_0}=\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)}$$
当$j\ge 1$时，$$\frac{\partial J(\theta )}{\partial {\theta }_j}=(\frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})x_j^{(i)})+\frac{\lambda }{m}{\theta }_j$$

【4】 优化算法：

判断优化算法优劣的方法是：看其是否收敛；参数是否达到了稳定值；是否还在波动变换；期望得到的参数避免有很大的波动，其能收敛到某个值且收敛速度加快。

除了梯度下降算法之外，还有下面三种算法：

1. Conjugate Gradient（共轭梯度法）
2. BFGS
3. L-BFGS

上面三点的优缺点是：不需要手动选择学习率$\alpha$，比梯度下降算法的运算速度更快；但是算法更加复杂。

（3）多元分类：一对多

针对多元分类问题，y={0,1,2,3,…,n}，总共有$n+1$个类别。其解决思路是：首先把问题转换为二元分类问题，即$y=0$是一个类别，y={1,2,3,…,n}作为另外一个类别，然后计算这两个类别的概率；接着，把$y=1$作为一个类别，把y={0,2,3,…,n}作为另外一个类别，再计算这两个类别的概率。由此类推，总共需要n+1个预测函数。

计算过程为：
$$h_{\theta }^{(0)}(x)=P(y=0|x;\theta )$$$$h_{\theta }^{(1)}(x)=P(y=1|x;\theta )$$$$\cdots$$$$h_{\theta }^{(n)}(x)=P(y=n|x;\theta )$$ $$prediction=ma{x}_i{h}_{\theta }^{(n)}(x)$$

预测出来的概率最高的那个类别就是样本所属的类别。

参考例子如下，要计算出Class1（三角形）的概率就把其他的类别看成一类，使用二分类方法求解；同理，要计算出Class2（正方形）和Class3（红叉）分别将其他两类看成一类，同样使用二分类计算概率。比较三个概率值，最高的概率值对应的类别就是该样本所属的类别。

（4）正则化

逻辑回归也会面临过拟合问题，所以我们也要考虑正则化。常见的有L1正则化和L2正则化。逻辑回归的L1正则化的损失函数表达式如下所述，相比普通的逻辑回归代价函数，增加了L1的范数做作为惩罚，超参数α作为惩罚系数，调节惩罚项的大小。

1. 欠拟合 & 过度拟合：

在线性回归算法中，

在逻辑回归算法中，

欠拟合是指忽略了训练集中某些数据，未能学习训练数据中的关系。
过拟合是指过分依赖训练数据，预测未知的数据的准确性较差。

【通俗来说，就是严格关注数据会过拟合，忽略数据会欠拟合。】

过度拟合的原因是：该问题是来源于过多的特征。

解决方法：
1）减少特征数量（减少特征会失去一些信息，即使特征选的很好）
• 可用人工选择要保留的特征。
• 模型选择算法。
2）正则化（特征较多时比较有效）
• 保留所有特征，但减少$\theta$的大小。

2. 正则化的方法：

正则化是结构风险最小化策略的实现，是在经验风险上加一个正则化项或惩罚项。正则化项一般是模型复杂度的单调递增函数，模型越复杂，正则化项就越大。

3. 线性回归模型的正则化：

正则化后的代价函数变为：$$J(\theta )=\frac{1}{2m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)}{)}^2+\lambda \sum _{j=1}^n{\theta }_j^2$$

该公式前半部分是线性回归模型的代价函数，后半部分是正则化项。加上这个部分有两个目的，第一是维持对训练样本的拟合，第二是避免对训练样本的过度拟合。从数学方面来看，代价函数增加了一个正则化项后，代价函数不再仅仅由预测值和真实值的误差所决定，还和参数$\theta$的大小有关。有了这个限制之后，要实现代价函数最小的目的，参数$\theta$就不能随便取值，通过调节参数$\lambda$，控制正则化项，避免过度拟合的情况。
正则化后的梯度下降变换为：$${\theta }_0:={\theta }_0-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_0^{(i)}$$ $${\theta }_j:={\theta }_j-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}+\frac{\lambda }{m}{\theta }_j$$ 其中，$j=1,2,...,n$
合并后的式子为：$${\theta }_j:={\theta }_j(1-\frac{\alpha \lambda }{m})-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}$$ 其中，$j=1,2,...,n$

用于拟合线性回归模型，这里使用的方法是 “正规方程” 。
（正规方程指的是将方程式数目正好等于未知数的个数，从而可求解出这些未知参数的方程。正规方程的推导 https://blog.csdn.net/mary_0830/article/details/96502269）
假设$X=\left[\begin{array}{c}(x^{(1)}{)}^T\\ ⋮\\ (x^{(m)}{)}^T\end{array}\right]$，$y=\left[\begin{array}{c}{y}^{(1)}\\ ⋮\\ y^{(m)}\end{array}\right]$，其中X中的每一行都代表一个单独的训练样本，y表示训练集里的所有标签。

为了最小化代价函数$J(\theta )$，就需要计算$\theta$的最小值，公式为：（假设以n=3为例）$$\theta =(X^{T}X+\lambda \left[\begin{array}{cccc}0& 0& 0& 0\\ 0& 1& 0& 0\\ 0& 0& 1& 0\\ 0& 0& 0& 1\end{array}\right]{)}^{-1}{X}^{T}y$$ 该矩阵为(n+1)×(n+1)维。通过上面公式的计算，可以计算出代价函数$J(\theta )$的最小值。

4. 逻辑回归模型的正则化：

对逻辑回归模型的代价函数进行正则化，其方法也是在原来的代价函数的基础上加上正则化项：$$J(\theta )=-\frac{1}{m}[\sum _{i=1}^m{y}^{(i)}\cdot log(h_{\theta }(x^{(i)}))+(1-y^{(i)})\cdot log(1-h_{\theta }(x^{(i)}))]+\frac{\lambda }{2m}\sum _{j=1}^n{\theta }_j^2$$
正则化后的梯度下降公式为：$${\theta }_j:={\theta }_j(1-\frac{\alpha \lambda }{m})-\alpha \frac{1}{m}\sum _{i=1}^m(h_{\theta }(x^{(i)})-y^{(i)})\cdot x_j^{(i)}$$ 其中，$j=1,2,...,n$。由于${\theta }_0$没有参与正则化。逻辑回归和线性回归看起来形式是一样的，但其实它们的算法并不一样，因为两个式子的预测函数$h_{\theta }(x)$是不一样的。

5. $\lambda$是正则项系数。

当它的值很大时，说明对模型的复杂度惩罚大，对拟合数据的惩罚小，这样就不会出现过拟合情况。 在训练数据中的偏差比较大，在未知数据上的方差比较小，可能会出现欠拟合的情况。
同样，如果它的值小时，说明比较注重对训练数据的拟合，使得在训练数据上的偏差比较小，在未知数据上的方差较大，从而导致过拟合。

（5）python代码实现——logistic回归算法作业

1. 根据要求，建立一个logistic回归模型来预测一名学生是否能合格进入大学，根据两次考试的结果来决定每个申请人的录取机会。有以前的申请人的历史数据， 可以用它作为逻辑回归的训练集。编写实验代码：

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

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
import matplotlib.ticker as ticker
import test
plt.style.use('fivethirtyeight') #样式美化
alpha = 0.004
iterations = 150000

def sigmoid(z):
    return 1 / (1 + np.exp(-z))

x_plot = np.linspace(-10., 10., 5000)
y_plot = sigmoid(x_plot)
plt.plot(x_plot, y_plot)
plt.title('Sigmoid')
plt.show()

# 代价函数和梯度的求解
def Cost_function(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    return np.sum(first - second) / (len(X))

def predict(theta, X):
    probability = sigmoid(X * theta.T)
    return [1 if x >= 0.5 else 0 for x in probability]

def draw():
    # 定义x y数据 x1 y1:未通过 x2 y2:通过
    x1 = []
    y1 = []
    x2 = []
    y2 = []

    # 导入训练数据
    train_data = open("ex2data1.txt")
    lines = train_data.readlines()
    for line in lines:
        scores = line.split(",")
        # 去除标记后面的换行符
        isQualified = scores[2].replace("\n", "")
        # 根据标记将两次成绩放到对应的数组
        if isQualified == "0":
            x1.append(float(scores[0]))
            y1.append(float(scores[1]))
        else:
            x2.append(float(scores[0]))
            y2.append(float(scores[1]))

    #设置样式参数,默认主题 darkgrid（灰色背景+白网格）,调色板 2色
    sns.set(context="notebook", style="darkgrid", palette=sns.color_palette("RdBu", 2),color_codes=False)
    
    #创建两个分数的散点图
    positive = data[data['Admitted'].isin([1])]
    negative = data[data['Admitted'].isin([0])]
    
    #使用颜色编码来可视化
    fig, ax = plt.subplots(figsize=(12,8))
    #                                                  点的大小  颜色   记号形状  设置标签
    ax.scatter(positive['Exam 1'], positive['Exam 2'], s=30, c='g', marker='o', label='Admitted')
    ax.scatter(negative['Exam 1'], negative['Exam 2'], s=30, c='r', marker='x', label='Not Admitted')
    ax.legend()
    ax.set_xlabel('Exam 1 Score')
    ax.set_ylabel('Exam 2 Score')

    # 设置坐标轴上刻度的精度为一位小数。因训练数据中的分数的小数点太多，若不限制坐标轴上刻度显示的精度，影响最终散点图的美观度
    plt.gca().xaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))
    plt.gca().yaxis.set_major_formatter(ticker.FormatStrFormatter('%.1f'))

    # 设置训练得到的模型对应的直线，即h(x)对应的直线
    # 设置x的取值范围：[30, 110]步长为10
    x = np.arange(30, 110, 10)
    y = (-result.x[0] - result.x[1] * x) / result.x[2]
    plt.plot(x, y)

    # 显示
    plt.show()

path = 'ex2data1.txt'
#读取CSV（逗号分割）文件到 DataFrame
data = pd.read_csv(path, header=None, names=['Exam 1', 'Exam 2', 'Admitted'])

#查看数据集的前 5 行
data.head()
print(data.head())

#查看数据的分布情况
data.describe()
print(data.describe())


#数据预处理
# add ones column
data.insert(0, 'Ones', 1)
# set X (training data) and y (target variable)
cols = data.shape[1]
X = data.iloc[:,0:cols-1]
y = data.iloc[:,cols-1:cols]
# convert to matrices and initialize theta
X = np.array(X.values)
y = np.array(y.values)
theta = np.matrix(np.array([0,0,0]))
#X.shape, theta.shape, y.shape
print('X.shape:',X.shape)
print('theta.shape:',theta.shape)
print('y.shape:',y.shape)
#计算一下初始损失值大小
Cost_function(theta, X, y)
print('Init_Cost_function:',Cost_function(theta, X, y))

# 梯度下降函数
def gradient(theta, X, y):
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    parameters = int(theta.ravel().shape[1])
    grad = np.zeros(parameters)
    #误差
    error = sigmoid(X * theta.T) - y
    for i in range(parameters):
        #点乘 元素相乘
        term = np.multiply(error, X[:,i])
        grad[i] = np.sum(term) / len(X)
    return grad

gradient(theta,X,y)
print('gradient:',gradient(theta,X,y))

#func为优化目标函数，x0为依赖的参数，fprime为梯度下降函数
import scipy.optimize as opt   #行速度通常远远超过梯度下降
result = opt.fmin_tnc(func=Cost_function, x0=theta, fprime=gradient, args=(X, y))
print('result:',result)

#result = opt.minimize(fun=cost, x0=theta, args=(X, y), method='Newton-CG', jac=gradient)
#print(result)

# 设置训练得到的模型对应的直线，即h(x)对应的直线
# 设置x的取值范围：[30, 110]步长为10
#X = np.arange(30, 110, 10)
#y = (-result.X[0] - result.X[1] * Xx) / result.X[2]
#plt.figure()
#plt.plot(X, y)

#代入参数值，计算损失函数的值
Cost_function(result[0], X, y)
print('Cost_function:',Cost_function(result[0], X, y))

#预测准确率评估
theta_min = np.matrix(result[0])
predictions = predict(theta_min, X)
print('predictions:',predictions)

correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]
#map()函数将correct转化为全部为int的列表
accuracy = sum(correct) / len(X)
print ('accuracy = {0}%'.format(accuracy))

# 绘图
draw()

实验结果：
sigmoid 函数图：

显示数据中前5行以及数据的分布情况：

绘制的散点图：

X，y，theta的维数：

计算得到初始的代价函数的值：

批量梯度下降：

梯度下降优化后的代价函数，寻找的参数$\theta$值：

带入参数值计算得到的代价函数的值：

用训练集进行预测：

预测结果的准确率：

寻找到的决策边界：

2. 逻辑回归正则化：
   假设您是工厂的产品经理，在两次不同的测试中获得了一些微芯片的测试结果，要从这两个测试中，想确定应该接受还是拒绝微芯片，就需要建立正则化的逻辑回归模型，编写实验代码：

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

import pandas as pd
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import seaborn
from sklearn.linear_model import LogisticRegression#调用sklearn的线性回归包

#sigmoid函数实现
def sigmoid(z):
    return 1 / (1 + np.exp(-z))

#计算原数据集的预测情况
def predict(theta, X):
    theta = np.matrix(theta)
    X = np.matrix(X)

    probability = sigmoid(X * theta.T)
    return [1 if i > 0.5 else 0 for i in probability]

def plot_data():
    positive = data2[data2['Accepted'].isin([1])]
    negative = data2[data2['Accepted'].isin([0])]

    fig, ax = plt.subplots(figsize=(8,5))
    ax.scatter(positive['Test 1'], positive['Test 2'], s=50, c='b', marker='o', label='Accepted')
    ax.scatter(negative['Test 1'], negative['Test 2'], s=50, c='r', marker='x', label='Rejected')
    ax.legend()
    ax.set_xlabel('Microchip Test 1 Score')
    ax.set_ylabel('Microchip Test 2 Score')
    plt.show()

def feature_mapping(x1, x2, power):#特征映射
    data = {}
    for i in np.arange(power + 1):
        for p in np.arange(i + 1):
            data["f{}{}".format(i - p, p)] = np.power(x1, i - p) * np.power(x2, p)

#     data = {"f{}{}".format(i - p, p): np.power(x1, i - p) * np.power(x2, p)
#                 for i in np.arange(power + 1)
#                 for p in np.arange(i + 1)
#             }
    return pd.DataFrame(data)

def Cost_function(theta, X, y):# 代价函数和梯度的求解
    theta = np.matrix(theta)
    X = np.matrix(X)
    y = np.matrix(y)
    first = np.multiply(-y, np.log(sigmoid(X * theta.T)))
    second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T)))
    return np.sum(first - second) / (len(X))

def costReg(theta, X, y, l=1):
    # 不惩罚第一项
    _theta = theta[1: ]
    reg = (l / (2 * len(X))) *(_theta @ _theta)  # _theta@_theta == inner product
    
    return Cost_function(theta, X, y) + reg

def gradient(theta, X, y):
    return (X.T @ (sigmoid(X @ theta) - y))/len(X)  
# the gradient of the cost is a vector of the same length as θ where the jth element (for j = 0, 1, . . . , n)

def gradientReg(theta, X, y, l=1):
    reg = (1 / len(X)) * theta
    reg[0] = 0  
    return gradient(theta, X, y) + reg

path2 = 'ex2data2.txt'
data2 = pd.read_csv(path2, header=None, names=['Test 1', 'Test 2', 'Accepted'])
#查看数据集的前 5 行
data2.head()
print(data2.head())
#查看数据的分布情况
data2.describe()
print(data2.describe())
#绘制散点图
plot_data()

x1 = data2['Test 1'].as_matrix()
x2 = data2['Test 2'].as_matrix()
_data2 = feature_mapping(x1, x2, power=6)
_data2.head()
# 这里因为做特征映射的时候已经添加了偏置项，所以不用手动添加了。
X = _data2.as_matrix()
y = data2['Accepted'].as_matrix()
theta = np.zeros(X.shape[1])
X.shape, y.shape, theta.shape  # ((118, 28), (118,), (28,))
print('X.shape:',X.shape)
print('theta.shape:',theta.shape)
print('y.shape:',y.shape)

cost_Reg =costReg(theta, X, y, l=1)  #     0.6931471805599454
print('Init_Cost_function:',cost_Reg)

gradient_Reg = gradientReg(theta, X, y, 1)
print('Init_gradient:',gradient_Reg)

result2 = opt.fmin_tnc(func=costReg, x0=theta, fprime=gradientReg, args=(X, y, 2))
print('result:',result2)

#Evaluating logistic regression
final_theta = result2[0]
predictions = predict(final_theta, X)
print('predictions:',predictions)

correct = [1 if a==b else 0 for (a, b) in zip(predictions, y)]
accuracy = sum(correct) / len(X)
print ('accuracy = {0}'.format(accuracy))

#Decision boundary
x = np.linspace(-1, 1.5, 250) #修改这个数值，会改变范围的平滑程度
xx, yy = np.meshgrid(x, x)

z = feature_mapping(xx.ravel(), yy.ravel(), 6).as_matrix()
z = z @ final_theta
z = z.reshape(xx.shape)

plot_data()
plt.contour(xx, yy, z, 0)
plt.ylim(-.8, 1.2)

实验结果：
显示数据中前5行以及数据的分布情况：

绘制的散点图：

观察上图发现，正负两类数据并没有线性的决策界限，因此直接用logistic回归算法在这个数据集上并不能表现良好。因此需要进行正则化。

X，y，theta的维数：

计算得到初始的代价函数$J(\theta )$的值：

初始的梯度下降值：

梯度下降优化后的代价函数，寻找的参数$\theta$值：

带入参数$\theta$值计算得到的代价函数$J(\theta )$的值：

用训练集进行预测：

预测结果的准确率：

寻找到的决策边界：

（LogisticRegression 这个模块需要更新才可以使用，更新的命令是：pip install -U scikit）