编程作业（python）| 吴恩达 机器学习（2）逻辑回归

8道编程作业及解析见：Coursera吴恩达机器学习编程作业

1. 线性可分

题目

根据学生的两门成绩，预测该学生是否会被大学录取，数据集：ex2data1.txt

（100行×3列，第一列为Exam1成绩，第二列为Exam2成绩，第三列0表示不及格，1表示及格）

1.1 导入数据

• 导入库：

  import numpy as np  			# 科学计算库，处理多维数组，进行数据分析
  import pandas as pd				# 基于 NumPy 的一种解决数据分析任务工具
  import matplotlib.pyplot as plt # 提供一个类似 Matlab 的绘图框架
  

• 导入数据集：

  path = 'ex2data1.txt'
  data = pd.read_csv(path, names=['Exam 1', 'Exam 2', 'Accepted'])
  data.head()
  

  

• 数据可视化

  fig,ax = plt.subplots()
  # 绘制标签为 0 的散点图
  ax.scatter(data[data['Accepted']==0]['Exam 1'], data[data['Accepted']==0]['Exam 2'],
  		   c='r',marker='x',label='y=0')
  # 绘制标签为 1 的散点图		   
  ax.scatter(data[data['Accepted']==1]['Exam 1'], data[data['Accepted']==1]['Exam 2'],
  		   c='b',marker='o',label='y=1')
  ax.legend()
  ax.set(xlabel='exam1', ylabel='exam2')
  plt.show()
  

1.2 构造数据集

• 封装一个get_Xy()函数来获取输入向量X和输出向量y以及初始化参数向量theta：

  	def get_Xy_theta(data):
  		data.insert(0,'$x_0$',1)
   		cols = data.shape[1]
   	
  		X_ = data.iloc[:,0:cols-1]
  		X = X_.values
  	
  		y_ = data.iloc[:,cols-1:cols]
  		y = y_.values   
  	
  		theta = np.zeros((X.shape[1], y.shape[1]))
  	
  		return X,y,theta
  
  	X,y,theta = get_Xy_theta(data)
  

  可以用 X.shape 、 y.shape 和theta.shape查看数组的维度 ，维度分别为（100，3）、（100，1）和（3，1）

1.3 代价函数

• 假设函数：$h_\theta(x)=g(\theta X)$
• 逻辑函数：$g(z)=\dfrac {1}{1+e^{-z}}$
• 代价函数：$J(θ)=-\displaystyle\frac{1}{m} SUM\left[y*\log(h_\theta(x))+(1-y)*\log(1-h_\theta(x))\right]$$SUM$ 表示对向量的所有项求和，最终得一个标量值。公式推导见 3 代价函数

  # 逻辑函数
  def sigmoid(z):    
  	return 1 / (1 + np.exp(-z))
  # 代价函数
  def costFunction(X,y,theta):
  	A = sigmoid(X@theta)# 假设函数
  	first = y * np.log(A)
  	second = (1-y) * np.log(1-A)
  	return -np.sum(first + second) / len(X) #SUM表示对向量的所有项求和
  

1.4 梯度下降函数

$\theta=\theta-\alpha\frac{1}{m}X^T\left [ g(X\theta)-y \right]$公式推导见： 1.4 梯度下降函数

• 定义梯度下降函数：

  def gradientDescent(X,y,theta,iters,alpha):
  	m = len(X)
  	costs = []
      
  	for i in range(iters):
      	A = sigmoid(X@theta)
      	theta  = theta -(alpha/m)  * X.T @ (A - y)
      	cost = costFunction(X,y,theta)
      	costs.append(cost)
      	if i % 10000 == 0:
         		print(cost) 
          
  	return costs, theta
  

• 给学习率alpha 和迭代次数iters赋上初值，就能得到迭代结束后的最终参数theta_final

  alpha = 0.004
  iters = 300000
  
  costs, theta_final = gradientDescent(X,y,theta,iters,alpha) 
  

1.5 绘制决策边界

相关概念见 3.决策边界

• 要绘制出决策线，先求出直线的两个参数$-\frac{\theta_0}{\theta_2}$与$-\frac{\theta_1}{\theta_2}$

  para1 = - theta_final[0,0] / theta_final[2,0]
  para2 = - theta_final[1,0] / theta_final[2,0]
  

• 结果可视化

  fig,ax = plt.subplots(figsize=(7,5))
  
  #绘制直线
  x = np.linspace(20,100,2)#在20到100之前绘制2个点确定直线，注意纵轴为Exam2，刻度在20-100
  y_ = para1 + para2 * x
  ax.plot(x,y_,c='g',label='Decision Boundary')
  
  #绘制原始数据
  ax.scatter(data[data['Accepted']==0]['Exam 1'], data[data['Accepted']==0]['Exam 2'],c='r',marker='x',label='y=0')		   
  ax.scatter(data[data['Accepted']==1]['Exam 1'], data[data['Accepted']==1]['Exam 2'],c='b',marker='o',label='y=1')
  
  ax.legend()#显示标签label
  ax.set(xlabel='exam1', ylabel='exam2')
  plt.show()
  

