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

В этом примере используются данные фонокардиограммы (ПХГ), полученные от людей с нормальной и аномальной сердечной функцией. Набор данных состоит из 3829 записей, 2575 записей от людей с нормальной сердечной функцией и 1254 записей от людей с аномальной сердечной функцией. Каждая запись имеет длину 10000 выборок и дискретизируется при частоте 2 кГц. Это представляет собой пять секунд данных фонокардиограммы. Набор данных составлен на основе данных обучения и валидации, используемых в PhysioNet Computing 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. Имеется 279 путей рассеяния и пять окон рассеяния для каждого из 2680 сигналов. Чтобы передать это в 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 к учебным данным мы выполняем пятикратную перекрестную проверку для оценки ошибки обобщения учебных данных. Здесь каждое окно рассеяния классифицируется отдельно.

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

В задаче классификации точность для класса - это количество правильных положительных результатов, деленное на число положительных результатов. Другими словами, из всех записей, которым классификатор присваивает данную метку, какая пропорция на самом деле принадлежит классу. Recall определяется как количество правильных меток, деленное на количество меток для данного класса. В частности, из всех записей, принадлежащих классу, какая доля сделала наш классификатор меткой этого класса. Оценивая точность вашего классификатора, вы в идеале хотите добиться успеха как на точности, так и на отзыве. Например, предположим, что у нас был классификатор, который помечал каждую запись 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