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

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

Загрузите данные

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

Файл physionet_phonocardiogram-master.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-master.zip'),tempdir)
unzip(fullfile(tempdir,'physionet_phonocardiogram-master','PCG_Data.zip'),...
    fullfile(tempdir,'PCG_Data'))

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

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

heartSoundData массив структур с двумя полями: Data и Classesданные 10000 3829 матрица, где каждый столбец является записью PCG. Classes 3829 1 категориальный массив диагностических меток, один для каждого столбца Data. Поскольку это - бинарная проблема классификации, классы являются "нормальными" и "аварийными". Как ранее утверждено, существует 2 575 нормальных записей и 1 254 аварийных записи. Эквивалентно, 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, делит 3 829 наблюдений так, чтобы 70% (2680) были в наборе обучающих данных с 1 802 нормальными и 878 аварийными. Остающиеся 1 149 записей (773 нормальных и 376 аварийных) протянуты в наборе тестов для предсказания. Генератор случайных чисел отобран в функции помощника, таким образом, результаты повторяемы. Код для partition_heartsounds и все другие функции помощника, используемые в этом примере, даны в разделе Supporting Functions в конце примера.

[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

Рассеивание функций

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

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

Для данных рассеивающихся параметров, scat_features_train 274 5 2 680 матрицами. Существует 274 рассеивающихся пути и пять рассеивающихся окон для каждого из 2 680 сигналов. Для пересылки/передачи это классификатору SVM, измените тензор в 13400 274 матрица, где каждая строка представляет одно окно рассеивания через 274 рассеивающихся пути. Общее количество строк равно продукту 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.98 percent

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

Accuracy is 99.02 percent

confmatCV = confusionchart(sequence_labels_train,predLabels);

Обратите внимание на то, что рассеивающаяся сеть приводит приблизительно к 99-процентной точности, когда каждый раз окно классифицируется отдельно. Однако эффективность на самом деле лучше, чем это значение, потому что у нас есть пять рассеивающихся окон на запись, и точность на 99,02 процентов основана на классификации всех окон отдельно. В этом случае используйте решение большинством голосов, чтобы получить присвоение единого класса на запись. Голосование класса соответствует режиму голосов за эти пять окон. Если никакой уникальный режим не найден, 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(ClassVotes,trainLabels);

Точность перекрестной проверки на обучающих данных - на самом деле 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 92.25 percent.

cmTest = confusionchart(ClassVotes,testLabels);

Из 1 149 тестовых записей 92,25% правильно классифицируются как "Нормальные" или "Аварийные". Из 773 нормальных записей PCG в наборе тестов, 732 правильно классифицируются. Из 376 аварийных записей в наборе тестов, 328 правильно классифицируются.

Точность, вспомните, и счет F1

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

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

                Precision    Recall    F1_Score
                _________    ______    ________

    Abnormal     87.234      88.889     88.054 
    Normal       94.696      93.846     94.269 

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

cmTest = confusionchart(ClassVotes,testLabels);

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

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

RecallAbnormal = 88.889
PrecisionAbnormal = 87.234
F1Abnormal = 88.054

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

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

RecallNormal = 93.846
PrecisionNormal = 94.696
F1Normal = 94.269

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

Этот пример использовал время вейвлета, рассеиваясь, чтобы надежно идентифицировать человеческие записи фонокардиограммы как нормальные или аварийные в бинарной проблеме классификации. Вейвлет, рассеивающийся, потребовал только спецификации одного параметра, длины инвариантного к масштабу, для того, чтобы произвести представления низкого отклонения данных 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(1,1)/sum(confmat(:,1))*100;
precisionNR = confmat(2,2)/sum(confmat(:,2))*100 ;
recallAB = confmat(1,1)/sum(confmat(1,:))*100;
recallNR = confmat(2,2)/sum(confmat(2,:))*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