1.6 计算分类正确率

从上图中可以看出来此线性的决策边界并不能将样本点完全区分开，部分点有误分情况。

我们知道在逻辑回归中，假设函数表示 “预测值=1” 的概率（见 2.假设函数）：

$h_θ(x)=g(\theta X)=P(y=1|x;\theta)$

• 定义假设函数，获取假设输出向量：

  def predict(X,theta):
  	prob = sigmoid(X@theta)  
  	return [1 if result >= 0.5 else 0 for result in prob]#返回一个列表
  
  h_x = predict(X,theta_final) 		#在参数为theta_final下的假设输出，为列表类型
  h_x = np.array(h_x)					#将列表转化为一维数组
  h_x = y_pre.reshape(len(h_x ),1)	#进而转化为二维数组，因为真实输出y也是二维数组，方便比较
  

  注意函数predict的返回值其实是一个将假设输出$h_\theta(x)$归类的过程：

  • 当 $h_θ(x)\geq0.5$ ，即 $\theta X\geq0$ 时，将其预测归为1 ，表现为上图决策边界上方的点
  • 当 $h_θ(x)<0.5$ ，即 $\theta X<0$ 时，将其预测归为0，表现为上图决策边界下方的点

• 计算正确率：

  acc = np.mean(h_x == y)#将假设输出和真实输出进行比较，求平均值
  print(acc)
  

  可以得到在参数为final_theta条件下，该模型的正确率为 0.89，由于总样本数为100个，因此有11个样本点分类错误，这时你再往上看决策界面，看看是不是这样。

2. 线性不可分

题目

作为工厂的生产主管，你要决定某芯片是否要被接受或抛弃，该芯片在两次测试中的结果保存在数据集：ex2data2.txt

（118行×3列，第一列为Test1结果，第二列为Test2结果，第三列0表示不合格，1表示合格）

2.1 导入数据

• 导入库：

  import numpy as np  			# 科学计算库，处理多维数组，进行数据分析
  import pandas as pd				# 基于 NumPy 的一种解决数据分析任务工具
  import matplotlib.pyplot as plt # 提供一个类似 Matlab 的绘图框架
  

• 导入数据集：

  path = 'ex2data2.txt'
  data = pd.read_csv(path, names=['Test 1', 'Test 2', 'Accepted'])
  data.head()
  

  

• 数据可视化

  fig,ax = plt.subplots()
  # 绘制标签为 0 的散点图
  ax.scatter(data[data['Accepted']==0]['Test 1'], data[data['Accepted']==0]['Test 2'],
  		   c='r',marker='x',label='y=0')
  # 绘制标签为 1 的散点图		   
  ax.scatter(data[data['Accepted']==1]['Test1'], data[data['Accepted']==1]['Test 2'],
  		   c='b',marker='o',label='y=1')
  ax.legend()
  ax.set(xlabel='exam1', ylabel='exam2')
  plt.show()
  

  

2.2 特征映射 Feature Mapping

从上图可以看出，该逻辑回归的样本数量多，且分布是非线性的，而原始特征只有$x1,x2$，给分类带来不小难度。

解决办法是用多项式创建更多的特征$x1、x2、x1x2、x1^2、x2^2、... x_1^nx_2^n$。因为更多的特征进行逻辑回归时，得到的决策边界可以是任意高阶函数的形状。

def feature_mapping(x1,x2,power):
    data = {} # 字典 保存特征多项式
    for i in np.arange(power+1):  # i的范围是(0,power]，整数
        for j in np.arange(i + 1):# j的范围是(0,i]，整数
            data['F{}{}'.format(i-j,j)] = np.power(x1,i-j) * np.power(x2,j)    
    return pd.DataFrame(data) 

以power = 2为例：

• i = 0, j = 0：字典的键值有$F00=x_1^0x_2^0=1$
• i = 1
  • j = 0：$F10=x_1^1x_2^0=x_1$
  • j = 1：$F10=x_1^0x_2^1=x_2$
• i = 2
  • j = 0：$F20=x_1^2x_2^0=x_1^2$
  • j = 1：$F11=x_1^1x_2^1=x_1x_2$
  • j = 2：$F02=x_1^0x_2^2=x_2^2$

2.3 构造数据集

• 获取特征映射后的特征集

  # 获取原始特征作为 feature_mapping 的参数
  x1 = data['Test 1']
  x2 = data['Test 2']
  # 得到特征映射后的特征集
  power = 6
  data_map = feature_mapping(x1,x2,power)
  

   以power = 6为例，得到的特征从 F00 到 F06，即 $x_1^0x_2^0$ 到 $x_2^6$，共28列：

