逻辑回归算法实验

【实验目的】

理解逻辑回归算法原理，掌握逻辑回归算法框架；
理解逻辑回归的sigmoid函数；
理解逻辑回归的损失函数；
针对特定应用场景及数据，能应用逻辑回归算法解决实际分类问题。

【实验内容】

1.根据给定的数据集，编写python代码完成逻辑回归算法程序，实现如下功能：

建立一个逻辑回归模型来预测一个学生是否会被大学录取。假设您是大学部门的管理员，您想根据申请人的两次考试成绩来确定他们的入学机会。您有来自以前申请人的历史数据，可以用作逻辑回归的训练集。对于每个培训示例，都有申请人的两次考试成绩和录取决定。您的任务是建立一个分类模型，根据这两门考试的分数估计申请人被录取的概率。
算法步骤与要求：

(1)读取数据；(2)绘制数据观察数据分布情况；(3)编写sigmoid函数代码；(4)编写逻辑回归代价函数代码；(5)编写梯度函数代码；(6)编写寻找最优化参数代码（可使用scipy.opt.fmin_tnc()函数）；(7)编写模型评估（预测）代码，输出预测准确率；(8)寻找决策边界，画出决策边界直线图。

2. 针对iris数据集，应用sklearn库的逻辑回归算法进行类别预测。

要求：

（1）使用seaborn库进行数据可视化；（2）将iri数据集分为训练集和测试集(两者比例为8:2)进行三分类训练和预测；（3）输出分类结果的混淆矩阵。

【实验报告要求】

对照实验内容，撰写实验过程、算法及测试结果；
代码规范化：命名规则、注释；
实验报告中需要显示并说明涉及的数学原理公式；
查阅文献，讨论逻辑回归算法的应用场景；

sigmoid函数

Sigmoid函数是一个在生物学中常见的S型函数，也称为S型生长曲线。在深度学习中，由于其单增以及反函数单增等性质，Sigmoid函数常被用作神经网络的激活函数，将变量映射到[0,1] [0,1][0,1]之间。

Sigmoid函数的特性与优缺点：

Sigmoid函数的输出范围是0到1。由于输出值限定在0到1，因此它对每个神经元的输出进行了归一化。

用于将预测概率作为输出的模型。由于概率的取值范围是0到1，因此Sigmoid函数非常合适梯度平滑，避免跳跃的输出值函数是可的。

这意味着可以找到任意两个点的Sigmoid曲线的斜率明确的预测，即非常接近1或0。函数输出不是以0为中心的，这会降低权重更新的效率Sigmoid函数执行指数运算，计算机运行得较慢。

import matplotlib.pyplot as plt
import numpy as np
import math
x = np.linspace(-10, 10, 100)
z = 1 / (1 + np.exp(-x))
plt.title("Sigmoid")
plt.plot(x, z)
plt.xlabel("x")
plt.ylabel("Sigmoid(X)")
plt.show()

损失函数

每一个样本经过模型后会得到一个预测值，然后得到的预测值和真实值的差值就成为损失。损失值越小证明模型越好。

那我们为什么需要损失函数呢！误差值越小，模型就越小。我们想让预测值无限接近于真实值，所以需要将误差值降到最低。因此在这个过程中就需要引入损失函数。

绝对值损失函数

平方损失函数

softmax函数

称归一化指数函数。它是二分类函数sigmod在多分类上的推广，目的是将多分类的结果以概率的形式展现出来。

softmax是怎么实现的呢！首先我们知道概率性质：非负数，并且概率和为1。softmax就是将在负无穷到正无穷上的预测结果按照概率的性质转换为概率。

1.首先先将结果转换成非负数，指数的值域是0到正无穷，softmax第一步就是把预测的结果映射到指数上。

2.为了确保各个预测结果的概率之和等于1。我们只需要将转换后的结果进行归一化处理。方法就是将转化后的结果除以所有转化后结果之和，可以理解为转化后结果占总数的百分比。这样就得到近似的概率。

$$\text{假如模型对一个三分类问题的预测结果为4、5、6。我们要用softmax将模型结果转为概率。步骤如下：}$$

$$\text{y1=exp(4)=20.08}$$

$$\text{y2=exp(5)=148.4}$$

$$\text{y3=exp(6)=403.4 sum=y1+y2+y3=571.88}$$

$$\text{z1=y1/sum=0.035}$$

