西瓜书课后题—第三章（线性模型）

3.3 对率回归问题

直接贴代码：

import numpy as np
#西瓜数据集
X = np.mat([[0.697,0.460,1],[0.774,0.376,1],[0.634,0.264,1],[0.608,0.318,1],[0.556,0.215,1],
     [0.403,0.237,1],[0.481,0.149,1],[0.437,0.211,1],[0.666,0.091,1],[0.243,0.267,1],
     [0.245,0.057,1],[0.343,0.099,1],[0.639,0.161,1],[0.657,0.198,1],[0.360,0.370,1],
     [0.593,0.042,1],[0.719,0.103,1]])          #X是一个矩阵 mat
Y = np.array([1,1,1,1,1,1,1,1,0,0,0,0,0,0,0,0,0])     #Y是一个数组 array
beta = np.random.rand(3,1)      #随机生成待优化参数，返回的是一个二维数组array
err = 0.000001       #设定允许的错误率

def getP1(X,beta):     #得到p1,对应于书上介绍的分为第一类的概率
    m,_ = X.shape
    P1 = []      #列表类型
    for i in range(m):
        P1.append((np.e ** np.dot(X[i],beta)[0,0])/(1+np.e ** np.dot(X[i],beta)[0,0]))
        #因dot结果还是一个矩阵mat，所以要将其转成一个标量数值，用[0,0]
    return np.array(P1)      #将P1转成一个一维数组返回

def getDbeta(X,Y,beta):    #得到一阶导数
    P1 = getP1(X,beta)
    m,_ = X.shape
    Dbeta = np.zeros((3,1))   #返回的是二维array
    for i in range(m):
        Dbeta += X[i].T*(Y[i]-P1[i])  #mat和标量相乘
    return -Dbeta      #返回二维array类型

def getD2beta(X,beta):       #得到二阶导数
    P1 = getP1(X,beta)
    m,_  = X.shape
    D2beta = np.zeros((3,3))    #得到二维array
    for i in range(m):
        D2beta += np.dot(X[i].T,X[i])*P1[i]*(1-P1[i])    #mat和标量相乘，mat和二维数组相加
    return np.mat(D2beta)    #转成mat返回


def main(beta):
    Errbeta = np.ones((3,1))
    k = 0
    while np.linalg.norm(Errbeta)>err:   #此处用误差的二阶范数来作为衡量收敛与否的标准
        D = getDbeta(X,Y,beta)
        D2 = getD2beta(X,beta)
        beta_ = beta - np.dot(D2.I,D)   #注意，只有mat类型才可以求逆，也就是用 .I 操作
        Errbeta  = beta_ - beta
        beta = beta_
        k += 1
    print(beta,k)

if __name__ == __main__:
    main(beta)

得到结果如下：

[[ 3.15832966] [ 12.52119579] [ -4.42886451]] 分别对应着w1、w2、b

习题3.4：选择两个UCI数据集，比较10折交叉验证法和留一法所估计出的对率回归的错误率

请参见这篇博客，个人觉得写的非常清晰。

习题3.5：实现线性判别分析，并给出西瓜数据集上的实现。

直接给代码：

import numpy as np
import matplotlib.pyplot as plt

X0 = np.array([[0.697,0.460],[0.774,0.376],[0.634,0.264],[0.608,0.318],[0.556,0.215],
     [0.403,0.237],[0.481,0.149],[0.437,0.211]])      # 第一类数据   8 X 2
X1 = np.array([[0.666,0.091],[0.243,0.267],
     [0.245,0.057],[0.343,0.099],[0.639,0.161],[0.657,0.198],[0.360,0.370],
     [0.593,0.042],[0.719,0.103]])                    # 第二类数据    9 X 2
mean0 = np.mean(X0,axis=0,keepdims=True)              # 第一类数据的均值  1 X 2
mean1 = np.mean(X1,axis=0,keepdims=True)              # 第二类数据的均值   1 X 2

Sw = (X0-mean0).T.dot(X0-mean0) + (X1-mean1).T.dot(X1-mean1)  # 类内散度矩阵  2 X 2
omega = np.linalg.inv(Sw).dot((mean0-mean1).T)                # 待求得参数值 2 X 1

plt.plot(X0[:,0],X0[:,1],b*)      # 第一类数据点
plt.plot(X1[:,0],X1[:,1],r+)      # 第二类数据点
left_density = 0                    # 密度属性的取值范围
right_density = 1
left_sugerrate = 0                  # 含糖率属性的取值范围
right_sugerrate = -(right_density*omega[0])/omega[1]

此处之所以这样进行计算含糖率的另一个端点值，是因为我们默认求得的直线是通过原点的。
为什么可以认为通过原点？ 这是因为我们的目的是要让原有数据集投影到该直线上得到最大类间差距，最小 类内差距。所以，我们关心的仅仅是该直线的斜率，对于截距的大小并不关心，因为其并不会造成任何影响。

plt.plot([left_density,right_density],[left_sugerrate,right_sugerrate],g-)   #绘制该直线

plt.xlabel(density)
plt.ylabel(sugerrate)
plt.title(LDA)
plt.show()

最终的结果如下图所示：

备注：这里给出一个介绍纠错输出码的文章链接：