逻辑回归算法实验
【实验目的】
理解逻辑回归算法原理，掌握逻辑回归算法框架；
理解逻辑回归的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）输出分类结果的混淆矩阵。
【实验报告要求】
1.对照实验内容，撰写实验过程、算法及测试结果；
2.代码规范化：命名规则、注释；
3.实验报告中需要显示并说明涉及的数学原理公式；
4.查阅文献，讨论逻辑回归算法的应用场景；
sigmoid函数
Sigmoid函数是一个在生物学中常见的S型函数，也称为S型生长曲线。在深度学习中，由于其单增以及反函数单增等性质，Sigmoid函数常被用作神经网络的激活函数，将变量映射到[0,1] [0,1][0,1]之间。
Sigmoid函数的特性与优缺点：
（1）Sigmoid函数的输出范围是0到1。由于输出值限定在0到1，因此它对每个神经元的输出进行了归一化。
（2）用于将预测概率作为输出的模型。由于概率的取值范围是0到1，因此Sigmoid函数非常合适梯度平滑，避免跳跃的输出值函数是可的。
（3）这意味着可以找到任意两个点的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。我们只需要将转换后的结果进行归一化处理。方法就是将转化后的结果除以所有转化后结果之和，可以理解为转化后结果占总数的百分比。这样就得到近似的概率。
假如模型对一个三分类问题的预测结果为4、5、6。我们要用softmax将模型结果转为概率。步骤如下：
y1=exp(4)=20.08
y2=exp(5)=148.4
y3=exp(6)=403.4 sum=y1+y2+y3=571.88
z1=y1/sum=0.035
z2=y2/sum=0.259
z3=y3/sum=0.705 # 注意这是近似值
(1)读取数据
import numpy as np import pandas as pd #读取数据 data=pd.read_csv("E:\homework/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
1 def sigmoid(x):
2 return 1/(1+np.exp(-x))
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时，计算此时的代价值
(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时，计算一下此时的梯度为多少
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]
在求得最优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))
(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()
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()
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()
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
（3）输出分类结果的混淆矩阵
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()
