Описание данных

Этот пример использует данные фонокардиограммы (PCG), полученные от людей с нормальной и аномальной функцией сердца. Набор данных состоит из 3829 записей, 2575 от лиц с нормальной функцией сердца и 1254 записей от лиц с аномальной функцией сердца. Каждая запись имеет длину 10 000 выборок и дискретизируется со скоростью 2 кГц. Это представляет пять секунд данных фонокардиограммы. Набор данных построен из данных обучения и валидации, используемых в Вычисление PhysioNet in Cardiology Challenge 2016 [1] [2].

Загрузка данных

Первый шаг - загрузка данных из репозитория GitHub. Чтобы загрузить данные, нажмите Code и выберите Download ZIP. Сохраните файл physionet_phonocardiogram-main.zip в папке, в которой у вас есть разрешение на запись. Инструкции для этого примера предполагают, что вы загрузили файл во временную директорию, (tempdir в MATLAB™). Измените последующие инструкции для распаковки и загрузки данных, если вы решите загрузить данные в папку, отличную от tempdir.

Файл physionet_phonocardiogram-main.zip содержит

и PCG_Data.zip содержит

heartSoundData.mat содержит данные и метки классов, используемые в этом примере. Файл .txt, Modified_physionet_data.txt, требуется политикой копирования PhysioNet и предоставляет атрибуты источника для данных, а также описание того, как каждый сигнал в heartSoundData.mat соответствует файлу в исходных данных PhysioNet. extrafiles.mat также содержит атрибуты исходного файла и поясняется в Modified_physionet_data.txt файле. Единственный файл, требуемый для запуска примера heartSoundData.mat.

Загрузка данных

Если вы следовали инструкциям по загрузке в предыдущем разделе, введите следующие команды, чтобы разархивировать два файла архива.

unzip(fullfile(tempdir,'physionet_phonocardiogram-main.zip'),tempdir)
unzip(fullfile(tempdir,'physionet_phonocardiogram-main','PCG_Data.zip'),...
    fullfile(tempdir,'PCG_Data'))

После того, как вы разархивируете файл PCG_Data.zip, загружайте данные в MATLAB.

load(fullfile(tempdir,'PCG_Data','heartSoundData.mat'))

heartSoundData массив структур с двумя полями: Data и Classes. Data является матрицей 10000 на 3829, где каждый столбец является записью PCG. Classes является 3829 на 1 категориальным массивом диагностических меток, по одному для каждого столбца Data. Поскольку это двоичная задача классификации, классы являются «нормальными» и «ненормальными». Как указывалось ранее, существует 2575 обычных записей и 1254 аномальных записей. Эквивалентно, 67,25% примеров в данных - от людей с нормальной функцией сердца, в то время как 32,75% - от людей с аномальной функцией сердца. Проверить это можно путем ввода:

summary(heartSoundData.Classes)

     normal        2575 
     abnormal      1254 

countcats(heartSoundData.Classes)./sum(countcats(heartSoundData.Classes))

ans = 2×1

    0.6725
    0.3275

Вейвлет

Использование waveletScattering для создания вейвлета рассеяния. Установите инвариантную шкалу так, чтобы она совпадала с длиной сигнала. Сеть рассеяния по умолчанию имеет два вейвлета (банки фильтров). Первый банк вейвлет-фильтров имеет восемь вейвлет на октаву. Вторая группа фильтров имеет по одному вейвлет на октаву. Установите 'OptimizePath' свойство к true.

N = 1e4;
sn = waveletScattering('SignalLength',N,'InvarianceScale',N,'OptimizePath',true);

Создайте наборы для обучения и тестирования

Функция помощника, partition_heartsoundsразделяет 3829 наблюдений так, что 70% (2680) находятся в наборе обучающих данных с 1802 нормальными и 878 ненормальными. Остальные 1149 записей (773 нормальных и 376 ненормальных) сохраняются в тестовом наборе для предсказания. Генератор случайных чисел посеян внутри вспомогательной функции, поэтому результаты повторяются. Код для partition_heartsounds и все другие вспомогательные функции, используемые в этом примере, приведены в разделе Вспомогательные функции в конце примера.

[trainData, testData, trainLabels, testLabels] = ...
    partition_heartsounds(70,heartSoundData.Data,heartSoundData.Classes);

Номера каждого класса можно проверить в обучающих и тестовых наборах.

summary(trainLabels)

     normal        1802 
     abnormal       878 

summary(testLabels)

     normal        773 
     abnormal      376 

Обратите внимание, что обучающие и тестовые наборы были разбиты так, что доля «нормальных» и «ненормальных» записей в обучающих и тестовых наборах совпадает с их пропорциями в общих данных. Подтвердить это можно следующими способами.

countcats(trainLabels)./sum(countcats(trainLabels))

ans = 2×1

    0.6724
    0.3276

