北太天元应用案例分享：用向量机做心脏病筛查

最新推荐文章于 2024-07-26 18:11:15 发布

北太天元

最新推荐文章于 2024-07-26 18:11:15 发布

阅读量1k

点赞数 27

文章标签： matlab 算法开发语言

本文链接：https://blog.csdn.net/BTZH2022/article/details/140472947

版权

数据

age,sex,cp,trestbps,chol,fbs,restecg,thalach,exang,oldpeak,slope,ca,thal,target

42,1,1,120,295,0,1,162,0,0,2,0,2,1

48,0,2,130,275,0,1,139,0,0.2,2,0,2,1

54,1,2,120,258,0,0,147,0,0.4,1,0,3,1

54,0,2,108,267,0,0,167,0,0,2,0,2,1

51,0,2,130,256,0,0,149,0,0.5,2,0,2,1

57,0,0,128,303,0,0,159,0,0,2,1,2,1

61,1,0,138,166,0,0,125,1,3.6,1,1,2,0

45,1,0,115,260,0,0,185,0,0,2,0,2,1

60,0,3,150,240,0,1,171,0,0.9,2,0,2,1

49,0,0,130,269,0,1,163,0,0,2,0,2,1

53,1,2,130,246,1,0,173,0,0,2,3,2,1

54,1,0,140,239,0,1,160,0,1.2,2,0,2,1

51,1,0,140,298,0,1,122,1,4.2,1,3,3,0

49,1,2,118,149,0,0,126,0,0.8,2,3,2,0

58,1,2,105,240,0,0,154,1,0.6,1,0,3,1

48,1,1,130,245,0,0,180,0,0.2,1,0,2,1

48,1,0,122,222,0,0,186,0,0,2,0,2,1

52,1,2,138,223,0,1,169,0,0,2,4,2,1

54,1,1,192,283,0,0,195,0,0,2,1,3,0

63,0,0,150,407,0,0,154,0,4,1,3,3,0

70,1,1,156,245,0,0,143,0,0,2,0,2,1

这是一个关于心脏病筛查的数据集片段，每一行代表一个病人的信息。以下是各列字段的解释：

age：年龄 - 病人的年龄（以岁为单位）。
sex：性别 - 通常为0（女性）或1（男性）。
cp：胸痛类型 - 这是一个分类变量，表示胸痛的类型。不同的数字可能代表不同的胸痛类型（例如，0可能表示无胸痛，1、2、3、4可能代表不同类型的胸痛）。
trestbps：静息血压 - 病人在静息状态下的血压（单位可能是mmHg）。
chol：血清胆固醇 - 病人的血清胆固醇水平（单位可能是mg/dL）。
fbs：空腹血糖 - 病人的空腹血糖水平（0可能表示正常，1可能表示高于正常）。
restecg：静息心电图结果 - 这可能是一个分类变量，表示静息状态下的心电图结果。
thalach：最大心率 - 病人能达到的最大心率（单位可能是次/分钟）。
exang：运动诱发的心绞痛 - 0可能表示没有，1可能表示有。
oldpeak：ST段下降 - 这是心电图的一个参数，可能表示心肌缺血的程度（单位可能是mV或其他相关单位）。
slope：ST段峰值斜率 - 这也是心电图的一个参数，可能与心脏的病理状况有关（例如，1、2、3可能表示不同的斜率类型）。
ca：冠状动脉造影的主要血管数量 - 这可能是一个分类变量，表示在冠状动脉造影中检测到的主要血管数量（0可能表示无血管受损，其他数字可能表示受损血管的数量）。
thal：缺陷类型 - 这可能是一个分类变量，表示某种心脏缺陷或病理状况的类型。
target：目标变量 - 这是我们要预测或分类的变量。通常，1可能表示病人有心脏病，0可能表示没有。

下面是北太天元软件的代码, 关于心脏病数据集的多类支持向量机（SVM）分类器的完整实现，包括数据预处理、主成分分析（PCA）降维、交叉验证划分、模型训练和测试、以及结果可视化。下面是对代码的详细解读：

导入插件和数据

北太天元代码

load_plugin("optimization");
data = readmatrix('heart disease.csv');

load_plugin("optimization"): 加载优化工具箱，用于求解二次规划问题。

