（西瓜书）对数几率回归代码详解

"""
训练样本（Training sample data.txt）
-0.017612	14.053064	0
-1.395634	4.662541	1
-0.752157	6.538620	0
-1.322371	7.152853	0
0.423363	11.054677	0
0.406704	7.067335	1
0.667394	12.741452	0
-2.460150	6.866805	1
0.569411	9.548755	0
-0.026632	10.427743	0
0.850433	6.920334	1
1.347183	13.175500	0
1.176813	3.167020	1
-1.781871	9.097953	0
-0.566606	5.749003	1
0.931635	1.589505	1
-0.024205	6.151823	1
-0.036453	2.690988	1
-0.196949	0.444165	1
1.014459	5.754399	1
1.985298	3.230619	1
-1.693453	-0.557540	1
-0.576525	11.778922	0
-0.346811	-1.678730	1
-2.124484	2.672471	1


测试样本（test data.txt）
0.074798	11.023650	0
-1.337472	0.468339	1
-0.102781	13.763651	0
-0.147324	2.874846	1
0.518389	9.887035	0
1.015399	7.571882	0
-1.658086	-0.027255	1
1.319944	2.171228	1
2.056216	5.019981	1
-0.851633	4.375691	1
-1.510047	6.061992	0
-1.076637	-3.181888	1
1.821096	10.283990	0
3.010150	8.401766	1
"""


import matplotlib.pyplot as plt
import  numpy as np

#函数功能：从文件中加载样本数据
def loadDataSet():
    dataMat = []  #创建数据列表
    labelMat = []  #创建标签列表
    fr = open('Training sample data.txt','r',encoding='UTF-8')  #打开文件
    for line in fr.readlines():  #逐行读取
        lineArr = line.strip().split()  #去回车，放入列表
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])  #向数据列表添加数据，其中float(lineArr[0])是一元线性方程中的x1，float(lineArr[1])是x2,1.0是为了把b添加到weights中，weights的第一个分量就是b
        labelMat.append(int(lineArr[2]))  #向标签列表添加标签
    fr.close()  #关闭文件
    return dataMat,labelMat #此时dataMat中的数据的下标是和labelMat中标签的下标相同的（一一对应的）


def sigmoid(Z):
    return 1.0 / (1 + np.exp(-Z))  #此处是e的矩阵次方


#梯度下降算法求weight向量
def gradient_descent(dataMat,labelMat):
    dataArr = np.mat(dataMat)  #将dataMat转换成矩阵形式，此处的mat()与array()作用相同
    labelArr = np.mat(labelMat).transpose()  #将labelMat转换成矩阵并转置（从行向量转置为列向量)
    m, n = np.shape(dataArr)  #矩阵的行（样本个数），列（属性个数）
    alpha = 0.001   #学习率（步长）
    maxCycles = 500 #迭代次数
    weights = np.ones((n,1))     ##初始化权重向量W（n*1的形状）,生成n个分量w，每个分量均为1，因为W的分量的要一一对应于样本的属性  https://blog.csdn.net/cunchi4221/article/details/107471968
    for i in range(maxCycles):
        h = sigmoid(dataArr * weights)  #预测值    和一元（多元）线性回归的不同之处就在于此处多套了一个sigmoid函数
        error = h - labelArr
        gradient = dataArr.transpose() * error   #梯度，此行加上两行为西瓜书p60式（3.30）
        weights = weights - alpha * gradient   #权值更新
    return weights