countcats(testLabels)./sum(countcats(testLabels))

ans = 2×1

    0.6728
    0.3272

Функции рассеяния

Получите преобразование рассеяния всех 2680 записей в набор обучающих данных. Для многомерных временных рядов преобразование рассеяния предполагает, что каждый столбец является отдельным сигналом. Используйте 'log' опция для получения естественного логарифма коэффициентов рассеяния.

scat_features_train = featureMatrix(sn,trainData,'Transform','log');

Для заданных параметров рассеяния scat_features_train является матрицей 279 на 5 на 2680. Для каждого из 2680 сигналов существует 279 путей рассеяния и пять окон рассеяния. В порядок, чтобы передать это в классификатор SVM, измените форму тензора в матрицу 13400 на 279, где каждая строка представляет одно окно рассеяния через 279 путей рассеяния. Общее количество строк равно продукту 5 и 2680 (количество записей в обучающих данных).

Nseq = size(scat_features_train,2);
scat_features_train = permute(scat_features_train,[2 3 1]);
scat_features_train = reshape(scat_features_train,...
    size(scat_features_train,1)*size(scat_features_train,2),[]);

Повторите процесс для тестовых данных.

scat_features_test = featureMatrix(sn,testData,'Transform','log');
scat_features_test = permute(scat_features_test,[2 3 1]);
scat_features_test = reshape(scat_features_test,...
    size(scat_features_test,1)*size(scat_features_test,2),[]);

Здесь мы тиражируем метки так, чтобы у нас была метка для каждого временного окна рассеяния.

[sequence_labels_train,sequence_labels_test] = ...
    createSequenceLabels_heartsounds(Nseq,trainLabels,testLabels);

Подбор SVM к обучающим данным. В этом примере мы используем кубическое полиномиальное ядро. После подбора кривой SVM к обучающим данным мы выполняем 5-кратную перекрестную валидацию, чтобы оценить ошибку обобщения обучающих данных. Здесь каждое окно рассеяния классифицируется отдельно.

rng default;
classificationSVM = fitcsvm(...
    scat_features_train, ...
    sequence_labels_train , ...
    'KernelFunction', 'polynomial', ...
    'PolynomialOrder', 3, ...
    'KernelScale', 'auto', ...
    'BoxConstraint', 1, ...
    'Standardize', true, ...
    'ClassNames', categorical({'normal','abnormal'}));
kfoldmodel = crossval(classificationSVM, 'KFold', 5);

Вычислите потерю в процентах и отобразите матрицу неточностей.

predLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100;
fprintf('Loss is %2.2f percent\n',loss);

Loss is 0.96 percent

accuracy = 100-loss;
fprintf('Accuracy is %2.2f percent\n',accuracy);

Accuracy is 99.04 percent

confmatCV = confusionchart(sequence_labels_train,predLabels);

Обратите внимание, что сеть рассеяния приводит к приблизительно 99-процентной точности, когда каждое временное окно классифицируется отдельно. Тем не менее, эффективность на самом деле лучше этого значения, потому что у нас есть пять окон рассеяния на запись, и 99-процентная точность основана на классификации всех окон отдельно. В этом случае используйте большинство голосов для получения одного назначения класса на запись. Классовое голосование соответствует режиму голосования за пять окон. Если уникальный режим не найден, helperMajorityVote классифицирует этот набор окон рассеяния как 'NoUniqueMode' для указания ошибки классификации. Это приводит к появлению дополнительного столбца в матрице неточностей.

classes = categorical({'abnormal','normal'});
ClassVotes = helperMajorityVote(predLabels,trainLabels,classes);
CVaccuracy = sum(eq(ClassVotes,trainLabels))./numel(trainLabels)*100;
fprintf('The true cross-validation accuracy is %2.2f percent.\n',CVaccuracy);

The true cross-validation accuracy is 99.89 percent.

Отобразите матрицу неточностей для классификаций голосов большинства.

cmCV = confusionchart(trainLabels,ClassVotes);

Точность перекрестной проверки обучающих данных на самом деле составляет 99,89 процента. Существует две нормальные записи, которые неправильно классифицируются как ненормальные. Одна ненормальная запись классифицируется как нормальная.

Используйте модель SVM подгонка с обучающими данными, чтобы делать предсказания классов на удерживаемых тестовых данных.

predTestLabels = predict(classificationSVM,scat_features_test);

Определите точность предсказаний на тестовом наборе с помощью большинства голосов.

ClassVotes = helperMajorityVote(predTestLabels,testLabels,classes);
testaccuracy = sum(eq(ClassVotes,testLabels))./numel(testLabels)*100;
fprintf('The test accuracy is %2.2f percent.\n',testaccuracy);

The test accuracy is 91.82 percent.

cmTest = confusionchart(testLabels,ClassVotes);

