吴恩达机器学习第二次编程作业-逻辑回归

机器学习

课时60编程作业：逻辑回归

工具：pycharm2019.1，Python3.5

参考资料： Logistic回归总结、梯度下降求解逻辑回归 、数据提取、Python 不同颜色可视化数据集、Dataframe筛选数据、DataFrame类型转换成array类型 、pandas数据读写

逻辑回归完整代码
正则化完整代码

1. 逻辑回归

任务：建立一个逻辑回归分类模型来预测一个学生是否能被大学录取

1.1 数据可视化

先利用pandas库提取数据，然后分类，绘制。
pandas需要手动添加一行表示列名，不然以第一行数据作为列名，数据类型自动识别。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
def pandas_getData():
    data = pd.read_csv('./ex2/ex2data1.txt')
    y = data['y']   # 按列名提取
    data1 = data[y == 0]  # 数据分类
    data2 = data[y == 1]  # 已录取
    return data,data1,data2,data.shape[0]

def plot_data(data1,data2):
    plt.plot(data1['x1'],data1['x2'],'yo',label='Not admitted')
    plt.plot(data2['x1'],data2['x2'],'k+',label='Admitted')
    plt.xlim(30,101)
    plt.ylim(30,100)
    plt.xlabel('Exam 1 score',fontsize=10)
    plt.ylabel('Exam 2 score', fontsize=10)
    plt.legend(loc=1,fontsize=8)
    plt.show()

Figure_1.png

1.2 Cost function and gradient

2019-07-16 14-44-13 的屏幕截图.png

关于为什么学习率α这里不用除以m均值化在参考资料中有回答说α是个常数，可以省掉m，本人不是很懂。
代价函数与梯度下降只要有线形回归的基础就很好计算：

1
2
3
4
5
6
7
8
9
10
def Cost_Function(x,y,theta,m):
    hypothesis=1/(1+np.exp(-np.dot(x,theta)))   # 假设函数
    Cost_J =-1/m*(y*np.log(hypothesis)+(1-y)*np.log(1-hypothesis)).sum()
    return Cost_J

def Gradient_descent(x,y,theta,m,alpha,iterations):
    for iteration in range(iterations):
        hypothesis=1/(1+np.exp(-np.dot(x,theta)))
        theta=theta-alpha*np.dot(x.T,hypothesis-y)
    return theta   # 我们只需要最终的参数theta来进行预测

用梯度下降的方法来计算参数θ还需要初始化参数init_theta、α以及迭代次数MaxIter，也别忘了X需要增加一列1表示x0。

1
2
3
4
5
6
7
8
9
data,data1,data2,m=pandas_getData()  # 获取数据
# plot_data(data1,data2)  # 绘制数据
x=np.array(data[['x1','x2']])  # 将DataFrame类型转换成array类型
y=np.array(data['y']).reshape(m,1)
x=np.hstack((np.ones((m,1)),x))   # 增加一列 X0=1
init_theta=np.zeros((x.shape[1],1)) # 初始化theta
MaxIter=2000000  # 迭代次数
alpha=0.00001 # 学习率
print('初始代价函数为：',Cost_Function(x,y,init_theta,m)) # 这里如果代价函数没写错计算出来的值应该是0.6931471805599453

  • 在源文档中给出的求参数的方法是用Octave调库直接得出，网上找到的很多博客也是用scipy库函数。毫无疑问库中封装了最简介高效的方法。但如果我就想使用梯度下降呢？
  • 本人已经试过了，需要不停的调参（α和MaxIter），这对于机器性能是个考验。用库函数得出的最优参数是：[-25.16131867, 0.20623159, 0.20147149]，用这个参数向量求出的代价函数值为：0.2034977015894744。
  • 本人在不停的调参过程中发现当我把学习率α设置为0.00001时，再一直增大迭代次数发现实验结果越来越接近文档给出的最优值，考虑到机器性能我只增加到了200000次，梯度下降得出的参数θ向量为：[-19.08468795 0.15767916 0.15228738]，用此参数求出的代价函数值为：0.21041324594328958
  • 有了参数θ就可以用来进行预测和绘制决策边界线了
预测和决策边界线
  1. 还记得之前说过hθ(x)=1/(1+e-θTX))，这个函数用来预测分类y=1概率。
  2. 我们可以设定阈值为0.5，只要h>=0.5，那么就判定这是属于1，否则属于0。
  3. 上面是具体计算方法，我们可以有两种形式来预测：一种是用给定的样本X看看其h函数概率；还有一种是直接用训练样本跑，看有多少个正确分类了，这个百分比越高说明预测越准确。
  4. 以上是预测方法，那么关于决策边界线，这个在吴恩达机器学习视频中已经给出。没有实践很容易忽略这点，我们要观察h函数，它是一个sigma函数，只要θTX的值大于等于0，那么h函数值就大于等于0.5，想想看是不是。而X向量是样本点，这条分界线就是若干个样本点，这些样本点所计算出来的θTX=θ0+θ1X1+… 值应该等于0，两侧的点要么大于0计算出来的h值就大于0.5，反之同理。
  5. 用文档给出的样本[1 45 85]预测的概率为0.7221802，文档给出的概率为0.776。看来梯度下降算法确实不如其他算法，唯一的优点可能就是简单了。

  6. 先用训练样本来看看实验效果如何： 89.00%

    1
    2
    3
    4
    5
    6
    7
    8
    def text_sample(x,y,theta,m):
        hypothesis=1/(1+np.exp(-np.dot(x,theta)))
        cnt_correct=0
        for i in range(m):
            if hypothesis[i][0]>=0.5 and y[i][0]==1:cnt_correct=cnt_correct+1 # 正确预测1
            if hypothesis[i][0]<0.5 and y[i][0]==0:cnt_correct=cnt_correct+1 # 正确预测0
        correct_percent=cnt_correct/m  # 计算准确率
        print('%.2f%%' %(correct_percent*100))

  7. 绘制决策边界，我们已经知道决策边界线的θTX=θ0+θ1X1+…θnXn=0，数据图的横轴是X1，纵轴是X2，那么随便造几百个X1，根据此公式计算X2再绘制直线即可。

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    def plotDecisionBoundary(data1,data2,x,y,theta):
        plt.plot(data1['x1'], data1['x2'], 'yo', label='Not admitted')
        plt.plot(data2['x1'], data2['x2'], 'k+', label='Admitted')
        plt.xlim(30, 101)
        plt.ylim(30, 100)
        plt.xlabel('Exam 1 score', fontsize=10)
        plt.ylabel('Exam 2 score', fontsize=10)
        plt.legend(loc=1, fontsize=8)

        x1 = np.arange(130, step=0.1)
        x2 = -(final_theta[0] + x1 * final_theta[1]) / final_theta[2]
        plt.plot(x1, x2)
        plt.show()

