逻辑回归算法实现鸢尾花数据集的二分类

1、介绍

前言
为了更好的理解本章内容，请参考下面的网址，自主学习一些矩阵求导运算。

https://zhuanlan.zhihu.com/p/158182135机器学习约定规则：
向量对标量求导，默认使用分子布局（本博客采用了这一条）
矩阵对标量求导，默认使用分子布局
标量对向量求导，默认使用分母布局
向量对向量求导，常使用分母布局

1、概述
逻辑回归（Logistic Regression）虽然是名词是回归，但是实际上它是一种分类模型。假设输入值函数（单值x1x_1x1）如下：
z=WTx1=w0+w1x11+w2x12+⋯+wpx1pW∈Rp+1,x1∈Rp+1z=W^Tx_1=w^0+w^1x_{11}+w^2x_{12}+\cdots+w^px_{1p} \\ W \in R^{p+1} , x_1 \in R^{p+1} z=WTx1=w0+w1x11+w2x12+⋯+wpx1pW∈Rp+1,x1∈Rp+1
逻辑回归模型通过sigmoid函数将输入值映射到[0,1]区间范围。其中，逻辑回归函数形式如下：
p=s(z)=11+e−zp=s(z)=\frac{1}{1+e^{-z}} p=s(z)=1+e−z1
逻辑回归函数图像如下：

