吴恩达《机器学习》课后测试Ex2:逻辑回归(详细Python代码注解)
基于吴恩达《机器学习》课程
参考黄海广的笔记
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import scipy.optimize as opt # 用于高级优化
Part1
在第一部分,我们将要构建一个逻辑回归模型来预测,某个学生是否被大学录取。
设想你是大学相关部分的管理者,想通过申请学生两次测试的评分,来决定他们是否被录取。现在你拥有之前申请学生的可以用于训练逻辑回归的训练样本集。对于每一个训练样本,你有他们两次测试的评分和最后是被录取的结果。为了完成这个预测任务,我们准备构建一个可以基于两次测试评分来评估录取可能性的分类模型。
sigmoid 函数
hθ(x)=11+e−θTXh_\theta \left( x \right)=\frac{1}{1+{{e}^{-\theta^{T}X}}}hθ(x)=1+e−θTX1
# 实现sigmoid函数
def sigmoid(z):return 1 / (1 + np.exp(-z))
做一个快速的检查,来确保它可以工作。
# 仅用于测试sigmoid函数运行
def sig_test():nums = np.arange(-10, 10, step=1) # 以-10为起点,10为终点,步长为1的列表fig, ax = plt.subplots(figsize=(12, 8)) # 以其他关键字参数**fig_kw来创建图ax.plot(nums, sigmoid(nums), 'r')plt.show()
代价函数
J(θ)=−1m∑i=1m[y(i)log(hθ(x(i)))+(1−y(i))log(1−hθ(x(i)))]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)]}J(θ)=−m1i=1∑m[y(i)log(hθ(x(i)))+(1−y(i))log(1−hθ(x(i)))]
def cost(theta, X, y):theta = np.matrix(theta) # 将theta转换为矩阵X = np.matrix(X) # 将X转换为矩阵y = np.matrix(y) # 将y转换为矩阵first = np.multiply(-y, np.log(sigmoid(X * theta.T))) # 公式第一项 np.multiply将数组或矩阵对应位置相乘。second = np.multiply((1 - y), np.log(1 - sigmoid(X * theta.T))) # 公式第二项return np.sum(first - second) / (len(X))
first
和second
都是一个矩阵。
计算梯度∂∂θjJ(θ)\frac{\partial }{\partial {\theta_j}}J\left( \theta \right)∂θj∂J(θ)
注意,我们实际上没有在这个函数中执行梯度下降,我们仅仅在计算一个梯度步长。梯度下降可以参考Ex1中的代码。
这次,我们不自己写代码实现梯度下降,我们会调用一个已有的库。这就是说,我们不用自己定义迭代次数和步长,功能会直接告诉我们最优解。Andrew NG在课程中用的是Octave的“fminunc”函数,由于我们使用Python,我们可以用scipy.optimize.fmin_tnc做同样的事情。
# 实现梯度计算(并没有更新θ),没有梯度下降
def gradient(theta, X, y):theta = np.matrix(theta) # 将theta转换为矩阵X = np.matrix(X) # 将X转换为矩阵y = np.matrix(y) # 将y转换为矩阵temp = np.matrix(np.zeros(theta.shape)) # np.zeros(theta.shape)=[0.,0.],然后将temp变为矩阵[0.,0.]parameters = int(theta.ravel().shape[1])# theta.ravel():将多维数组theta降为一维,.shape[1]是统计这个一维数组有多少个元素。parameters表示参数grad = np.zeros(parameters)error = sigmoid(X * theta.T) - yfor i in range(parameters):term = np.multiply(error, X[:, i]) # 将误差与训练数据相乘,得到偏导数grad[i] = np.sum(term) / len(X)return grad
这里的error
就是δ=a−y=∂C∂z\delta=a-y=\frac{\partial C}{\partial {z}}δ=a−y=∂z∂C,a=g′(z)a=g'(z)a=g′(z),X[:, j]
则是∂z∂w=a\frac{\partial z}{\partial {w}}=a∂w∂z=a,所以term
为∂C∂w=∂z∂w∂C∂z\frac{\partial C}{\partial {w}}=\frac{\partial z}{\partial {w}}\frac{\partial C}{\partial {z}}∂w∂C=∂w∂z∂z∂C,即为代价函数偏导数。可参考笔记反向传播详解。
评价逻辑回归模型的预测
predictA
:预测学生是否录取。如果一个学生exam1得分45,exam2得分85,那么他录取的概率应为0.776
predictB
:看模型在训练集上的正确率怎样。给出数据以及参数后,predictB的函数会返回“1”或者“0”。然后再把这个predict函数用于训练集上,看准确率怎样。
def predictA(theta, X):return sigmoid(X @ theta)def predictB(theta, X):probability = sigmoid(X * theta.T)return [1 if x >= 0.5 else 0 for x in probability] #大于0.5的数取值为1,小于0.5的为0
这里注意 *
和@
的区别,*
星乘表示矩阵内各对应位置相乘,@
等价于.dot()
点乘表示求矩阵内积。
逻辑回归代码
此处代码是连续的,中途输出和一些说明单独拿出来了。
def ex2_1():# 从检查数据开始data1 = pd.read_csv('ex2data1.txt', names=['text1', 'text2', 'admitted'])print(data1.head())print(data1.describe()) # 查看数据的基本情况,包括:count非空值数、mean平均值、std 标准差、max、min、(25%、50%、75%)分位数# 让我们创建两个分数的散点图,并使用颜色编码来可视化,正的(被接纳)或负的(未被接纳)positive = data1[data1['admitted'].isin([1])] # isin()为筛选函数,令positive为数据集中为1的数组negative = data1[data1['admitted'].isin([0])] # isin()为筛选函数,令positive为数据集中为0的数组fig, ax = plt.subplots(figsize=(12, 8)) # 以其他关键字参数**fig_kw来创建图# figsize=(a,b):figsize 设置图形的大小,b为图形的宽,b为图形的高,单位为英寸ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='admitted')# x为Exam 1,y为Exam 2,散点大小s设置为50,颜色参数c为蓝色,散点形状参数marker为圈,以关键字参数**kwargs来设置标记参数labele是Admittedax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='not admitted')# x为Exam 1,y为Exam 2,散点大小s设置为50,颜色参数c为红色,散点形状参数marker为叉,以关键字参数**kwargs来设置标记参数labele是Not Admittedax.legend() # legend为显示图例函数ax.set_xlabel('text 1 Score') # 设置x轴变量ax.set_ylabel('text 2 Score') # 设置y轴变量plt.show()
此处输出如下:
text1 text2 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
text1 text2 admitted
count 100.000000 100.000000 100.000000
mean 65.644274 66.221998 0.600000
std 19.458222 18.582783 0.492366
min 30.058822 30.603263 0.000000
25% 50.919511 48.179205 0.000000
50% 67.032988 67.682381 1.000000
75% 80.212529 79.360605 1.000000
max 99.827858 98.869436 1.000000
# 代码接上# 初始化X,y,thetadata1.insert(0, 'Ones', 1) # 在第0列插入表头为“ONE”的列,数值为1cols = data1.shape[1] # 获取表格df的列数X = data1.iloc[:, 0:cols - 1] # 除最后一列外,取其他列的所有行,即X为O和人口组成的列表y = data1.iloc[:, cols - 1:cols] # 取最后一列的所有行,即y为利润X = np.array(X.values) # 将X的各个值转换为数组y = np.array(y.values) # 将y的各个值转换为数组theta = np.zeros(3) # 初始化theta数组为(0,0,0)# 检查矩阵的维度print(X.shape, theta.shape, y.shape)# 用初始θ计算代价print(cost(theta, X, y))# 看看梯度计算结果print(gradient(theta, X, y))
输出如下:
(100, 3) (3,) (100, 1)
0.6931471805599453
[ -0.1 -12.00921659 -11.26284221]
# 代码接上# 用SciPy's truncated newton(TNC)实现寻找最优参数。result = opt.fmin_tnc(func=cost, x0=theta, fprime=gradient, args=(X, y))print(result)# 看看在这个结论下代价函数计算结果是什么个样子print(cost(result[0], X, y)) # result[0]:取数组result中的第一个元素
result
返回值说明:数组:返回优化问题目标值;nfeval整数:function evaluations的数目;rc。
在进行优化的时候,每当目标优化函数被调用一次,就算一个function evaluation。在一次迭代过程中会有多次function evaluation,这个参数不等同于迭代次数,而往往大于迭代次数。
输出结果如下:
(array([-25.16131863, 0.20623159, 0.20147149]), 36, 0)
0.20349770158947458
# 代码接上# 画出决策曲线plotting_x1 = np.linspace(30, 100, 100)plotting_h1 = (- result[0][0] - result[0][1] * plotting_x1) / result[0][2]# theta^T*x=0为判定边界,x0theta0+x1theta1+x2theta2=0可得:x2=(-x0theta0-x1theta1)/theta2fig, ax = plt.subplots(figsize=(12, 8))ax.plot(plotting_x1, plotting_h1, 'y', label='Prediction')ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='Admitted')ax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='Not Admitted')ax.legend()ax.set_xlabel('Exam 1 Score')ax.set_ylabel('Exam 2 Score')plt.show()
θTx=0{\theta^{T}}x=0θTx=0 这条线便是判定边界(decision boundary)。
输出图像:
# 代码接上# predictA预测学生是否录取。print(predictA(result[0], [1, 45, 85]))# predictB预测准确率。theta_min = np.matrix(result[0]) # 将最优化最优参数中优化问题目标值数组x转化为矩阵predictions = predictB(theta_min, X)correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y)]# predictions为预测值,y为实际值.zip为打包对象函数,将(predictions, y)打包成元组,并赋值为(a,b),然后判断正确的则返回1,判断错误的返回0accuracy = (sum(map(int, correct)) % len(correct))# map(int, correct):将crrect列表内容的类型映射成int型,原来应该布尔型.sum为求和函数,%为求模运算,计算除法的余数。# 实际是求1占数据总数的比例print('accuracy = {0}%'.format(accuracy)) # format为格式化字符串的函数,{0}为设置指定位置.
predictB
这是训练集的准确性。我们没有保持住了设置或使用交叉验证得到的真实逼近,所以这个数字有可能高于其真实值。
输出结果如下:
0.7762906240463825
accuracy = 89%
Part2
在第二部分,我们将实现加入正则项提升逻辑回归算法。
设想你是工厂的生产主管,你有一些芯片在两次测试中的测试结果,测试结果决定是否芯片要被接受或抛弃。你有一些历史数据,帮助你构建一个逻辑回归模型。
正则化的代价函数
J(θ)=1m∑i=1m[−y(i)log(hθ(x(i)))−(1−y(i))log(1−hθ(x(i)))]+λ2m∑j=1nθj2J\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}}J(θ)=m1i=1∑m[−y(i)log(hθ(x(i)))−(1−y(i))log(1−hθ(x(i)))]+2mλj=1∑nθj2
# 正则化的代价函数
def costreg(theta, X, y, learningRate):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)))reg = (learningRate / (2 * len(X)) * np.sum(np.power(theta[:, 1:theta.shape[1]], 2)))return np.sum(first - second) / len(X) + reg
learningRate
就是正则化参数λ\lambdaλ。
正则化的梯度计算
我们未对θ0\theta_0θ0进行正则化,所以梯度下降算法应分两种情形:
RepeatRepeatRepeat untiluntiluntil convergenceconvergenceconvergence{
θ0:=θ0−a1m∑i=1m((hθ(x(i))−y(i))x0(i)){\theta_0}:={\theta_0}-a\frac{1}{m}\sum\limits_{i=1}^{m}{(({h_\theta}({{x}^{(i)}})-{{y}^{(i)}})x_{0}^{(i)}})θ0:=θ0−am1i=1∑m((hθ(x(i))−y(i))x0(i))
θj:=θj−a[1m∑i=1m((hθ(x(i))−y(i))xj(i)+λmθj]{\theta_j}:={\theta_j}-a[\frac{1}{m}\sum\limits_{i=1}^{m}{(({h_\theta}({{x}^{(i)}})-{{y}^{(i)}})x_{j}^{\left( i \right)}}+\frac{\lambda }{m}{\theta_j}]θj:=θj−a[m1i=1∑m((hθ(x(i))−y(i))xj(i)+mλθj] (forforfor j=1,2,...nj=1,2,...nj=1,2,...n)
}
# 正则化的梯度计算,没有梯度下降
def gradientReg(theta, X, y, learningRate):theta = np.matrix(theta)X = np.matrix(X)y = np.matrix(y)parameters = int(theta.ravel().shape[1])#theta.ravel():将多维数组转换为一维数组,并不会产生源数据的副本. shape[1]:获取列数grad = np.zeros(parameters) #创建零数组error = sigmoid(X * theta.T) - yfor i in range(parameters):term = np.multiply(error, X[:,i]) # term即偏导数if (i == 0): # 不对第0个参数正则化grad[i] = np.sum(term) / len(X)else:grad[i] = (np.sum(term) / len(X)) + ((learningRate / len(X)) * theta[:,i])return grad
这里的error
就是δ=a−y=∂C∂z\delta=a-y=\frac{\partial C}{\partial {z}}δ=a−y=∂z∂C,a=g′(z)a=g'(z)a=g′(z),X[:, j]
则是∂z∂w=a\frac{\partial z}{\partial {w}}=a∂w∂z=a,所以term
为∂C∂w=∂z∂w∂C∂z\frac{\partial C}{\partial {w}}=\frac{\partial z}{\partial {w}}\frac{\partial C}{\partial {z}}∂w∂C=∂w∂z∂z∂C,即为代价函数偏导数。可参考笔记反向传播详解。
正则化的逻辑回归代码
此处代码是连续的,中途输出和一些说明单独拿出来了。
def ex2_2():# 第一部分很像,从数据可视化开始data2 = pd.read_csv('ex2data2.txt', names=['text1', 'text2', 'accepted'])print(data2.head())positive = data2[data2['accepted'].isin([1])] # isin()为筛选函数,令positive为数据集中为1的数组negative = data2[data2['accepted'].isin([0])] # isin()为筛选函数,令positive为数据集中为0的数组fig, ax = plt.subplots(figsize=(12, 8)) # 以其他关键字参数**fig_kw来创建图ax.scatter(positive['text1'], positive['text2'], s=50, c='b', marker='o', label='accepted')# x为Exam 1,y为Exam 2,散点大小s设置为50,颜色参数c为蓝色,散点形状参数marker为圈,以关键字参数**kwargs来设置标记参数labele是Admittedax.scatter(negative['text1'], negative['text2'], s=50, c='r', marker='x', label='not accepted')# x为Exam 1,y为Exam 2,散点大小s设置为50,颜色参数c为红色,散点形状参数marker为叉,以关键字参数**kwargs来设置标记参数labele是Not Admittedax.legend() # 显示图例ax.set_xlabel('text1 Score') # 设置x轴变量ax.set_ylabel('text2 Score') # 设置y轴变量plt.show()
输出如下:
text1 text2 accepted
0 0.051267 0.69956 1
1 -0.092742 0.68494 1
2 -0.213710 0.69225 1
3 -0.375000 0.50219 1
4 -0.513250 0.46564 1
以上图片显示,这个数据集不能像之前一样使用直线将两部分分割。而逻辑回归只适用于线性的分割,所以,这个数据集不适合直接使用逻辑回归。
一种更好的使用数据集的方式是为每组数据创造更多的特征(特征映射):
# 代码接上degree = 5x1 = data2['text1']x2 = data2['text2']data2.insert(3, 'Ones', 1) # 在第4列插入标题为One,各行值为1的列# 这里是把输入x1,x2重新映射成5-1=4阶了,x0=1,x1,x1^2,x1x2,x1^3,x1^2x2,x1x2^2...for i in range(1, degree):for j in range(0, i):data2['F' + str(i) + str(j)] = np.power(x1, i - j) * np.power(x2, j)# np.power(A,B) ## 对A中的每个元素求B次方data2.drop('text1', axis=1, inplace=True) # 删除表头为Test 1的列,axis=1表示默认删除行或列,inplace=True表示原数组被data2替换.data2.drop('text2', axis=1, inplace=True) # 删除表头为Test 2的列,axis=1表示默认删除行或列,inplace=True表示原数组被data2替换.print(data2.head())
输出如下:
accepted Ones F10 F20 ... F40 F41 F42 F43
0 1 1 0.051267 0.002628 ... 0.000007 0.000094 0.001286 0.017551
1 1 1 -0.092742 0.008601 ... 0.000074 -0.000546 0.004035 -0.029801
2 1 1 -0.213710 0.045672 ... 0.002086 -0.006757 0.021886 -0.070895
3 1 1 -0.375000 0.140625 ... 0.019775 -0.026483 0.035465 -0.047494
4 1 1 -0.513250 0.263426 ... 0.069393 -0.062956 0.057116 -0.051818
[5 rows x 12 columns]
接下来就像在第一部分中做的一样,初始化变量。
# 代码接上cols = data2.shape[1] # 获取data2的列数X2 = data2.iloc[:, 1:cols] # 取所有列每一行的数据y2 = data2.iloc[:, 0:1] # 取第一和第二列每一行的数据X2 = np.array(X2.values) # 将数值转换为数组y2 = np.array(y2.values)theta2 = np.zeros(11) # 创建长度为11的零数组# 正则化参数取合理值,我们可以之后再折腾这个。learningRate = 1print(costreg(theta2, X2, y2, learningRate)) # 用初始参数计算代价print(gradientReg(theta2, X2, y2, learningRate)) # 看看梯度计算结果
输出如下:
0.6931471805599454
[0.00847458 0.01878809 0.05034464 0.01150133 0.01835599 0.007323930.00819244 0.03934862 0.00223924 0.01286005 0.00309594]
# 代码接上
# 使用和第一部分相同的优化函数来计算优化后的结果。result2 = opt.fmin_tnc(func=costreg, x0=theta2, fprime=gradientReg, args=(X2, y2, learningRate))print(result2)# 使用第1部分中的预测函数predictB来查看我们的方案在训练数据上的准确度。theta_min = np.matrix(result2[0])predictions = predictB(theta_min, X2)correct = [1 if ((a == 1 and b == 1) or (a == 0 and b == 0)) else 0 for (a, b) in zip(predictions, y2)]# predictions为预测值,y为实际值.zip为打包对象函数,将(predictions, y)打包成元组,并赋值为(a,b),然后判断正确的则返回1,判断错误的返回0accuracy = (sum(map(int, correct)) % len(correct))# map(int, correct):将crrect列表内容的类型映射成int型,原来应该布尔型.sum为求和函数,%为求模运算,计算除法的余数。# 实际是求1占数据总数的比例print('accuracy = {0}%'.format(accuracy))
输出如下:
(array([ 0.53010247, 0.29075567, -1.60725764, -0.58213819, 0.01781027,-0.21329508, -0.40024142, -1.3714414 , 0.02264304, -0.9503358 ,0.0344085 ]), 22, 1)
accuracy = 78%
result
返回值说明:数组:返回优化问题目标值;nfeval整数:function evaluations的数目;rc。
后面还有“画出决策的曲线”和“改变λ,观察决策曲线”两部分,但是好像和**支持向量机(SVM)**有关,就先不深究了。
吴恩达《机器学习》课后测试Ex2:逻辑回归(详细Python代码注解)相关推荐
- 吴恩达机器学习作业2:逻辑回归(Python实现)
逻辑回归 在训练的初始阶段,将要构建一个逻辑回归模型来预测,某个学生是否被大学录取.设想你是大学相关部分的管理者,想通过申请学生两次测试的评分,来决定他们是否被录取.现在你拥有之前申请学生的可以用于训 ...
- 逻辑回归python sigmoid(z)_python实现吴恩达机器学习练习2(逻辑回归)-data1
python实现吴恩达机器学习练习2(逻辑回归)-data1 这篇是第一个数据集:这部分练习中,你将建立一个预测学生是否被大学录取的逻辑回归模型. 假如一所大学会每个报名学生进行两项入学考试,根据两项 ...
- 吴恩达机器学习课程-作业1-线性回归(python实现)
Machine Learning(Andrew) ex1-Linear Regression 椰汁学习笔记 最近刚学习完吴恩达机器学习的课程,现在开始复习和整理一下课程笔记和作业,我将陆续更新. Li ...
- 吴恩达机器学习(四)逻辑回归(二分类与多分类)
目录 0. 前言 1. 假设函数(Hypothesis) 2. 决策边界(Decision Boundary) 3. 代价函数(Cost Funciton) 4. 梯度下降(Gradient Desc ...
- 吴恩达机器学习第二次作业——逻辑回归
逻辑回归 一.逻辑回归 1,数据可视化 2,sigmoid函数,逻辑回归模型 3,代价函数以及梯度 4,评价逻辑回归 二.正规化逻辑回归 1,数据可视化 2,特征映射(Feature Mapping) ...
- 吴恩达机器学习课后习题ex2
ex2:Logistic_regression import numpy as np import pandas as pd import matplotlib.pyplot as plt impor ...
- 吴恩达机器学习笔记三之逻辑回归
本节目录: 分类问题 假说表示 判定边界 代价函数 高级优化 多类别分类 1.分类问题 在分类问题中,我们尝试预测的是结果是否属于某一个类(例如正确或错误).分类问 题的例子有:判断一封电子邮件是否是 ...
- 目录:吴恩达机器学习课后作业
简单介绍 本博客为作者自行完成的吴恩达机器学习课后练习题目录,均使用PyTorch进行实现,无法保证对错,仅供参考. 作业题目以及源代码 百度云盘连接 提取码:3dvb 题目的命名方式与下表中的作业名 ...
- 2.吴恩达机器学习课程-作业2-逻辑回归
fork了别人的项目,自己重新填写,我的代码如下 https://gitee.com/fakerlove/machine-learning/tree/master/code 代码原链接 文章目录 2. ...
- 1. 吴恩达机器学习课程-作业1-线性回归
fork了别人的项目,自己重新填写,我的代码如下 https://gitee.com/fakerlove/machine-learning/tree/master/code 代码原链接 文章目录 1. ...
最新文章
- Java项目:前后端分离疫情防疫平台设计和实现(java+springmvc+VUE+node.js+mybatis+mysql+springboot+redis+jsp)
- python菜鸟工具-Python3 教程
- [其实有加强版的]校门外的树
- ffmpeg为AVPacket添加解码头信息
- Dubbo的RPC原理
- openai-gpt_为什么到处都看到GPT-3?
- 家庭上网用路由器和ADSL的连接
- 每年考证时间表(绝对会用得到的一天,怕到时候很难找到有用) ——自己留住,哦!!!!
- 【数理统计】卡方检验
- 微软升级网页版Skype 没有帐户也能拨打网络电话
- dolphinscheduler 3.0.1 资源中心
- 图解Topo拓扑排序
- 清橙A1210. 光棱坦克
- STM32内部ADC测量时产生噪声的原因与消除的方法
- Jenkins打包部署项目到Windows或Linux运行
- 比Siri更厉害的个人助理Viv 能否一统江湖?
- C语言顺序表:1、顺序表的存储、2、顺序表的实现.
- Jaccard和Levenshtein
- 拼写检查器的编写[转]
- StreamWriter