Logistic回归-从氙气病症预测病马的死亡率

摸鱼中qwq

于 2024-04-23 19:26:48 发布

阅读量1.2k

点赞数 23

文章标签：人工智能机器学习逻辑回归 python

本文链接：https://blog.csdn.net/weixin_62266832/article/details/138135975

版权

Logistic回归-从氙气病症预测病马的死亡率

机器学习课程作业，网上找了很多代码基本都是复制粘贴，原理也没有写的很清楚，这里根据个人理解整理一下。

原始数据链接：Horse Colic - UCI Machine Learning Repository

一、问题介绍

疝气病是描述马胃肠痛的术语，然而，这种病并不一定源自马的胃肠问题，其他问题也可能引发疝气病。数据一共包括368个样本和21个特征。该数据集中包含了医院检测马疝气病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。另外，除了部分指标主观和难以测量之外，该数据还存在一个问题，数据集中有30%的值是缺失的。

二、数据处理

将样本分为两组，其中300个作为训练用，记为horseColicTraining.txt，另外68个作为测试用，记为horseColicTest.txt。对数据集中缺失的数据，以0代替。处理后数据中每行为一个样本，前21列为数据，最后一列为0或1，代表分类的类别。

三、理论分析

该问题是一个二分类问题，对于测试集中每一组数据，我们需要输出0或1的预测值，因此采用logistic回归算法。

在logistic回归中，需要寻找参数 $\in R^{n}$ 和 $\in R$ ，对于给定特征 $\in R^{n}$ ，计算 $\sigma(x) = \frac{1}{1 + e^{- w^{T}x + b}}$ ，其中 $\sigma$ 为Sigmoid函数。若 $\sigma(x) > 0.5$ 则输出1，否则输出0。

为了方便起见，我们可以记新的 $\binom{x}{1}$ ， $\binom{w}{b}$ ，从而 $\sigma(x) = \frac{1}{1 + e^{- w^{T}x}}$

为了寻找合适的参数，可以采用极大似然估计法：

设训练数据集为 $\left\{ \left( {x^{1},y_{1}} \right),\left( {x^{2},y_{2}} \right),\ldots,\left( x^{k},y_{k} \right) \right\}$ 。其中每个 $x^{i} \in R^{n}$ ， $y_{i} \in \left\{ 0,1 \right\}$

由logistic回归模型，对于给定的 $x$ ， $P\left( {Y = 1} \middle| x \right) = \frac{e^{- w^{T}x}}{1 + e^{- w^{T}x}}$ ， $P\left( {Y = 0} \middle| x \right) = \frac{1}{1 + e^{- w^{T}x}}$

记 $P\left( {Y = 1} \middle| x \right)$

则似然函数为

${\prod\limits_{i = 1}^{n}{P\left( Y = y_{i} \middle| x^{i} \right)}} = {\prod\limits_{i = 1}^{n}{{f\left( x^{i} \right)}^{y_{i}}\left( 1 - f\left( x^{i} \right) \right)^{(1 - y_{i})}}}$

我们的目标是使似然函数最大，可以对似然函数取对数，利用 $y_{i} \in \left\{ 0,1 \right\}$ 得