readmatrix('heart disease.csv'): 读取心脏病数据集。

数据预处理

北太天元代码

if any(any(ismissing(data)))
    data = rmmissing(data);
end

检查数据中是否有缺失值，并使用rmmissing函数删除包含缺失值的行。

分离特征和标签

北太天元代码

X = data(:,1:end-1);
y = data(:,end);

将数据集中的最后一列作为标签，其余列作为特征。

主成分分析降维

北太天元代码

[coeff,score,latent] = pca(X);
X_pca = score(:,1:2);

使用自定义的pca函数对特征进行降维，保留前两个主成分。

划分训练集和测试集

北太天元代码

cv = my_cvpartition(size(X_pca,1),'HoldOut',0.3);
% ...（省略部分代码）

使用自定义的my_cvpartition函数进行HoldOut交叉验证划分，将数据集分为70%的训练集和30%的测试集。

训练SVM模型

北太天元代码

SVMModel = my_fitcecoc(X_train,Y_train,C);

使用自定义的my_fitcecoc函数训练多类SVM模型。该函数内部采用一对一（One-vs-One）策略，为每个类别对训练一个SVM分类器。

预测测试集并计算准确率

北太天元代码

[label,score] = predict(SVMModel,X_test);
accuracy = sum(label-1 == Y_test) / length(Y_test);

使用predict函数对测试集进行预测，并计算准确率。注意这里假设标签是从1开始的，所以计算准确率时需要将预测的标签减1。

可视化

matlab复制代码

gscatter(X_pca(:,1),X_pca(:,2),y);
% ...（省略绘制决策边界的代码）

使用gscatter函数绘制原始数据点的散点图，并使用轮廓线绘制决策边界。

自定义函数解读

pca: 实现主成分分析算法，对数据进行降维。

my_cvpartition: 自定义的交叉验证划分函数，这里只实现了HoldOut方法。

my_fitcecoc: 实现多类SVM分类器的训练，采用一对一策略。

my_fitcsvm和my_fitcsvm_soft: 实现硬间隔和软间隔SVM分类器的训练，使用二次规划求解。

predict: 对测试集进行预测，并返回预测的标签和得分。

predictOneVsAll: 使用单个SVM模型进行预测，并返回预测结果和决策函数的值。

注意事项

代码中存在一些假设和简化，例如标签从1开始、直接计算决策函数的值作为得分等，这些在实际应用中可能需要根据具体情况进行调整。
在使用二次规划求解SVM问题时，需要注意目标函数和约束条件的设置，以确保求解的正确性。
可视化部分仅绘制了前两个主成分的散点图和决策边界，对于更高维度的数据可能需要进行其他形式的可视化。

下面是全部代码

load_plugin("optimization");
% 导入数据
data = readmatrix('heart disease.csv');
% 检查并处理缺失值
if any(any(ismissing(data)))
    data = rmmissing(data);
end
% 分离特征和标签
X = data(:,1:end-1);
y = data(:,end);
% 主成分分析降维
[coeff,score,latent] = pca(X);
X_pca = score(:,1:2);
% 划分训练集和测试集
cv = my_cvpartition(size(X_pca,1),'HoldOut',0.3);
idx = cv.test;
train_idx = setdiff(1: size(X_pca,1), idx);
X_train = X_pca(train_idx,:);
Y_train = y(train_idx,:);
X_test = X_pca(idx,:);
Y_test = y(idx,:);
% 训练SVM模型
C = 1e2; % 设置较大的C值以确保硬间隔分类（对于线性可分数据）
SVMModel = my_fitcecoc(X_train,Y_train,C);
% 预测测试集
[label,score] = predict(SVMModel,X_test);
% 计算准确率
accuracy = sum(label-1 == Y_test) / length(Y_test);
% 可视化
gscatter(X_pca(:,1),X_pca(:,2),y);
hold on;
% 绘制决策边界c
d = 0.15;
[x1Grid,x2Grid] = meshgrid(min(X_pca(:,1)):d:max(X_pca(:,1)),...
min(X_pca(:,2)):d:max(X_pca(:,2)));
xGrid = [x1Grid(:),x2Grid(:)];
[~,scores] = predict(SVMModel,xGrid);
contour(x1Grid,x2Grid,reshape(scores(:,2),size(x1Grid)),[0 0],'k');
% 输出准确率
fprintf('The accuracy of the SVM model is %.2f%%\n', accuracy * 100);
hold off;
function [coeff, score, latent] = pca(X)
    % X: 数据矩阵，每一列是一个特征，每一行是一个样本
    % coeff: 主成分系数
    % score: 表示主成分得分
    % latent: 主成分对应的特征值
    % 标准化数据（均值为0，方差为1）
    X = (X - mean(X)) ./ std(X);