$$\text{z2=y2/sum=0.259}$$

$$\text{z3=y3/sum=0.705 注意这是近似值}$$

逻辑回归算法原理，掌握逻辑回归算法框架

建立一个逻辑回归模型来预测一个学生是否会被大学录取

(1)读取数据

 

import numpy as np
import pandas as pd
#读取数据
data=pd.read_csv("E:\jupyter_root_directory\data/ex2data1.txt",header=None,names=['grade1','grade2','Admitted'])
data

	grade1	grade2	Admitted
0	34.623660	78.024693	0
1	30.286711	43.894998	0
2	35.847409	72.902198	0
3	60.182599	86.308552	1
4	79.032736	75.344376	1
...	...	...	...
95	83.489163	48.380286	1
96	42.261701	87.103851	1
97	99.315009	68.775409	1
98	55.340018	64.931938	1
99	74.775893	89.529813	1

100 rows × 3 columns

(2)绘制数据观察数据分布情况

from pyecharts.charts import Scatter
import pyecharts.options as opts
admittedData=data[data['Admitted'].isin([1])]
noAdmittedData=data[data['Admitted'].isin([0])]
c=(
    Scatter()
    .add_xaxis(admittedData['grade1'])
    .add_yaxis("grade1",noAdmittedData['grade1'],label_opts=opts.LabelOpts(is_show=False))
    .add_xaxis(admittedData['grade2'])
    .add_yaxis("grade2",noAdmittedData['grade2'],label_opts=opts.LabelOpts(is_show=False))
)
c.render_notebook()

#在逻辑回归模型中，x0=1,即训练数据应该添加一列，值为1
data.insert(0, 'ones',1)
loc=data.shape[1]   
X=np.array(data.iloc[:,0:loc-1])     #取前三列
Y=np.array(data.iloc[:,loc-1:loc])   #取最后一列
T=np.zeros(X.shape[1]) 
X.shape,Y.shape,T.shape

 

(3)编写sigmoid函数代码

 def sigmoid(x):

return 1/(1+np.exp(-x)) 

(4)编写逻辑回归代价函数代码

代价函数的公式

 

def computeCost(theta,X,Y):
    theta = np.matrix(theta) #矩阵对象, 不能缺少，因为参数theta是一维数组，进行矩阵想乘时要把theta先转换为矩阵 创建了一个新的相同的矩阵。当修改新矩阵时，原来的矩阵不会改变。
    h=sigmoid(np.dot(X,(theta.T)))
    a=np.multiply(-Y,np.log(h)) #矩阵对应元素相乘  此处的np.log为数学上的ln 
    b=np.multiply((1-Y),np.log(1-h)) 
    return np.sum(a-b)/len(X)
computeCost(T,X,Y)  #当theta值为0时，计算此时的代价值

结果为：0.6931471805599453

(5)编写梯度函数代码

def gradient(theta,X,Y):
    theta = np.matrix(theta) #要先把theta转化为矩阵
    h=sigmoid(np.dot(X,(theta.T)))
    grad=np.dot(((h-Y).T),X)/len(X)
    return np.array(grad).flatten()  #因为下面寻找最优化参数的函数（opt.fmin_tnc())要求传入的gradient函返回值需要是一维数组，因此需要利用flatten（）将grad进行转换以下

gradient(T,X,Y) #测试一下，当T值都为为0时，计算一下此时的梯度为多少

结果：array([-12.00921659, -11.26284221])

(6)编写寻找最优化参数代码（可使用scipy.opt.fmin_tnc()函数）

在实现线性回归时，是利用梯度下降的方式来寻找最优参数。在此处使用scipy.optimize包下的fmin_tnc函数来求解最优参数，该函数利用截断牛顿算法中的梯度信息，最小化具有受边界约束的变量的函数。

import scipy.optimize as opt
result = opt.fmin_tnc(func=computeCost, x0=T, fprime=gradient, args=(X, Y)) #func：优化的目标函数 （在这里要优化的是代价函数）
# x0：初始值，必须是一维数组  fprime：提供优化函数func的梯度函数，不然优化函数func必须返回函数值和梯度，或者设置approx_grad=True （在这里梯度函数是gradient函数，并且要求返回的是一维数组）
# args：元组，是传递给优化函数的参数
# 返回值： x ： 数组，返回的优化问题目标值 （在这里即优化后，theta的最终取值） 
# nfeval：整数，功能评估的数量。在进行优化的时候，每当目标优化函数被调用一次，就算一个function evaluation。在一次迭代过程中会有多次function evaluation。这个参数不等同于迭代次数，而往往大于迭代次数。
# rc: int，返回码
print(result)
theta=result[0]