Figure_2.png

至此，简单逻辑回归分类模型已经完成，下面是第二部分-正则化。

2. 正则化

芯片质检问题，两轮检测，一个结果，建立逻辑回归模型。

2.1 绘制数据

和1.1一样，使用Pandas库提取数据，然后按y=0/1分类，传参绘制

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
import matplotlib
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt

def Get_data():
    data=pd.read_csv('./ex2/ex2data2.txt',names=['x1','x2','y']) # setting header or default setting header=None
    data1=data[data['y']==1] # accepted
    data2=data[data['y']==0] # rejected
    print(data1.shape,data2.shape)
    return data,data1,data2

def Plot_Data(data1,data2):
    plt.plot(data1['x1'],data1['x2'],'k+',label='y=1')
    plt.plot(data2['x1'],data2['x2'],'yo',label='y=0')
    plt.xlim(-1,1.5)
    plt.ylim(-0.8,1.2)
    plt.xlabel('Microchip Test 1',fontsize=10)
    plt.ylabel('Microchip Test 2',fontsize=10)
    plt.legend(loc='upper right',fontsize=8)
    plt.show()
data,data1,data2=Get_data()
Plot_Data(data1,data2)

Figure_3.png

2.2 Feature mapping 特征扩展

特征越多，维数越高，所作出的曲线拟合得也更好（第8章课程总结第8页）。

1
2
3
4
5
6
7
def mapFeature(x1, x2,m):
    out=np.ones((m,1))
    degree = 6
    for i in range(1,degree+1):
        for j in range(0,i+1):
            out=np.hstack((out,((x1**(i-j))*(x2**j)).reshape(m,1)))
    return out

2.3 Cost function and gradient

2019-07-19 15-57-31 的屏幕截图.png
2019-07-19 15-57-35 的屏幕截图.png
根据上面的公式就可以计算了。
先初始化：

1
2
3
4
5
6
7
8
9
data,data1,data2,m=Get_data()
# Plot_Data(data1,data2)
X=np.array(data[['x1','x2']])
X=mapFeature(X.T[0],X.T[1],m)
Y=np.array(data['y']).reshape(m,1)
init_theta=np.zeros((X.shape[1],1))
lamda=1
alpha=0.0001
MaxIter=2000000

代价函数：

1
2
3
4
5
def Cost_Function(X,Y,init_theta,m,lamda):
    hypothesis=1/(1+np.exp(-X.dot(init_theta)))
    Cost_J=-1/m*(Y*np.log(hypothesis)+(1-Y)*np.log(1-hypothesis)).sum()
    Cost_J=Cost_J+lamda/(2*m)*(init_theta**2).sum()
    return Cost_J

梯度下降：

1
2
3
4
5
6
def Gradient_descent(X,Y,theta,m,alpha,iterations,lamda):
    for i in range(iterations):
        tmp=theta[0][0] # theta0单独处理
        theta=theta-alpha/m*(np.dot(X.T,(1/(1+np.exp(-X.dot(theta)))-Y))+lamda*theta)
        theta[0][0]=tmp+alpha*lamda/m
    return theta

2.4 绘制决策边界线

这里是用等高线的方法作决策边界线，随机生成x1和x2，然后再用特征扩展方法扩展维数与θ向量相乘。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
def plot_Decision_boundary(data1,data2,theta):
    plt.plot(data1['x1'], data1['x2'], 'k+', label='y=1')
    plt.plot(data2['x1'], data2['x2'], 'yo', label='y=0')
    plt.xlim(-1, 1.5)
    plt.ylim(-0.8, 1.2)
    plt.xlabel('Microchip Test 1', fontsize=10)
    plt.ylabel('Microchip Test 2', fontsize=10)
    plt.title('lambda = 1',fontsize=15)
    point_num=50
    x1=np.linspace(-1,1.5,point_num)
    x2=np.linspace(-1,1.5,point_num)
    z=np.zeros((point_num,point_num))
    for i in range(point_num):
        for j in range(point_num):
            z[i][j]=mapFeature(x1[i],x2[j],1).dot(theta)
    plt.contour(x1,x2,z)
    plt.legend(loc='upper right', fontsize=8)
    plt.show()

Figure_4.png

由于这里的参数θ以及学习率分别用梯度下降跑出来的和自己设置的，与用库函数跑出来的最优值不同，所以作出的效果与文档也有差别。