%    X = (X - mean(X));
    % 计算协方差矩阵
    CovMat = cov(X);
    % 对协方差矩阵进行特征值分解
    [V, D] = eig(CovMat);
    % 将特征值按降序排序，并获取对应的特征向量
    [latent, order] = sort(diag(D), 'descend');
    coeff = V(:, order);
    % 计算主成分得分
    score = X * coeff;
end
function cv = my_cvpartition(n, method, param)
    % n: 总样本数
    % method: 分区方法，目前仅支持 'HoldOut'
    % param: 分区方法的参数，对于 'HoldOut'，该参数为留出的比例
    % 初始化输出结构体
    cv = struct('train', [], 'test', []);
    if strcmp(method, 'HoldOut')
        % 生成一个从 1 到 n 的整数数组
        idx = 1:n;
        % 随机打乱 idx 数组
        idx = idx(randperm(n));
        % 计算留出样本的数量
        numHoldOut = floor(n * param);
        % 分配训练和测试索引
        cv.test = idx(1:numHoldOut);
        cv.train = idx(numHoldOut+1:end);
    else
        error('Unsupported partition method.');
    end
end
function SVMModel = my_fitcecoc(X_train, Y_train, C)
    % X_train: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
    % Y_train: 标签向量 (n x 1)
    % C: 正则化参数，控制对错分样本的惩罚程度
    % 获取类别数量
    uniqueClasses = unique(Y_train);
    numClasses = length(uniqueClasses);
    % 初始化 SVM 模型结构体数组
    SVMModel = struct('Classifiers', cell(numClasses*(numClasses-1)/2, 1), ...
    'ClassPairs', [], 'ClassPairsIndex', [], 'UniqueClasses', uniqueClasses);
    % 训练一对一 SVM 分类器
    classifierIndex = 1;
    for i = 1:numClasses
        for j = (i+1):numClasses
            % 提取当前类别对的训练数据
            classI = Y_train == uniqueClasses(i);
            classJ = Y_train == uniqueClasses(j);
            X_train_ij = [X_train(classI, :); X_train(classJ, :)];
            Y_train_ij = [ones(sum(classI), 1); -1*ones(sum(classJ), 1)];
            % 训练 SVM 分类器
            SVMModel.Classifiers{classifierIndex} = my_fitcsvm_soft(X_train_ij, Y_train_ij, C);
            % 记录当前分类器对应的类别对
            SVMModel.ClassPairs(classifierIndex, :) = [uniqueClasses(i), uniqueClasses(j)];
            SVMModel.ClassPairsIndex(classifierIndex, :) = [i,j];
            classifierIndex = classifierIndex + 1;
        end
    end