结果：(array([0.01044469, 0.00043394]), 11, 1)

(7)编写模型评估（预测）代码，输出预测准确率

在求得最优theta值后，利用得到的模型在训练数据中进行预测，并求准确率。
由逻辑回归的假设模型可知：
当hθ(x)>=0.5时，预测y=1;
当hθ(x)<0.5时，预测y=0;

predict函数：通过训练数据以及T值进行预测，并且把预测结果使用列表返回；
hypothesis=[1 if a==b else 0 for (a,b)in zip(predictValues,Y)] 目的是将预测值与实际值进行比较，如果二者相等，则为1，否则为0；
accuracy=hypothesis.count(1)/len(hypothesis) 计算hypothesis中1的个数然后除以总的长度，得到准确率

def predict(theta, X):
    theta = np.matrix(theta)
    temp = sigmoid(X * theta.T)
    #print(temp)
    return [1 if x >= 0.5 else 0 for x in temp]
    
predictValues=predict(theta,X)
hypothesis=[1 if a==b else 0 for (a,b)in zip(predictValues,Y)]
accuracy=hypothesis.count(1)/len(hypothesis)
print ('accuracy = {0}%'.format(accuracy*100))

结果：accuracy = 60.0%

(8)寻找决策边界，画出决策边界直线图。

 

#决策边界
def find_x2(x1,theta):
    return [(-theta[0]-theta[1]*x_1)/theta[2] for x_1 in x1]
x1 = np.linspace(30, 100, 1000)
x2=find_x2(x1,theta)

#数据可视化
admittedData=data[data['isAdmitted'].isin([1])]
noAdmittedData=data[data['isAdmitted'].isin([0])]
fig,ax=plt.subplots(figsize=(12,8))
ax.scatter(admittedData['exam1'],admittedData['exam2'],marker='+',label='addmitted')
ax.scatter(noAdmittedData['exam2'],noAdmittedData['exam1'],marker='o',label="not addmitted")
ax.plot(x1,x2,color='r',label="decision boundary")
ax.legend(loc=1)
ax.set_xlabel('Exam1 score')
ax.set_ylabel('Exam2 score')
ax.set_title("Training data with decision boundary")
plt.show()

 

 

2.针对iris数据集，应用sklearn库的逻辑回归算法进行类别预测

（1）使用seaborn库进行数据可视化

from sklearn.datasets import load_iris    #数据模块
import seaborn as sns
import matplotlib.pyplot as plt
iris = load_iris()   #加载数据
data = pd.DataFrame(iris.data, columns=iris.feature_names)   #整理数据
g = sns.pairplot(data, diag_kind="kde",kind="reg")
#kind:用于控制非对角线上图的类型，可选'scatter'与'reg'
#diag_kind:用于控制对角线上的图分类型，可选'hist'与'kde'
#kind='scatter'时就相当于原图
#可以看到对角线上是各个属性的直方图（分布图），而非对角线上是两个不同属性之间的相关图
plt.show()

（2）将iri数据集分为训练集和测试集(两者比例为8:2)进行三分类训练和预测

from sklearn.linear_model import LogisticRegression
from sklearn.model_selection import train_test_split
from sklearn.datasets import load_iris
data = load_iris()
iris_target = data.target  # 获取y
iris_features = pd.DataFrame(data=data.data, columns=data.feature_names)
# 测试集大小为20%
x_train, x_test, y_train, y_test = train_test_split(iris_features, iris_target, test_size = 0.2, random_state = 2020)
## 定义 逻辑回归模型 
clf = LogisticRegression(random_state=0, solver='lbfgs')  #拟牛顿法的一种，利用损失函数二阶导数矩阵即海森矩阵来迭代优化损失函数
# 在训练集上训练逻辑回归模型
clf.fit(x_train, y_train)
## 查看其对应的w
print('the weight of Logistic Regression:\n',clf.coef_)
 
## 查看其对应的w0
print('the intercept(w0) of Logistic Regression:\n',clf.intercept_)
 
## 由于这个是3分类，所有我们这里得到了三个逻辑回归模型的参数，其三个逻辑回归组合起来即可实现三分类。