s(z)函数的定义域范围为(−∞,+∞)(-\infty,+\infty)(−∞,+∞)，sigmoid函数是一个（0，1）上的函数，从而可以将连续数据进行分类，这可能就是“回归”二字的含义吧。此外，我们也可以将s(z)函数看成概率p，假设p=0.8，此时概率值离1很近，可以将数据类别判别为1；假设p=0.2，概率值离0很近，可以将数据类别判别为0。
因此，假设我们已经训练好了一组权值W ，只要把我们需要预测的数据X代入到方程，输出的y值就是这个标签为0/1的概率，进而判断输入数据是属于哪个类别。下面我们分析使用梯度下降法来分析权值W的求解过程。
2、梯度下降法求解权值W
（1）损失函数的求解过程
假设二分类问题的输出值是一个包含2个离散值函数，即y∈{0,1}y \in \{0,1\}y∈{0,1}，令p等于y=1y=1y=1发生的概率，将每个样本看成一个事件，每个事件发生的概率如下：
P(yi∣xi)={p,y=11−p，y=0P(y_i|x_i)= \begin{cases} p,y=1 \\ 1-p，y=0 \end{cases} P(yi∣xi)={p,y=11−p，y=0
上面的表达式等价于下面的形式：
P(yi∣xi)=pyi(1−p)1−yiP(y_i|x_i)=p^{y_i}(1-p)^{1-y_i} P(yi∣xi)=pyi(1−p)1−yi
假设采集到服从独立同分布的样本N个数据，其联合概率就是每个样本的概率的乘积，此时似然函数有：
L(W)=∏i=1Npyi(1−p)1−yiL(W)=\prod_{i=1}^N{p^{y_i}(1-p)^{1-y_i}} L(W)=i=1∏Npyi(1−p)1−yi
为了计算简便，对上面的似然函数求对数：
F(W)=ln(L(W))=∑i=1N(yiln(p)+(1−yi)ln(1−p))F(W)=ln(L(W))=\sum_{i=1}^N({{y_i}ln(p)+({1-y_i})ln(1-p)}) F(W)=ln(L(W))=i=1∑N(yiln(p)+(1−yi)ln(1−p))
上面的表达式求出最大值，对应的似然函数的求解，最大化似然函数对应最小化代价函数，因此可以在前面加上负号，从而将上面的表达式看成我们的损失函数。
（2）对损失函数求梯度
假设输入值的特征向量有n个，可以对输入值特征向量进行下面的扩展，此时输出值是一个n维的向量，其中向量的每一个元素表示一个输入值的特征向量对应的输出。
[x1T，1x2T，1⋯xnT，1]\left[ \begin{array}{c} x_1^T，1 \\ x_2^T ，1\\ \cdots \\ x_n^T，1 \end{array} \right] ⎣⎢⎢⎡x1T，1x2T，1⋯xnT，1⎦⎥⎥⎤
假设逻辑回归函数形式使用下面的形式:
p′=s(z)′=11+e−XWp′∈Rnp^{'}=s(z)^{'}=\frac{1}{1+e^{-XW}} \\ p^{'} \in R^{n} p′=s(z)′=1+e−XW1p′∈Rn
p′p^{'}p′进行wiw_iwi求导：
(p′)′=(11+e−XW)′=−(11+e−XW)−2.(1+e−XW)′=−(11+e−XW)−2.e−XW.(−XW)′=(11+e−XW)−2.e−XW.Xi=p′(1−p′)Xi(p^{'})^{'}=(\frac{1}{1+e^{-XW}})^{'} \\ =-(\frac{1}{1+e^{-XW}})^{-2}.(1+e^{-XW})^{'} \\ =-(\frac{1}{1+e^{-XW}})^{-2}.e^{-XW}.(-XW)^{'} \\ =(\frac{1}{1+e^{-XW}})^{-2}.e^{-XW}.X_i \\ =p^{'}(1-p^{'})X_i (p′)′=(1+e−XW1)′=−(1+e−XW1)−2.(1+e−XW)′=−(1+e−XW1)−2.e−XW.(−XW)′=(1+e−XW1)−2.e−XW.Xi=p′(1−p′)Xi
当我们采用整个矩阵进行求导的时候，我发现上面得出的矩阵结果维度好像不多，但是有一点就是通过下面的代码可以看出(−XW)′(-XW)^{'}(−XW)′对wiw_iwi求导确实只剩下XiX_iXi，此时向量对于变量进行求导，采用分子布局方式。
[1x11x12⋯x1p1x21x22⋯x2p⋯1xn1xn2⋯xnp]∗[w0w1⋯wi⋯wp]=[y0y1⋯yi⋯yn]\left[ \begin{array}{c} 1 & x_{11} & x_{12} & \cdots & x_{1p} \\ 1 & x_{21} & x_{22} & \cdots & x_{2p}\\ \cdots \\ 1 & x_{n1} & x_{n2} & \cdots & x_{np} \end{array} \right] * \left[ \begin{array}{c} w_0 \\ w_1 \\ \cdots \\ w_i \\ \cdots \\ w_p \end{array} \right] = \left[ \begin{array}{c} y_0 \\ y_1 \\ \cdots \\ y_i \\ \cdots \\ y_n \end{array} \right] ⎣⎢⎢⎡11⋯1x11x21xn1x12x22xn2⋯⋯⋯x1px2pxnp⎦⎥⎥⎤∗⎣⎢⎢⎢⎢⎢⎢⎡w0w1⋯wi⋯wp⎦⎥⎥⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎡y0y1⋯yi⋯yn⎦⎥⎥⎥⎥⎥⎥⎤
因此，我们不对整个输入值矩阵进行求导，我们只针对单个输入值进行求导，我们采用下面的求导方式，然后观察出规律即可。假设输入值为x(i)x^{(i)}x(i)，逻辑回归函数对wiw_iwi进行求导

(p)′=(11+e−WTx(i))′=−(11+e−WTx(i))−2.(1+e−WTx(i))′=−(11+e−WTx(i))−2.e−WTx(i).(−WTx(i))′=(11+e−WTx(i))−2.e−WTx(i).xi(i)=p(1−p)xi(i)(p)^{'}=(\frac{1}{1+e^{-W^Tx^{(i)}}})^{'} \\ =-(\frac{1}{1+e^{-W^Tx^{(i)}}})^{-2}.(1+e^{-W^Tx^{(i)}})^{'} \\ =-(\frac{1}{1+e^{-W^Tx^{(i)}}})^{-2}.e^{-W^Tx^{(i)}}.(-W^Tx^{(i)})^{'} \\ =(\frac{1}{1+e^{-W^Tx^{(i)}}})^{-2}.e^{-W^Tx^{(i)}}.x^{(i)}_i \\ =p(1-p)x^{(i)}_i (p)′=(1+e−WTx(i)1)′=−(1+e−WTx(i)1)−2.(1+e−WTx(i))′=−(1+e−WTx(i)1)−2.e−WTx(i).(−WTx(i))′=(1+e−WTx(i)1)−2.e−WTx(i).xi(i)=p(1−p)xi(i)

损失函数F(W)对W进行求导：
F(W)′=(∑i=1N(yiln(p)+(1−yi)ln(1−p)))′=∑i=1N(yiln(p)′+(1−yi)ln(1−p)′)=∑i=1N(yi(1pp′)+(1−yi)ln(1−p)′)=∑i=1N(yi(1pp′)+(1−yi)(11−p(1−p)′))=∑i=1N(yi−p)xi(i){F(W)}^{'}=(\sum_{i=1}^N({{y_i}ln(p)+({1-y_i})ln(1-p)}))^{'} \\ =\sum_{i=1}^N({{y_i}ln(p)^{'}+({1-y_i})ln(1-p)^{'}}) \\ =\sum_{i=1}^N({{y_i}(\frac{1}{p}p^{'})+({1-y_i})ln(1-p)^{'}}) \\ =\sum_{i=1}^N({{y_i}(\frac{1}{p}p^{'})+({1-y_i}) (\frac{1}{1- p}{(1-p)^{'}})}) \\ =\sum_{i=1}^N(y_i - p)x^{(i)}_i F(W)′=(i=1∑N(yiln(p)+(1−yi)ln(1−p)))′=i=1∑N(yiln(p)′+(1−yi)ln(1−p)′)=i=1∑N(yi(p1p′)+(1−yi)ln(1−p)′)=i=1∑N(yi(p1p′)+(1−yi)(1−p1(1−p)′))=i=1∑N(yi−p)xi(i)
最后采用下面的公式，在代码中更新权重即可：
wik+1=wik+α∗∂F(w)∂wiw^{k+1}_i = w^{k}_i + \alpha*\frac{\partial {F(w)}}{\partial {w_i}} wik+1=wik+α∗∂wi∂F(w)

2、代码

1、算法描述
本博客采用逻辑回归算法实现鸢尾花数据集的二分类。
2、数据csv准备
（1）数据格式

sepallength  sepalwidth  petallength petalwidth  class
5.1 3.5 1.4 0.2 Iris-setosa
4.9 3   1.4 0.2 Iris-setosa
4.7 3.2 1.3 0.2 Iris-setosa
4.6 3.1 1.5 0.2 Iris-setosa
5   3.6 1.4 0.2 Iris-setosa
5.4 3.9 1.7 0.4 Iris-setosa
4.6 3.4 1.4 0.3 Iris-setosa
5   3.4 1.5 0.2 Iris-setosa
4.4 2.9 1.4 0.2 Iris-setosa
...

为了大家能够更好的学习KNN算法，我将这个csv文件放到我博客的资源中，通过下面的链接即可下载：

https://download.csdn.net/download/weixin_43334389/13130463

3、代码

# @Time : 2020/11/22 20:45
# @Description : 使用鸢尾花数据集采用逻辑回归算法实现二分类
import numpy as np
import pandas as pddata = pd.read_csv(r"dataset/iris.arff.csv", header=0)
data.drop_duplicates(inplace=True)
# 名称映射为数字
data["class"] = data["class"].map({"Iris-setosa": 0, "Iris-versicolor": 1, "Iris-virginica": 2})
# 筛选数据，只选择0和1
data = data[data["class"] != 2]class LogisticRegression:"""@desc:逻辑回归算法实现"""def __init__(self, learning_rate, times):"""@desc:初始化算法@param learning_rate: 学习率@param times: 迭代次数"""self.learning_rate = learning_rateself.times = timesdef sigmoid(self, z):"""@desc:sigmoid函数实现@param z: 自变量： z=w.T*x@return: float, 属于[0,1]区间。s>=0.5(z>0)时，类别为1，否则为0"""return 1.0 / (1.0 + np.exp(-z))def fit(self, X, y):"""@desc:训练模型@param X:训练集特征矩阵，shape[0]为样本数量（即特征矩阵），shape[1]为特征数量@param y:标签数组"""X = np.asarray(X)print(X.shape)print(type(X))y = np.asarray(y)# 创建权重向量，初始化为0.多出来的1是截距self.w_ = np.zeros(1 + X.shape[1])  # shape[0]为样本数量，shape[1]为特征数量# 创建损失列表，保存每次迭代后的损失值self.loss_ = []for i in range(self.times):z = np.dot(X, self.w_[1:]) + self.w_[0]# 计算概率值(判定为1的概率值)p = self.sigmoid(z)# 根据逻辑回归的代价函数（目标函数或叫损失函数）计算损失值# 通过计算逻辑回归函数的最大似然函数得到下面的公式，np.log()函数以e为底# J(W) = - sum(yi * log(s(zi))+(1-yi)log(1-s(zi)))cost = - np.sum(y * np.log(p) + (1 - y) * np.log(1 - p))self.loss_.append(cost)# 调整权重值，上面的损失函数对w求导得到# 公式 权重(j行) = 权重(j行) + 学习率 * sum( y - s(z)) * X(j列)self.w_[0] += self.learning_rate * np.sum(y - p) * 1self.w_[1:] += self.learning_rate * np.dot(X.T, y - p)def predict_proba(self, X):"""@desc: 根据参数传递的样本，对样本数据进行预测属于某一个类别的概率@param X:测试集特征矩阵@return:预测（概率值）"""X = np.asarray(X)z = np.dot(X, self.w_[1:]) + self.w_[0]p = self.sigmoid(z)# 将预测结果转换为二维数组，方便后续拼接p = p.reshape(-1, 1)# 将两个数组拼接，方向是横向拼接test1 = np.concatenate([1 - p, p], axis=1)print("--------test1----------")# (18, 2)print(test1.shape)# <class 'numpy.ndarray'>print(type(test1))# [[0.99066619 0.00933381]#  [0.99297624 0.00702376]#  [0.97928108 0.02071892]print(test1)return np.concatenate([1 - p, p], axis=1)  # 横向拼接为为1def predict(self, X):"""@desc:根据参数预测样本属于哪个类别@param X:测试集特征矩阵@return: 返回测试数据属于那个类别标签对应的较大值的索引axis=1 表示每一行的最大索引，当下面函数输出为0的时候，表示第一列数据大，刚好对应分类为0的标签。"""return np.argmax(self.predict_proba(X), axis=1)t1 = data[data["class"] == 0]
t2 = data[data["class"] == 1]
t1 = t1.sample(len(t1), random_state=0)
t2 = t2.sample(len(t2), random_state=0)
train_X = pd.concat([t1.iloc[:40, :-1], t2.iloc[:40, :-1]], axis=0)
train_y = pd.concat([t1.iloc[:40, -1], t2.iloc[:40, -1]], axis=0)
test_X = pd.concat([t1.iloc[40:, :-1], t2.iloc[40:, :-1]], axis=0)
test_y = pd.concat([t1.iloc[40:, -1], t2.iloc[40:, -1]], axis=0)# 鸢尾花的特征列都在同一个数量级，我们可以不用标准化处理
lr = LogisticRegression(learning_rate=0.01, times=20)
lr.fit(train_X, train_y)
# 预测的概率值
# lr.predict_proba(test_X)
# lr.predict(test_X)预测的类别
result = lr.predict(test_X)
# 计算准确性
accuracy = np.sum(result == test_y) / len(test_y)
# 1.0
print(accuracy)import matplotlib as mpl
import matplotlib.pyplot as pltmpl.rcParams['font.family'] = 'SimHei'
mpl.rcParams['axes.unicode_minus'] = False# 绘制预测值
plt.plot(result, 'ro', ms=15, label="预测值")  # ms指定圆圈大小
plt.plot(test_y.values, 'go', label="真实值")  # pandas读取时serise类型，我们需要转为ndarray
plt.title('逻辑回归预测')
plt.xlabel('样本序号')
plt.ylabel('预测值')
plt.legend()
plt.show()# 绘制目标函数随着迭代次数的损失值图像
plt.plot(range(1, lr.times + 1), lr.loss_, 'go-')
plt.xlabel('迭代次数')
plt.ylabel('损失值loss')
plt.show()

4、结果展示
（1）预测值和真实值之间的拟合

（2）损失函数图像

3、分析

1、argmax函数

 def predict(self, X):#索引值和分类类别相对应return np.argmax(self.predict_proba(X), axis=1)# 举例
matrix2_2 = np.array([[1, 2], [3, 4]])
# [1 1] 即第一行的2和第二行的4比较大，作为索引返回。
# axis=1表示返回每1行的最大索引(0/1).
# 参数axis,默认是0.
print(np.argmax(matrix2_2, axis=1))

2、数组拼接函数concatenate

#拼接之后成为一个n*2的数组矩阵，其索引值（0/1）和分类类别相对应
return np.concatenate([1 - p, p], axis=1)  # 横向拼接为1# 举例 按轴axis连接array组成一个新的array
# 一维情况进行连接 axis=1
# 我们对其进行转换尺寸
a1 = a.reshape(-1, 1)
b1 = b.reshape(-1, 1)
ab1 = np.concatenate((a1, b1), axis=1)
# [[ 1 11]
#  [ 2 22]
#  [ 3 33]]
print(ab1)
(3, 2)
print((ab1).shape)# 二维情况进行连接 axis=0
c = np.array([[1, 2], [3, 4]])
d = np.array([[5, 6]])
cd = np.concatenate((c, d), axis=0)
# [[1 2]
#  [3 4]
#  [5 6]]
print(cd)# 二维情况进行连接 axis=0
cd2 = np.concatenate((c, d.T), axis=1)
# [[1 2 5]
#  [3 4 6]]
print(cd2)