end
function [wb] = my_fitcsvm(X, Y)
    % X: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
    % Y: 标签向量 (n x 1)，取值为 +1 或 -1
    load_plugin("optimization");
    % 假设X和Y已经定义，如之前的示例
    N = size(X, 1); % 数据点的数量
    D = size(X,2); % 数据的维度
    % 将w和b组合成一个向量，以便使用quadprog
    % 注意：这里我们将w放在前面，b放在最后
    p = rand(D + 1, 1); % 初始猜测解（通常设为0）
    Aeq = []; % 没有等式约束
    beq = [];
    % 构造二次规划的目标函数和线性不等式约束
    H = eye(D+1); % Hessian矩阵（目标函数的二次项系数）
    f = zeros(D+1,1); % 目标函数的一次项系数（对于SVM原始问题，通常为负）
    % 线性不等式约束 Ax <= b
    % 对于每个数据点 (x_i, y_i)，我们有 y_i * (w' * x_i + b) >= 1
    A = [-Y.*X, -Y]; % 不等式约束的系数矩阵
    b = -ones(N, 1); % 不等式约束的右侧向量
    % 使用quadprog求解二次规划问题
    options = optimoptions('quadprog','Algorithm','interior-point');
    [w_b, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, [], [], p, options);
    % 分离出w和b
    w = w_b(1:D);
    b = w_b(D+1);
    wb = struct('w',w,'b',b);
end
function [wb] = my_fitcsvm_soft(X, Y, C)
    % X: 特征矩阵 (n x d)，其中 n 是样本数，d 是特征维度
    % Y: 标签向量 (n x 1)，取值为 +1 或 -1
    % C: 正则化参数，控制对错分样本的惩罚程度
    % 假设X和Y已经定义，如之前的示例
    N = size(X, 1); % 数据点的数量
    D = size(X, 2); % 数据的维度
    % 初始化松弛变量（slack variables）
    xi = zeros(N, 1);
    % 将w和b以及松弛变量组合成一个向量，以便使用quadprog
    p = zeros(D + 1 + N, 1); % 初始猜测解
    % 构造二次规划的目标函数和线性不等式约束
    H = [diag( [ones(1,D),0] ), zeros(D+1, N); zeros(N, D+1), zeros(N,N)]; % Hessian矩阵
    f = [zeros(D+1, 1); C*ones(N, 1)]; % 目标函数的一次项系数
    % 线性不等式约束 Ax <= b
    % 对于每个数据点 (x_i, y_i)，我们有 y_i * (w' * x_i + b) >= 1 - xi_i
    A = [-Y.*X, -Y, diag( -ones(1,N) )]; % 不等式约束的系数矩阵
    b = -ones(N, 1); % 不等式约束的右侧向量
    Aeq = []; % 没有等式约束
    beq = [];
    lb = [ -inf*ones(D+1,1); zeros(N,1)];
    ub = inf*ones(D+1+N,1);
    % 使用quadprog求解二次规划问题
    options = optimoptions('quadprog','Algorithm','interior-point');
    [w_b_xi, fval, exitflag, output] = quadprog(H, f, A, b, Aeq, beq, lb, ub, p, options)
    % 分离出w, b和xi
    w = w_b_xi(1:D);
    b = w_b_xi(D+1);
    xi = w_b_xi(D+2:end);
    % 构造并返回结构体wb，包含w和b
    wb = struct('w', w, 'b', b);
end
function [label, score] = predict(SVMModel, X_test)
    % SVMModel: 训练好的多类SVM模型，由my_fitcecoc函数返回
    % X_test: 测试集特征数据
    % label: 预测的类别标签
    % score: 预测的得分（可选，这里简化为投票得分）
    % 初始化预测标签和得分
    [numTest, ~] = size(X_test);
    numClasses = length(SVMModel.Classifiers);
    label = zeros(numTest, 1);
    score = zeros(numTest, numClasses);
    % 对每个测试样本进行预测
    for i = 1:numTest
        sample = X_test(i, :);
        votes = zeros(1, numClasses+1);
        % 使用每个SVM分类器进行预测，并进行投票
        for j = 1:numClasses
            svmModel = SVMModel.Classifiers{j};
            [~, ~, ~, output] = predictOneVsAll(svmModel, sample);
            if output == 1
                classIdx1 = SVMModel.ClassPairsIndex(j, 1);
                votes(classIdx1) = votes(classIdx1) + 1;
            else
                classIdx2 = SVMModel.ClassPairsIndex(j, 2);
                votes(classIdx2) = votes(classIdx2) + 1;
            end
        end
        % 找到得票最多的类别作为预测标签
        [~, maxVoteIdx] = max(votes);
        label(i) = SVMModel.UniqueClasses(maxVoteIdx);
        score(i, maxVoteIdx) = votes(maxVoteIdx); % 将最高票数作为得分
    end
end
function [predictedLabel, predictedScore, decisionValues, output] = predictOneVsAll(svmModel, sample)
    % 使用单个SVM模型进行预测
    % 这里简化为直接计算决策函数的值，实际应用中可能需要更复杂的处理
    w = svmModel.w;
    b = svmModel.b;
    decisionValues = dot(w, sample) + b;
if decisionValues >= 0
        predictedLabel = 1; % 属于正类
    else
        predictedLabel = -1; % 属于负类
    end
    predictedScore = abs(decisionValues); % 决策函数的绝对值作为得分
    output = predictedLabel; % 用于投票的输出
end

处理非连续或非标准化标签

处理非连续或非标准化标签（例如，不是从1开始的整数或不是连续整数）是机器学习中的一个常见问题。在我的SVM分类器实现中，为了处理这种情况，我引入了SVMModel.UniqueClasses和SVMModel.ClassPairsIndex两个字段。下面我将详细解释这两个字段的作用以及如何在预测过程中使用它们。

SVMModel.UniqueClasses

SVMModel.UniqueClasses存储了数据集中所有唯一的类别标签。当训练数据中的标签不是从1开始的连续整数时，这个字段特别有用。例如，如果标签是[-1, 0, 2]，那么UniqueClasses就会是[-1, 0, 2]。

在训练过程中，这个字段用于记录所有不同的标签，以便在预测时能够正确地将预测的类别索引映射回原始标签。

SVMModel.ClassPairsIndex

SVMModel.ClassPairsIndex存储了训练一对一SVM分类器时所使用的类别对的索引。由于我们可能不能直接使用原始标签作为索引（尤其是当它们不是从1开始的连续整数时），我们需要一个从1开始的连续整数索引来标识不同的类别。这个索引在训练SVM分类器和投票过程中都是必要的。

例如，如果UniqueClasses是[-1, 0, 2]，我们可以为这些类别分配一个从1开始的索引：-1对应索引1，0对应索引2，2对应索引3。然后，ClassPairsIndex就会存储这些索引对应的类别对，例如[1, 2]、[1, 3]和[2, 3]。

在预测过程中使用这些字段

在预测过程中，SVMModel.UniqueClasses和SVMModel.ClassPairsIndex一起使用来确保正确的类别映射。以下是一个简化的流程：

对于测试集中的每个样本，使用所有的一对一SVM分类器进行预测。
根据每个分类器的预测结果（正类或负类）进行投票。但是，由于我们使用了从1开始的连续整数索引来标识类别，因此投票结果也是基于这些索引的。
找到得票最多的索引（例如，索引2）。
使用SVMModel.ClassPairsIndex（如果必要的话）和SVMModel.UniqueClasses将得票最多的索引转换回原始标签。在这个例子中，索引2对应UniqueClasses中的第二个元素，即标签0。
返回预测的原始标签作为预测结果。

通过这种方式，即使训练数据中的标签不是从1开始的连续整数，我的SVM分类器实现也能够正确地预测类别并返回原始标签。

软间隔的SVM（支持向量机）算法

软间隔的SVM（支持向量机）算法是为了解决线性不可分或存在噪声数据的问题而引入的。在硬间隔SVM中，所有样本都必须被正确地分类，这在实际应用中往往很难满足。软间隔SVM则允许部分样本被错误分类，通过在目标函数中引入一个惩罚项（或称为损失函数）来控制这种错误分类的程度。

软间隔SVM的算法可以通过二次规划问题来描述。二次规划问题是一个优化问题，它涉及最小化或最大化一个二次函数，同时满足一组线性不等式约束。在软间隔SVM中，这个二次函数通常与超平面的参数（如权重向量w和偏置项b）有关，而约束条件则与样本点到超平面的距离和分类正确性有关。

下面是软间隔SVM二次规划问题的标准形式：

灵感来自卢朓老师《数值方法：原理、算法及应用》课堂学生大作业~

学生作业的原代码是用MATLAB写的，以上为卢朓老师修改之后可以在北太天元上运行的代码分享

原视频可查看：机器学习预测心脏病风险-北大王暄_哔哩哔哩_bilibili

北太天元

关注

27
点赞
踩
23

收藏

觉得还不错? 一键收藏
0
评论
北太天元应用案例分享：用向量机做心脏病筛查

pca: 实现主成分分析算法，对数据进行降维。my_cvpartition: 自定义的交叉验证划分函数，这里只实现了HoldOut方法。my_fitcecoc: 实现多类SVM分类器的训练，采用一对一策略。my_fitcsvm和my_fitcsvm_soft: 实现硬间隔和软间隔SVM分类器的训练，使用二次规划求解。predict: 对测试集进行预测，并返回预测的标签和得分。predictOneVsAll: 使用单个SVM模型进行预测，并返回预测结果和决策函数的值。
复制链接

扫一扫