Из 1149 тестовых записей примерно 92% правильно классифицируются как «нормальные» или «ненормальные». Из 773 обычных записей PCG в тестовом наборе 732 правильно классифицированы. Из 376 аномальных записей в тестовом наборе 323 правильно классифицированы.

Точность, отзыв и F1 счет

В задаче классификации точностью для класса является количество правильных положительных результатов, разделенных на количество положительных результатов. Другими словами, из всех записей, которые классификатор присваивает заданной метке, какая пропорция на самом деле принадлежит классу. Вызов определяется как количество правильных меток, разделенных на количество меток для данного класса. В частности, из всех записей, принадлежащих классу, какая пропорция сделала нашу метку классификатора этим классом. Оценивая точность вашего классификатора, вы идеально хотите сделать хорошо как на точности, так и на отзыве. Например, предположим, что у нас был классификатор, который пометил каждую запись PCG как ненормальную. Тогда наш отзыв для аномального класса был бы на 100%. Все записи, принадлежащие ненормальному классу, будут помечены как ненормальные. Однако точность будет низкой. Поскольку наш классификатор пометил все записи как ненормальные, в этом случае будет 2575 ложных срабатываний для точности 1254/3829, или 32,75%. F1 счета является гармоническим средним значением точности и отзыва и предоставляет одну метрику, которая суммирует эффективность классификатора с точки зрения как отзыва, так и точности. Функция помощника, helperF1heartSounds, вычисляет точность, отзыв и F1 счетов для результатов классификации на наборе тестов и возвращает эти результаты в таблицу.

PRTable = helperF1heartSounds(cmTest.NormalizedValues);
disp(PRTable)

                Precision    Recall    F1_Score
                _________    ______    ________
    Abnormal     88.736      85.904     87.297 
    Normal       93.248      94.696     93.967 

В этом случае F1 счетов для ненормальных и нормальных групп подтверждают, что наша модель имеет и хорошую точность, и отзыв. В двоичной классификации просто определить как точность, так и отзыв непосредственно из матрицы неточностей. Чтобы увидеть это, снова постройте график матрицы неточностей для удобства.

cmTest = confusionchart(testLabels,ClassVotes);

Отзывом для ненормального класса является количество ненормальных записей, идентифицированных как ненормальные, что является записью во второй строке и втором столбце матрицы неточностей, разделенной на сумму записей во второй строке. Точность для ненормального класса - это доля истинных ненормальных записей в общем количестве, определенном классификатором как ненормальное. Это соответствует записи во второй строке и втором столбце матрицы неточностей, разделенной на сумму записей во втором столбце. Счет F1 является гармоническим средним для двух.

RecallAbnormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(2,:));
PrecisionAbnormal = cmTest.NormalizedValues(2,2)/sum(cmTest.NormalizedValues(:,2));
F1Abnormal = harmmean([RecallAbnormal PrecisionAbnormal]);
fprintf('RecallAbnormal = %2.3f\nPrecisionAbnormal = %2.3f\nF1Abnormal = %2.3f\n',...
    100*RecallAbnormal,100*PrecisionAbnormal,100*F1Abnormal);

RecallAbnormal = 85.904
PrecisionAbnormal = 88.736
F1Abnormal = 87.297

Повторите вышесказанное для нормального класса.

RecallNormal = cmTest.NormalizedValues(1,1)/sum(cmTest.NormalizedValues(1,:));
PrecisionNormal = cmTest.NormalizedValues(1,1)/sum(cmTest.NormalizedValues(:,1));
F1Normal = harmmean([RecallNormal PrecisionNormal]);
fprintf('RecallNormal = %2.3f\nPrecisionNormal = %2.3f\nF1Normal = %2.3f\n',...
    100*RecallNormal,100*PrecisionNormal,100*F1Normal);

RecallNormal = 94.696
PrecisionNormal = 93.248
F1Normal = 93.967

Сводные данные

Этот пример использовал вейвлет рассеяние, чтобы надежно идентифицировать записи фонокардиограммы человека как нормальные или ненормальные в двоичной задаче классификации. Вейвлет рассеяние требовал только спецификации одного параметра, длины инварианта шкалы, в порядок для создания представлений данных PCG с малыми дисперсиями, которые позволили классификатору машины опорных векторов точно смоделировать различие между этими двумя группами. Классификатор машины опорных векторов с вейвлет смог достичь наилучшего решения как в точности, так и в отзыве для обеих групп, несмотря на значительно несбалансированное количество нормальных и ненормальных записей PCG как в обучающем, так и в тестовом наборе.

Ссылки

Вспомогательные функции

partition_heartsounds создает обучающие и тестовые наборы, состоящие из заданных пропорций данных. Функция также сохраняет долю ненормальных и нормальных записей PCG в каждом наборе.