结果：

the weight of Logistic Regression:
 [[-0.45928925  0.83069887 -2.26606531 -0.99743981]
 [ 0.33117319 -0.72863424 -0.06841147 -0.9871103 ]
 [ 0.12811606 -0.10206464  2.33447678  1.98455011]]
the intercept(w0) of Logistic Regression:
 [  9.4388067    3.93047364 -13.36928034]

from sklearn import metrics
## 在训练集和测试集上分布利用训练好的模型进行预测
train_predict = clf.predict(x_train)
test_predict = clf.predict(x_test)
 
## 由于逻辑回归模型是概率预测模型（前文介绍的 p = p(y=1|x,\theta)）,所有我们可以利用 predict_proba 函数预测其概率
train_predict_proba = clf.predict_proba(x_train) # 第一个值表示预测为0的概率，第二个值表示预测为1的概率，并且每行的概率值之和都为1
test_predict_proba = clf.predict_proba(x_test)    
 
print('The test predict Probability of each class:\n',test_predict_proba)
# 其中第一列代表预测为0类的概率，第二列代表预测为1类的概率，第三列代表预测为2类的概率。
 
# 利用accuracy（准确度）预测正确的样本数目占总预测样本数目的比例评
print('The accuracy of the Logistic Regression is:',metrics.accuracy_score(y_train,train_predict))
print('The accuracy of the Logistic Regression is:',metrics.accuracy_score(y_test,test_predict))

The test predict Probability of each class:
 [[1.03461737e-05 2.33279477e-02 9.76661706e-01]
 [9.69926591e-01 3.00732874e-02 1.21677000e-07]
 [2.09992549e-02 8.69156616e-01 1.09844129e-01]
 [3.61934872e-03 7.91979966e-01 2.04400686e-01]
 [7.90943209e-03 8.00605299e-01 1.91485269e-01]
 [7.30034956e-04 6.60508053e-01 3.38761912e-01]
 [1.68614211e-04 1.86322045e-01 8.13509341e-01]
 [1.06915331e-01 8.90815532e-01 2.26913671e-03]
 [9.46928071e-01 5.30707288e-02 1.20016060e-06]
 [9.62346385e-01 3.76532228e-02 3.91897297e-07]
 [1.19533386e-04 1.38823469e-01 8.61056998e-01]
 [8.78881880e-03 6.97207359e-01 2.94003822e-01]
 [9.73938143e-01 2.60617342e-02 1.22613839e-07]
 [1.78434056e-03 4.79518177e-01 5.18697483e-01]
 [5.56924345e-04 2.46776840e-01 7.52666235e-01]
 [9.83549842e-01 1.64500666e-02 9.13617272e-08]
 [1.65201476e-02 9.54672748e-01 2.88071041e-02]
 [8.99853722e-03 7.82707575e-01 2.08293888e-01]
 [2.98015029e-05 5.45900069e-02 9.45380192e-01]
 [9.35695863e-01 6.43039522e-02 1.85301368e-07]
 [9.80621190e-01 1.93787398e-02 7.00125265e-08]
 [1.68478817e-04 3.30167227e-01 6.69664294e-01]
 [3.54046168e-03 4.02267804e-01 5.94191734e-01]
 [9.70617284e-01 2.93824735e-02 2.42443971e-07]
...
 [9.64848137e-01 3.51516747e-02 1.87917886e-07]
 [9.70436779e-01 2.95624021e-02 8.18591621e-07]]
The accuracy of the Logistic Regression is: 0.9833333333333333
The accuracy of the Logistic Regression is: 0.8666666666666667

关于函数predict和predict_proba的区别可以参考这个博客

https://blog.csdn.net/qq_43468807/article/details/105740396

（3）输出分类结果的混淆矩阵

混淆矩阵的理解可以参考

https://blog.csdn.net/SartinL/article/details/105844832

 

from sklearn import metrics
## 查看混淆矩阵
confusion_matrix_result = metrics.confusion_matrix(test_predict,y_test)
print('The confusion matrix result:\n',confusion_matrix_result)
 
# 利用热力图对于结果进行可视化
plt.figure()
sns.heatmap(confusion_matrix_result, annot=True, cmap='Blues')
plt.xlabel('Predicted labels')
plt.ylabel('True labels')
plt.show()

The confusion matrix result:
 [[10  0  0]
 [ 0  8  2]
 [ 0  2  8]]