#绘制样本点图像和训练函数图像
def plotDataSet(weights):
    dataMat, labelMat = loadDataSet() #加载数据集
    dataArr = np.array(dataMat)  #将dataMat转换成numpy的array矩阵
    n = np.shape(dataArr)[0] #样本总个数,此时是100
    x1cord1 = []  #存放正样本的x1值,即横坐标值
    x2cord1 = []  #存放正样本的x2值，即纵坐标值
    """x1和x2是样本的两个属性特征，而样本的y值并非是labelMat向量，而是样本为正例的概率"""
    x1cord2 = []  #存放负样本的x1值
    x2cord2 = []  #存放负样本的x2值
    for i in range(n):  #根据数据集标签进行分类(有n个样本就有n个标签)
        if int(labelMat[i])  == 1:  #如果该标签是正样本
            x1cord1.append(dataArr[i, 1])  #将矩阵dataArr中i行1列的数存入xcord1列表中(x1值,即横坐标值)
            x2cord1.append(dataArr[i, 2])  #将矩阵dataArr中i行2列的数存入xcord1列表中（x2值，即纵坐标值）
        else:   #如果该标签是负样本
            x1cord2.append(dataArr[i, 1])  # 将矩阵dataArr中i行1列的数存入xcord2列表中(x1值,即横坐标值)
            x2cord2.append(dataArr[i, 2])  # 将矩阵dataArr中i行2列的数存入xcord2列表中（x2值，即纵坐标值）

    fig = plt.figure()
    ax = fig.add_subplot(111)   #添加subplot
    ax.scatter(x1cord1, x2cord1, s=20, c='red', marker='s', alpha=0.5)    #绘制正样本
    ax.scatter(x1cord2, x2cord2, s=20, c='green', alpha=0.5)  #绘制负样本

    #绘制训练函数图像
    x1 = np.arange(-3.0, 3.0, 0.1)  #x1是一个列表    arange()函数用法：https://blog.csdn.net/island1995/article/details/90179076
    x2 = (-weights[0,0] - weights[1,0] * x1) / weights[2,0]  #由w[1]*x1+w[2]*x2+w[0]=0推出的x2公式,注意weights是一个3*1的向量（矩阵），要取里面的数字要用[i,j]的方式取，不可直接用[i]的方式取数字！！！
    #print("x1:",x1)
    #print("x2:",x2)
    ax.plot(x1, x2)  #x1,x2维数要相同，同一维的分量构成一个坐标   plot()函数用法：https://blog.csdn.net/feng98ren/article/details/79392747

    plt.title('Training sample data')   #绘制title
    plt.xlabel('x1');plt.ylabel('x2')     #绘制坐标轴
    plt.show()  #显示


#测试测试样本
def test(weights):
    def loadTestDataSet():  #载入测试样本数据
        dataMatTest = []  # 创建数据列表
        fr = open('test data.txt', 'r', encoding='UTF-8')  # 打开文件
        for line in fr.readlines():  # 逐行读取
            lineArr = line.strip().split()  # 去回车，放入列表
            dataMatTest.append([1.0,float(lineArr[0]), float(lineArr[1])])  # 向数据列表添加数据，其中float(lineArr[0])是一元线性方程中的x1，float(lineArr[1])是x2,1.0是为了把b添加到weights中，weights的第一个分量就是b
        fr.close()  # 关闭文件
        dataArrTest = np.array(dataMatTest)
        return dataArrTest  # 此时dataMat中的数据的下标是和labelMat中标签的下标相同的（一一对应的）

    dataArrTest=loadTestDataSet()
    p1 = 1.0/(1 + np.exp(-dataArrTest * weights))  #西瓜书p59式（3.23）   测试样本为正例的概率
    p2 = 1 - p1     #测试样本为反例的概率

    labelMatTest = []  #测试样本的分类情况
    for p in p1 : #判断测试样本是正例还是反例
        if p > 0.5 :
            labelMatTest.append([1])  #要添加列表形式，否则后面无法与向量p1,p2合并
        else:
            labelMatTest.append([0])

    return p1, p2 ,np.array(labelMatTest)


if __name__ == '__main__' :
    # 计算weights
    dataMat, labelMat = loadDataSet()
    weights = gradient_descent(dataMat, labelMat)
    print("w为",weights[1:])
    print("b为",weights[0])

    # 绘制样本点和weights
    plotDataSet(weights)

    #测试
    p1, p2 ,labelMatTest= test(weights)
    p1_p2_labelMatTest = np.hstack((p1, p2,labelMatTest))  #合并3个向量
    print("[正例的概率 反例的概率 分类结果]：\n", p1_p2_labelMatTest)