function [trainData, testData, trainLabels, testLabels] = partition_heartsounds(percent_train_split,Data,Labels)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.

% Labels in heart sound data are not sequential.
percent_train_split = percent_train_split/100;
% Each column is an observation
NormalData = Data(:,Labels == 'normal');
AbnormalData = Data(:,Labels == 'abnormal');
LabelsNormal = Labels(Labels == 'normal');
LabelsAbnormal = Labels(Labels == 'abnormal');
Nnormal = size(NormalData,2);
Nabnormal = size(AbnormalData,2);
num_train_normal = round(percent_train_split*Nnormal);
num_train_abnormal = round(percent_train_split*Nabnormal);
rng default;
Pnormal = randperm(Nnormal,num_train_normal);
Pabnormal = randperm(Nabnormal,num_train_abnormal);
notPnormal = setdiff(1:Nnormal,Pnormal);
notPabnormal = setdiff(1:Nabnormal,Pabnormal);
trainNormalData = NormalData(:,Pnormal);
trainNormalLabels = LabelsNormal(Pnormal);
trainAbnormalData = AbnormalData(:,Pabnormal);
trainAbnormalLabels = LabelsAbnormal(Pabnormal);
testNormalData = NormalData(:,notPnormal);
testNormalLabels = LabelsNormal(notPnormal);
testAbnormalData = AbnormalData(:,notPabnormal);
testAbnormalLabels = LabelsAbnormal(notPabnormal);
trainData = [trainNormalData trainAbnormalData];
trainData = (trainData-mean(trainData))./std(trainData,1);
trainLabels = [trainNormalLabels;  trainAbnormalLabels];
testData = [testNormalData  testAbnormalData];
testData = (testData-mean(testData))./std(testData,1);
testLabels = [testNormalLabels; testAbnormalLabels];
   
end

createSequenceLabels_heartsounds создает метки классов для вейвлета временных последовательностей рассеяния.

function [sequence_labels_train,sequence_labels_test] = createSequenceLabels_heartsounds(Nseq,trainLabels,testLabels)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.
Ntrain = numel(trainLabels);
trainLabels = repmat(trainLabels',Nseq,1);
sequence_labels_train = reshape(trainLabels,Nseq*Ntrain,1);
Ntest = numel(testLabels);
testLabels = repmat(testLabels',Nseq,1);
sequence_labels_test = reshape(testLabels,Ntest*Nseq,1);
end

helperMajorityVote реализует большинство голосов для классификации на основе режима. Если уникальный режим отсутствует, голосование NoUniqueMode возвращается для проверки записи ошибки классификации.

function [ClassVotes,ClassCounts] = helperMajorityVote(predLabels,origLabels,classes)
% This function is in support of Wavelet Toolbox examples. It may
% change or be removed in a future release.

% Make categorical arrays if the labels are not already categorical
predLabels = categorical(predLabels);
origLabels = categorical(origLabels);
% Expects both predLabels and origLabels to be categorical vectors
Npred = numel(predLabels);
Norig = numel(origLabels);
Nwin = Npred/Norig;
predLabels = reshape(predLabels,Nwin,Norig);
ClassCounts = countcats(predLabels);
[mxcount,idx] = max(ClassCounts);
ClassVotes = classes(idx);
% Check for any ties in the maximum values and ensure they are marked as
% error if the mode occurs more than once
modecnt = modecount(ClassCounts,mxcount);
ClassVotes(modecnt>1) = categorical({'NoUniqueMode'});
ClassVotes = ClassVotes(:);

%-------------------------------------------------------------------------
function modecnt = modecount(ClassCounts,mxcount)
modecnt = Inf(size(ClassCounts,2),1);
for nc = 1:size(ClassCounts,2)
    modecnt(nc) = histc(ClassCounts(:,nc),mxcount(nc));
end

end

% EOF
end

helperF1heartSounds вычисление точности, отзыва и F1 счетов для результатов классификатора.

function PRTable = helperF1heartSounds(confmat)
% This function is only in support of the Wavelet Time Scattering 
% Classification of Phonocardiogram Data example. It may change or be
% removed in a future release.

precisionAB = confmat(2,2)/sum(confmat(:,2))*100;
precisionNR = confmat(1,1)/sum(confmat(:,1))*100 ;
recallAB = confmat(2,2)/sum(confmat(2,:))*100;
recallNR = confmat(1,1)/sum(confmat(1,:))*100;
F1AB = 2*(precisionAB*recallAB)/(precisionAB+recallAB);
F1NR = 2*(precisionNR*recallNR)/(precisionNR+recallNR);
% Construct a MATLAB Table to display the results.
PRTable = array2table([precisionAB recallAB F1AB;...
    precisionNR recallNR F1NR],...
    'VariableNames',{'Precision','Recall','F1_Score'},'RowNames',...
    {'Abnormal','Normal'});

end