• 获取输入输出向量 $X,y$ 及参数向量 $\theta$

  X = data_map .values
  
  cols = data.shape[1]
  y = data.iloc[:,cols-1:cols]
  y = y.values   
  	
  theta = np.zeros((X.shape[1], y.shape[1]))
  

 用 X.shape 、 y.shape 和theta.shape查看数组的维度 ，维度分别为（118，28）、（118，1）和（28，1）

注意：经过特征映射后的$X$已经具有$x_1^0x_2^0=1=x_0$这一项，因此不需要再单独插入$x_0=1$

2.4 代价函数

由于特征映射后维数较高，模型更复杂，容易产生过拟合现象，因此需要对代价函数正则化(Regularize)==（或称 加入惩罚项），$\lambda$ 为正则化系数：

注意惩罚项中不包括$\theta_0$，公式推导见 2.2 正则化的逻辑回归

• 包含正则化项的代价函数如下：

  # 逻辑函数
  def sigmoid(z):    
  	return 1 / (1 + np.exp(-z))
  
  # 代价函数
  def costFunction(X,y,theta):
  	A = sigmoid(X@theta)# 假设函数
  	
  	first = y * np.log(A)
  	second = (1-y) * np.log(1-A)
  	reg = np.sum(np.power(theta[1:],2)) * (lamda / (2 * len(X)))#注意正则化项的θ是从j=1开始
  	
  	return -np.sum(first + second) / len(X) + reg
  

2.5 梯度下降函数

梯度下降公式分两种情况：

可以用向量统一表示为（最右边一项即为正则化项）：
$\theta=\theta-\alpha\frac{1}{m}X^T\left [ g(X\theta)-y \right]-\alpha\frac{\lambda}{m}\theta^{reg}，其中\theta^{reg}_0=0，其余项与\theta相等$

• 定义梯度下降函数：

  def gradientDescent(X,y,theta,alpha,iters,lamda):
  
  costs = []
  
  for i in range(iters):
      
      reg = alpha * theta[1:] * (lamda / len(X)) # 先截取除了 θ_0 的 θ，此时reg维度为(27,1)
      reg = np.insert(reg,0,values=0,axis=0)     # 再在首项插入 0，使reg维度为(28,1)
      
      theta = theta - (X.T@(sigmoid(X@theta) - y)) * alpha / len(X) -reg
      cost = costFunction(X,y,theta,lamda)
      costs.append(cost)
      
      if i % 10000 == 0:
          print(cost)
          
  return theta,costs
  

• 给学习率alpha 、迭代次数iters和正则化系数lamda赋上初值，就能得到迭代结束后的最终参数theta_final，例如：

  alpha = 0.001
  iters = 200000
  lamda = 0.001
  costs, theta_final = gradientDescent(X,y,theta,iters,alpha) 
  

2.6 绘制决策边界（欠拟合与过拟合）

• 结果可视化

  fig,ax = plt.subplots()
  
  # 原始数据
  ax.scatter(data[data['Accepted']==0]['Test 1'],data[data['Accepted']==0]['Test 2'],c='r',marker='x',label='y=0')
  ax.scatter(data[data['Accepted']==1]['Test 1'],data[data['Accepted']==1]['Test 2'],c='b',marker='o',label='y=1')
  
  # 绘制决策边界
  a = np.linspace(-1.2,1.2,200)
  xx,yy = np.meshgrid(a,a) #生成200×200的网格采样点
  z = feature_mapping(xx.ravel(),yy.ravel(),6).values # 特征映射得到多维特征 z ,ravel()将多维数组转换为一维数组
  
  z = x @ theta_final
  zz = zz.reshape(xx.shape)
  plt.contour(xx,yy,zz,0)# Xθ = 0即为决策边界，contour()用于绘制等高线
  
  ax.legend()
  ax.set(xlabel='Test1',ylabel='Test2')
  
  plt.show()
  

  

• alpha 和iters不变，改变 正则化系数$\lambda$ ：

  • lamda = 10得到分类结果如下，分类正确率为：acc = 0.746（代码与#1.6一样，就不码了，注意此时样本个数为118个）（欠拟合）
    
  • lamda = 5得到分类结果如下，分类正确率为：acc = 0.814（欠拟合）
    
  • lamda = 1得到分类结果如下，分类正确率为：acc = 0.831（分类效果：优）
    
  • lamda = 0得到分类结果如下，分类正确率为：acc = 0.813（过拟合）
    

由上述结果可以看出：

• 当$\lambda$过大时，容易产生欠拟合
• 当$\lambda$过小时，容易产生过拟合（本例不太好展示，若样本点不多或者正负样本交叉较多，过拟合效果会很明显）

转载：https://blog.csdn.net/m0_37867091/article/details/104887979