${\sum\limits_{i = 1}^{n}\left( y_{i}{\ln{f\left( x^{i} \right)}} + \left( {1 - y_{i}} \right)ln\left( 1 - f\left( x^{i} \right) \right) \right.}$

于是只需要求解使 $l n L (w)$ 最大，这里可以采用梯度上升法(也可以利用梯度下降法求 $- l n L (w)$ 的最小值，原理相同)。

下面计算梯度。记 $z = w^{T}x$ ，注意到Sigmoid函数的性质：

$\sigma(z) = \frac{1}{1 + e^{- z}}$ ， ${\sigma(z)}^{'} = \sigma(z)\left( 1 - \sigma(z) \right)$

故

$\frac{\partial ln\sigma(z)}{\partial w_{i}} = \frac{\partial ln\sigma(z)}{\partial z}\frac{\partial z}{\partial w_{i}} = \frac{1}{\sigma(z)}\frac{d\sigma(z)}{dz}x_{i} = \frac{1}{\sigma(z)}\sigma(z)\left( {1 - \sigma(z)} \right)x_{i} = \left( {1 - \sigma(z)} \right)x_{i}$

类似可得 $\frac{\partial ln\left( 1 - \sigma(z) \right)}{\partial w_{i}} = - \sigma(z)x_{i}$

所以

$\frac{\partial lnL(w)}{\partial w_{j}} = {\sum\limits_{i = 1}^{n}\left( y_{i}\left( {1 - f\left( x^{i} \right)} \right)x_{j}^{i} - \left( {1 - y_{i}} \right)f\left( x^{i} \right)x_{j}^{i} \right)} = {\sum\limits_{i = 1}^{n}\left( y_{i} - f\left( x^{i} \right)x_{j}^{i} \right.}$

(注：这里的 $x_{j}^{i}$ 指第i个训练数据的第j个分量)

于是梯度上升算法可表示为 $\left. w_{j}\leftarrow w_{j} + \alpha{\sum\limits_{i = 1}^{n}\left( y_{i} - f\left( x^{i} \right)x_{j}^{i} \right.} \right.$ ，其中 $\alpha$ 为学习率。

具更体地，我们可以把训练数据中的 $x$ 写成矩阵形式：

$\left( \left( x^{1} \right),\ldots,\left( x^{k} \right) \right)^{T} = \begin{pmatrix} x_{1}^{1} & \cdots & x_{n}^{1} \\ \vdots & \ddots & \vdots \\ x_{1}^{k} & \cdots & x_{n}^{k} \end{pmatrix}$

把 $\left( y_{i} - f\left( x^{i} \right) \right.$ 记为 $e_{i}$ ， $\left( e_{1},\ldots,e_{k} \right)^{T}$

于是 ${\sum\limits_{i = 1}^{n}\left( y_{i} - f\left( x^{i} \right)x_{j}^{i} \right.} = {\sum\limits_{i = 1}^{n}e_{i}}x_{j}^{i}$

记 $\left( {w_{1},\ldots,w_{n}} \right)$ ，则梯度上升每次迭代可表示为 $\left. w\leftarrow w + \alpha x^{T}e \right.$

这便是代码中梯度上升的实现方式。

四、代码实现

Sigmoid()函数：每次处理一个numpy数组

def sigmoid(arr): # sigmoid阶跃函数  
 m, n = np.shape(arr) #m行1列  
 res = np.empty((m, 1))  
 for i in range(m):  
 inx = arr[i, 0]  
 if inx>=0: #对sigmoid函数的优化，避免了出现极大的数据溢出  
 res[i, 0] = (1.0/(1+np.exp(-inx)))  
 else:  
 res[i, 0] = (np.exp(inx)/(1+np.exp(inx)))  
 return res

gradAscent()函数：梯度上升算法实现

# dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。  
# classLabels 是类别标签，它是一个 1*300 的向量。为了便于矩阵计算，需要转置  
def gradAscent(dataMatIn, classLabels):  
 dataMatrix = np.mat(dataMatIn) # m*n矩阵，行数m=样本数，列数n=特征数  
 labelMat = np.mat(classLabels).transpose() # 行数m=样本数，列数=1  
 m,n = np.shape(dataMatrix)  
 alpha = 0.1 # alpha代表学习率  
 maxCycles = 5000 # 迭代次数  
 weights = np.ones((n, 1)) #n*1矩阵，代表回归系数，初始化为全1  
 for k in range(maxCycles):  
 # m*n 的矩阵 * n*1 的单位矩阵 ＝ m*1的矩阵  
 h = sigmoid(dataMatrix*weights) # 得到m*1矩阵，其中每行代表对每个样本计算sigmoid函数的值  
 error = (labelMat - h) # 均为m*1矩阵 每行为这个样本的误差  
 weights = weights + alpha * dataMatrix.transpose() * error # 矩阵乘法，最后得到回归系数  
 return weights

classifyVector()函数：分类

def classifyVector(inX, weights): # 分类函数，根据回归系数和特征向量来计算 Sigmoid的值  
 prob = sigmoid(sum(inX * weights))  
 if prob > 0.5:  
 return 1.0  
 else:  
 return 0.0

colicTest()函数：数据处理，统计，模型训练和测试

# 打开测试集和训练集,并对数据进行格式化处理  
def colicTest():  
 frTrain = open('horseColicTraining.txt')  
 frTest = open('horseColicTest.txt')  
 trainingSet = []  
 trainingLabels = []  
 # 解析训练数据集中的数据特征和Labels  
 # trainingSet 中存储训练数据集的特征(前22列)，trainingLabels 存储训练数据集的样本对应的分类标签(0或1)  
 for line in frTrain.readlines():  
 currLine = ["1"] + line.strip().split('\t')  
 lineArr = []  
 for i in range(22):  
 lineArr.append(float(currLine[i]))  
 trainingSet.append(lineArr)  
 trainingLabels.append(float(currLine[22]))  
 # 使用梯度上升算法 求得在此数据集上的最佳回归系数 trainWeights  
 trainWeights = gradAscent(np.array(trainingSet), trainingLabels)  
 errorCount = 0  
 numTestVec = 0.0  
 # 读取 测试数据集 进行测试，计算分类错误的样本条数和最终的错误率  
 for line in frTest.readlines():  
 numTestVec += 1.0  
 currLine = ["1"] + line.strip().split('\t')  
 lineArr = []  
 for i in range(22):  
 lineArr.append(float(currLine[i]))  
 if int(classifyVector(np.array(lineArr), trainWeights)) != int(currLine[22]):  
 errorCount += 1  
 errorRate = (float(errorCount) / numTestVec)  
 print("the error rate of this test is: %f" % errorRate)  
 return errorRate

五、测试结果

在测试中，选择学习率 $\alpha = 0.1$ ，迭代次数=500，进行训练。

在这里插入图片描述

测试结果显示，错误率为25.37%。考虑到训练数据有30%缺失，这个结果可以接受。

附录(完整代码)

import numpy as np


def sigmoid(arr):  # sigmoid阶跃函数
    m, n = np.shape(arr)  #m行1列
    res = np.empty((m, 1))
    for i in range(m):
        inx = arr[i, 0]
        if inx>=0:      #对sigmoid函数的优化，避免了出现极大的数据溢出
            res[i, 0] = (1.0/(1+np.exp(-inx)))
        else:
            res[i, 0] = (np.exp(inx)/(1+np.exp(inx)))
    return res


# dataMatIn 是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。
# classLabels 是类别标签，它是一个 1*300 的向量。为了便于矩阵计算，需要转置
def gradAscent(dataMatIn, classLabels):
    dataMatrix = np.mat(dataMatIn)  # m*n矩阵，行数m=样本数，列数n=特征数+1(22)
    labelMat = np.mat(classLabels).transpose()  # 行数m=样本数，列数=1
    m,n = np.shape(dataMatrix)
    alpha = 0.1  # alpha代表学习率
    maxCycles = 5000  # 迭代次数
    weights = np.ones((n, 1))  #n*1矩阵，代表回归系数，初始化为全1
    for k in range(maxCycles):
        # m*n 的矩阵 * n*1 的单位矩阵 ＝ m*1的矩阵
        h = sigmoid(dataMatrix*weights)     # 得到m*1矩阵，其中每行代表对每个样本计算sigmoid函数的值
        # labelMat是实际值
        error = (labelMat - h)              # 均为m*1矩阵 每行为这个样本的误差
        weights = weights + alpha * dataMatrix.transpose() * error  # 矩阵乘法，最后得到回归系数
    return weights



def classifyVector(inX, weights):  # 分类函数，根据回归系数和特征向量来计算 Sigmoid的值
    prob = sigmoid(sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0


# 打开测试集和训练集,并对数据进行格式化处理
def colicTest():
    frTrain = open('horseColicTraining.txt')
    frTest = open('horseColicTest.txt')
    trainingSet = []
    trainingLabels = []
    # 解析训练数据集中的数据特征和Labels
    # trainingSet 中存储训练数据集的特征(前22列)，trainingLabels 存储训练数据集的样本对应的分类标签(0或1)
    for line in frTrain.readlines():
        currLine = ["1"] + line.strip().split('\t')
        lineArr = []
        for i in range(22):
            lineArr.append(float(currLine[i]))
        trainingSet.append(lineArr)
        trainingLabels.append(float(currLine[22]))
    # 使用梯度上升算法 求得在此数据集上的最佳回归系数 trainWeights
    trainWeights = gradAscent(np.array(trainingSet), trainingLabels)
    errorCount = 0
    numTestVec = 0.0
    # 读取 测试数据集 进行测试，计算分类错误的样本条数和最终的错误率
    for line in frTest.readlines():
        numTestVec += 1.0
        currLine = ["1"] + line.strip().split('\t')
        lineArr = []
        for i in range(22):
            lineArr.append(float(currLine[i]))
        if int(classifyVector(np.array(lineArr), trainWeights)) != int(currLine[22]):
            errorCount += 1
    errorRate = (float(errorCount) / numTestVec)
    print("the error rate of this test is: %f" % errorRate)
    return errorRate


# 调用 colicTest() 1次并求结果的平均值(确定性算法无需重复)
def multiTest():
    numTests = 1
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" % (numTests, errorSum / float(numTests)))


multiTest()
均值(确定性算法无需重复)
def multiTest():
    numTests = 1
    errorSum = 0.0
    for k in range(numTests):
        errorSum += colicTest()
    print("after %d iterations the average error rate is: %f" % (numTests, errorSum / float(numTests)))


multiTest()

摸鱼中qwq

关注

23
点赞
踩
11

收藏

觉得还不错? 一键收藏
1
评论
Logistic回归-从氙气病症预测病马的死亡率

疝气病是描述马胃肠痛的术语，然而，这种病并不一定源自马的胃肠问题，其他问题也可能引发疝气病。该数据集中包含了医院检测马疝气病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。另外，除了部分指标主观和难以测量之外，该数据还存在一个问题，数据集中有30%的值是缺失的。将样本分为两组，其中300个作为训练用，记为horseColicTraining.txt，另外68个作为测试用，记为horseColicTest.txt。对数据集中缺失的数据，以0代替。，迭代次数=500，进行训练。
复制链接

扫一扫