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

В этом примере используются данные ЭКГ, полученные из трех групп или классов людей: людей с сердечной аритмией, людей с застойной сердечной недостаточностью и людей с нормальным синусовым ритмом. В примере используются 162 записи ЭКГ из трех баз данных PhysioNet: база данных аритмии MIT-BIH [3] [7], база данных нормального синусового ритма MIT-BIH [3] и база данных застойной сердечной недостаточности BIDMC [1] [3]. Всего насчитывается 96 записей от лиц с аритмией, 30 записей от лиц с застойной сердечной недостаточностью и 36 записей от лиц с нормальными синусовыми ритмами. Цель состоит в том, чтобы обучить классификатор различать аритмию (ARR), застойную сердечную недостаточность (CHF) и нормальный синусовый ритм (NSR).

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

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

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

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

ECGData.mat содержит данные, используемые в этом примере. Файл .txt, Modified_physionet_data.txt, требуется политикой копирования PhysioNet и предоставляет исходные атрибуты для данных, а также описание этапов предварительной обработки, применяемых к каждой записи ЭКГ.

Загрузить файлы

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

unzip(fullfile(tempdir,'physionet_ECG_data-main.zip'),tempdir)
unzip(fullfile(tempdir,'physionet_ECG_data-main','ECGData.zip'),...
    fullfile(tempdir,'ECGData'))

После распаковки файла ECGData.zip загрузите данные в MATLAB.

load(fullfile(tempdir,'ECGData','ECGData.mat'))

ECGData - структурный массив с двумя полями: Data и Labels. Data является матрицей 162 на 65536, где каждая строка представляет собой запись ЭКГ, дискретизированную при 128 герц. Labels представляет собой массив диагностических меток типа 162 на 1, по одной для каждой строки Data. Три диагностические категории: «ARR» (аритмия), «CHF» (застойная сердечная недостаточность) и «NSR» (нормальный синусовый ритм).

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

Случайное разделение данных на два набора - обучающие и тестовые наборы данных. Вспомогательная функция helperRandomSplit выполняет случайное разделение. helperRandomSplit принимает требуемый процент разделения для данных обучения и ECGData. helperRandomSplit функция выводит два набора данных вместе с набором меток для каждого. Каждая строка trainData и testData - сигнал ЭКГ. Каждый элемент trainLabels и testLabels содержит метку класса для соответствующей строки матриц данных. В этом примере мы случайным образом присваиваем обучающему набору 70% данных в каждом классе. Оставшиеся 30% удерживаются для тестирования (прогнозирования) и назначаются тестовому набору.

percent_train = 70;
[trainData,testData,trainLabels,testLabels] = ...
    helperRandomSplit(percent_train,ECGData);

Есть 113 записей в trainData набор и 49 рекордов в testData. По проекту данные обучения содержат 69,75% (113/162) данных. Напомним, что класс ARR представляет 59,26% данных (96/162), класс CHF - 18,52% (30/162), а класс NSR - 22,22% (36/162). Проверьте процент каждого класса в обучающих и тестовых наборах. Проценты в каждом из них согласуются с общими процентами классов в наборе данных.

Ctrain = countcats(categorical(trainLabels))./numel(trainLabels).*100

Ctrain = 3×1

   59.2920
   18.5841
   22.1239

Ctest = countcats(categorical(testLabels))./numel(testLabels).*100

Ctest = 3×1

   59.1837
   18.3673
   22.4490

Образцы графика

Постройте первые несколько тысяч образцов четырех случайно выбранных записей из ECGData. Вспомогательная функция helperPlotRandomRecords делает это. helperPlotRandomRecords принимает ECGData и случайное начальное число в качестве входных данных. Начальное начальное число устанавливается равным 14, так что по меньшей мере одна запись из каждого класса наносится на график. Можно выполнить helperPlotRandomRecords с ECGData в качестве единственного входного аргумента столько раз, сколько вы хотите получить представление о многообразии форм сигналов ЭКГ, связанных с каждым классом. Исходный код для этой вспомогательной функции можно найти в разделе «Вспомогательные функции» в конце этого примера.

helperPlotRandomRecords(ECGData,14)

Извлечение элементов

Извлеките характеристики, используемые в классификации сигналов для каждого сигнала. В этом примере используются следующие признаки, извлеченные из 8 блоков каждого сигнала длительностью приблизительно одну минуту (8192 выборки):

Дополнительно для каждого сигнала по всей длине данных [6] извлекают многомасштабные оценки вейвлет-дисперсии. Используется несмещенная оценка вейвлет-дисперсии. Это требует, чтобы в оценках дисперсии использовались только уровни с, по меньшей мере, одним вейвлет-коэффициентом, не затронутым граничными условиями. Для длины сигнала 2 ^ 16 (65 536) и вейвлета «db2» это приводит к 14 уровням.

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

Коэффициенты AR для каждого окна оцениваются с использованием метода Бурга, arburg. В [8] авторы использовали методы выбора порядка моделей, чтобы определить, что модель AR (4) обеспечивает наилучшее соответствие для форм сигналов ЭКГ в аналогичной проблеме классификации. В [5] вычислена информационная теоретическая мера, энтропия Шеннона, на терминальных узлах вейвлет-дерева пакетов и использована со случайным классификатором леса. Здесь мы используем недекимированное вейвлет-пакетное преобразование ,modwpt, до уровня 4.

Определение энтропии Шеннона для недекимированного вейвлет-пакетного преобразования после [5] дано следующим образом: $${}_{}\underset{\mathrm{}}{\overset{}{}}{}_{}\mathrm{}{}_{\mathrm{SEj=-\sum k=1Npj,k*logpj,k}}$$, где $$N$$ - число соответствующих коэффициентов в j-м узле и $${}_{\mathrm{pj},}$$k - нормализованные квадраты коэффициентов вейвлет-пакета в j-м терминальном узле.

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

Дисперсию импульса для всего сигнала получают, используя modwtvar. Вейвлет-дисперсия измеряет изменчивость сигнала по масштабу или эквивалентную изменчивость сигнала по частотным интервалам октавной полосы.

Всего имеется 190 признаков: 32 признака AR (4 коэффициента на блок), 128 значений энтропии Шеннона (16 значений на блок), 16 фрактальных оценок (2 на блок) и 14 оценок вейвлет-дисперсии.

helperExtractFeatures функция вычисляет эти признаки и объединяет их в вектор признаков для каждого сигнала. Исходный код для этой вспомогательной функции можно найти в разделе «Вспомогательные функции» в конце этого примера.

timeWindow = 8192;
ARorder = 4;
MODWPTlevel = 4;
[trainFeatures,testFeatures,featureindices] = ...
    helperExtractFeatures(trainData,testData,timeWindow,ARorder,MODWPTlevel);

trainFeatures и testFeatures представляют собой матрицы 113 на 190 и 49 на 190 соответственно. Каждая строка этих матриц является характерным вектором для соответствующих данных ЭКГ в trainData и testDataсоответственно. При создании векторов признаков данные уменьшаются с 65536 выборок до 190 векторов элементов. Это значительное сокращение данных, но цель не просто сокращение данных. Цель состоит в том, чтобы свести данные к гораздо меньшему набору признаков, который фиксирует разницу между классами, чтобы классификатор мог точно разделять сигналы. Индексы для элементов, которые составляют оба trainFeatures и testFeatures содержатся в массиве структуры, featureindices. Эти индексы можно использовать для просмотра элементов по группам. В качестве примера рассмотрим диапазон степеней Холдера в спектрах сингулярности для первого временного окна. Постройте график данных для всего набора данных.

allFeatures = [trainFeatures;testFeatures];
allLabels = [trainLabels;testLabels];
figure
boxplot(allFeatures(:,featureindices.HRfeatures(1)),allLabels,'notch','on')
ylabel('Holder Exponent Range')
title('Range of Singularity Spectrum by Group (First Time Window)')
grid on

Для этой функции можно выполнить односторонний анализ дисперсии и подтвердить то, что появляется в boxplot, а именно, что группы ARR и NSR имеют значительно больший диапазон, чем группа CHF.

[p,anovatab,st] = anova1(allFeatures(:,featureindices.HRfeatures(1)),...
    allLabels);

c = multcompare(st,'display','off')

c = 3×6

    1.0000    2.0000    0.0176    0.1144    0.2112    0.0155
    1.0000    3.0000   -0.1591   -0.0687    0.0218    0.1764
    2.0000    3.0000   -0.2975   -0.1831   -0.0687    0.0005

В качестве дополнительного примера рассмотрим дисперсию в поддиапазоне вейвлетов второй и самой низкой частот (второй и самой большой шкалы) для трех групп.

boxplot(allFeatures(:,featureindices.WVARfeatures(end-1)),allLabels,'notch','on')
ylabel('Wavelet Variance')
title('Wavelet Variance by Group')
grid on

При выполнении анализа дисперсии для этой функции обнаруживается, что группа NSR имеет значительно меньшую дисперсию в этом вейвлет-поддиапазоне, чем группы ARR и CHF. Эти примеры предназначены только для иллюстрации того, как отдельные особенности служат для разделения классов. Хотя одного элемента недостаточно, цель состоит в том, чтобы получить достаточно богатый набор элементов, позволяющий классификатору разделить все три класса.

Классификация сигналов

Теперь, когда данные сведены к вектору признаков для каждого сигнала, следующим шагом является использование этих векторов признаков для классификации сигналов ЭКГ. Приложение Classification Learner позволяет быстро оценить большое количество классификаторов. В этом примере используется многоклассный SVM с квадратичным ядром. Проводят два анализа. Сначала мы используем весь набор данных (наборы обучения и тестирования) и оцениваем частоту неправильной классификации и матрицу путаницы, используя пятикратную перекрестную проверку.

features = [trainFeatures; testFeatures];
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial',...
    'PolynomialOrder',2,...
    'KernelScale','auto',...
    'BoxConstraint',1,...
    'Standardize',true);
model = fitcecoc(...
    features,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100

loss = 8.0247

[confmatCV,grouporder] = confusionmat([trainLabels;testLabels],classLabels);

Пятикратная ошибка классификации является 8,02% (91,98% верной). Матрица путаницы, confmatCVпоказывает, какие записи были неправильно классифицированы. grouporder дает порядок групп. Две из группы ARR были неправильно классифицированы как CHF, восемь из группы CHF были неправильно классифицированы как ARR и одна как NSR, а две из группы NSR были неправильно классифицированы как ARR.

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

В задаче классификации точность для класса - это количество правильных положительных результатов, деленное на число положительных результатов. Другими словами, из всех записей, которым классификатор присваивает данную метку, какая пропорция на самом деле принадлежит классу. Recall определяется как количество правильных меток, деленное на количество меток для данного класса. В частности, из всех записей, принадлежащих классу, какая доля сделала наш классификатор меткой этого класса. Оценивая точность вашей системы машинного обучения, вы в идеале хотите добиться успеха как на точности, так и на отзыве. Например, предположим, что у нас был классификатор, который помечал каждую запись как ARR. Тогда наш отзыв для класса ARR был бы 1 (100%). Все записи, принадлежащие классу ARR, будут помечены как ARR. Однако точность будет низкой. Поскольку наш классификатор помечал все записи как ARR, в этом случае было бы 66 ложных срабатываний для точности 96/162 или 0,5926. Оценка F1 представляет собой среднее гармоническое точности и отзыва и поэтому обеспечивает единую метрику, которая суммирует производительность классификатора как с точки зрения отзыва, так и с точки зрения точности. Следующая вспомогательная функция вычисляет показатели точности, отзыва и F1 для трех классов. Вы можете видеть, как helperPrecisionRecall вычисляет точность, отзыв и оценку F1 на основе матрицы путаницы путем анализа кода в разделе «Вспомогательные функции».

CVTable = helperPrecisionRecall(confmatCV);

Можно просмотреть таблицу, возвращенную helperPrecisionRecall с помощью следующей команды.

disp(CVTable)

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR     90.385      97.917         94 
    CHF     91.304          70     79.245 
    NSR     97.143      94.444     95.775 

Точность и повторный вызов хороши для классов ARR и NSR, в то время как повторный вызов значительно ниже для класса CHF.

Для следующего анализа мы подгоняем многоклассный квадратичный SVM только к учебным данным (70%), а затем используем эту модель для прогнозирования 30% данных, задержанных для тестирования. В наборе тестов имеется 49 записей данных.

model = fitcecoc(...
     trainFeatures,...
     trainLabels,...
     'Learners',template,...
     'Coding','onevsone',...
     'ClassNames',{'ARR','CHF','NSR'});
predLabels = predict(model,testFeatures);

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

correctPredictions = strcmp(predLabels,testLabels);
testAccuracy = sum(correctPredictions)/length(testLabels)*100

testAccuracy = 97.9592

[confmatTest,grouporder] = confusionmat(testLabels,predLabels);

Точность классификации в наборе тестовых данных составляет приблизительно 98%, и матрица путаницы показывает, что одна запись CHF была неправильно классифицирована как NSR.

Аналогично тому, что было сделано в анализе перекрестной проверки, получите точность, отзыв и оценки F1 для тестового набора.

testTable = helperPrecisionRecall(confmatTest);
disp(testTable)

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR        100         100        100 
    CHF        100      88.889     94.118 
    NSR     91.667         100     95.652 

Классификация по необработанным данным и кластеризация

Из предыдущего анализа возникают два естественных вопроса. Необходима ли экстракция признаков для достижения хороших результатов классификации? Необходим ли классификатор или эти функции могут разделить группы без классификатора? Для решения первого вопроса повторите результаты перекрестной проверки для необработанных данных временных рядов. Следует отметить, что следующий этап является дорогостоящим в вычислительном отношении, поскольку мы применяем SVM к матрице 162 на 65536. Если вы не хотите выполнять этот шаг самостоятельно, результаты будут описаны в следующем абзаце.

rawData = [trainData;testData];
Labels = [trainLabels;testLabels];
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial', ...
    'PolynomialOrder',2, ...
    'KernelScale','auto', ...
    'BoxConstraint',1, ...
    'Standardize',true);
model = fitcecoc(...
    rawData,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100

loss = 33.3333

[confmatCVraw,grouporder] = confusionmat([trainLabels;testLabels],classLabels);
rawTable = helperPrecisionRecall(confmatCVraw);
disp(rawTable)

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR        64          100     78.049 
    CHF       100       13.333     23.529 
    NSR       100       22.222     36.364 

Коэффициент неправильной классификации для необработанных данных временных рядов составляет 33,3%. Повторение точности, напомним, и F1 анализа баллов показывает очень плохие оценки F1 как для групп CHF (23,52), так и для групп NSR (36,36). Получить коэффициенты дискретного преобразования Фурье (DFT) амплитуды для каждого сигнала, чтобы выполнить анализ в частотной области. Поскольку данные являются действительными, мы можем добиться некоторого уменьшения данных, используя DFT, используя тот факт, что величины Фурье являются равномерной функцией.

rawDataDFT = abs(fft(rawData,[],2));
rawDataDFT = rawDataDFT(:,1:2^16/2+1);
rng(1)
template = templateSVM(...
    'KernelFunction','polynomial',...
    'PolynomialOrder',2,...
    'KernelScale','auto',...
    'BoxConstraint',1,...
    'Standardize',true);
model = fitcecoc(...
    rawDataDFT,...
    [trainLabels;testLabels],...
    'Learners',template,...
    'Coding','onevsone',...
    'ClassNames',{'ARR','CHF','NSR'});
kfoldmodel = crossval(model,'KFold',5);
classLabels = kfoldPredict(kfoldmodel);
loss = kfoldLoss(kfoldmodel)*100

loss = 19.1358

[confmatCVDFT,grouporder] = confusionmat([trainLabels;testLabels],classLabels);
dftTable = helperPrecisionRecall(confmatCVDFT);
disp(dftTable)

           Precision    Recall    F1_Score
           _________    ______    ________

    ARR     76.423      97.917     85.845 
    CHF        100      26.667     42.105 
    NSR     93.548      80.556     86.567 

Использование значений DFT снижает скорость неправильной классификации до 19,13%, но это все еще более чем в два раза превышает частоту ошибок, полученную с нашими 190 характеристиками. Эти анализы показывают, что классификатор получил пользу от тщательного выбора признаков.

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

rng default
eva = evalclusters(features,'kmeans','gap','KList',[1:6]);
eva

eva = 
  GapEvaluation with properties:

    NumObservations: 162
         InspectedK: [1 2 3 4 5 6]
    CriterionValues: [1.2777 1.3539 1.3644 1.3570 1.3591 1.3752]
           OptimalK: 3

Статистика разрыва показывает, что оптимальное количество кластеров - три. Тем не менее, если вы посмотрите на количество записей в каждом из трех кластеров, вы увидите, что кластеризация k-означает, основанная на векторах признаков, сделала плохую работу по разделению трех диагностических категорий.

countcats(categorical(eva.OptimalY))

ans = 3×1

    61
    74
    27

Напомним, что в классе ARR 96 человек, в классе CHF - 30, в классе СМП - 36 человек.

Резюме

В этом примере использовалась обработка сигналов для извлечения вейвлет-признаков из ЭКГ-сигналов и эти признаки для классификации ЭКГ-сигналов по трем классам. Извлечение признаков не только привело к значительному сокращению данных, но и захватило различия между классами ARR, CHF и NSR, что продемонстрировано результатами перекрестной проверки и производительностью классификатора SVM в тестовом наборе. Пример далее продемонстрировал, что применение классификатора SVM к необработанным данным приводило к низкой производительности, как и кластеризация векторов признаков без использования классификатора. Ни классификатор, ни только признаки не были достаточными для разделения классов. Однако, когда извлечение признаков использовалось в качестве этапа сокращения данных до использования классификатора, три класса были хорошо разделены.

Ссылки

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

helperPlotRandomRecords строит графики четырех сигналов ЭКГ, случайно выбранных из ECGData.

function helperPlotRandomRecords(ECGData,randomSeed)
% This function is only intended to support the XpwWaveletMLExample. It may
% change or be removed in a future release.

if nargin==2
    rng(randomSeed)
end

M = size(ECGData.Data,1);
idxsel = randperm(M,4);
for numplot = 1:4
    subplot(2,2,numplot)
    plot(ECGData.Data(idxsel(numplot),1:3000))
    ylabel('Volts')
    if numplot > 2
        xlabel('Samples')
    end
    title(ECGData.Labels{idxsel(numplot)})
end

end

helperExtractFeatures Извлекает вейвлет-признаки и коэффициенты AR для блоков данных заданного размера. Элементы объединяются в векторы элементов.

function [trainFeatures, testFeatures,featureindices] = helperExtractFeatures(trainData,testData,T,AR_order,level)
% This function is only in support of XpwWaveletMLExample. It may change or
% be removed in a future release.
trainFeatures = [];
testFeatures = [];

for idx =1:size(trainData,1)
    x = trainData(idx,:);
    x = detrend(x,0);
    arcoefs = blockAR(x,AR_order,T);
    se = shannonEntropy(x,T,level);
    [cp,rh] = leaders(x,T);
    wvar = modwtvar(modwt(x,'db2'),'db2');
    trainFeatures = [trainFeatures; arcoefs se cp rh wvar']; %#ok<AGROW>
    
end

for idx =1:size(testData,1)
    x1 = testData(idx,:);
    x1 = detrend(x1,0);
    arcoefs = blockAR(x1,AR_order,T);
    se = shannonEntropy(x1,T,level);
    [cp,rh] = leaders(x1,T);
    wvar = modwtvar(modwt(x1,'db2'),'db2');
    testFeatures = [testFeatures;arcoefs se cp rh wvar']; %#ok<AGROW>
    
end

featureindices = struct();
% 4*8
featureindices.ARfeatures = 1:32;
startidx = 33;
endidx = 33+(16*8)-1;
featureindices.SEfeatures = startidx:endidx;
startidx = endidx+1;
endidx = startidx+7;
featureindices.CP2features = startidx:endidx;
startidx = endidx+1;
endidx = startidx+7;
featureindices.HRfeatures = startidx:endidx;
startidx = endidx+1;
endidx = startidx+13;
featureindices.WVARfeatures = startidx:endidx;
end


function se = shannonEntropy(x,numbuffer,level)
numwindows = numel(x)/numbuffer;
y = buffer(x,numbuffer);
se = zeros(2^level,size(y,2));
for kk = 1:size(y,2)
    wpt = modwpt(y(:,kk),level);
    % Sum across time
    E = sum(wpt.^2,2);
    Pij = wpt.^2./E;
    % The following is eps(1)
    se(:,kk) = -sum(Pij.*log(Pij+eps),2);
end
se = reshape(se,2^level*numwindows,1);
se = se';
end


function arcfs = blockAR(x,order,numbuffer)
numwindows = numel(x)/numbuffer;
y = buffer(x,numbuffer);
arcfs = zeros(order,size(y,2));
for kk = 1:size(y,2)
    artmp =  arburg(y(:,kk),order);
    arcfs(:,kk) = artmp(2:end);
end
arcfs = reshape(arcfs,order*numwindows,1);
arcfs = arcfs';
end


function [cp,rh] = leaders(x,numbuffer)
y = buffer(x,numbuffer);
cp = zeros(1,size(y,2));
rh = zeros(1,size(y,2));
for kk = 1:size(y,2)
    [~,h,cptmp] = dwtleader(y(:,kk));
    cp(kk) = cptmp(2);
    rh(kk) = range(h);
end
end

helperPrecisionRecall возвращает оценки точности, отзыва и F1 на основе матрицы путаницы. Вывод результатов в виде таблицы MATLAB.

function PRTable = helperPrecisionRecall(confmat)
% This function is only in support of XpwWaveletMLExample. It may change or
% be removed in a future release.
precisionARR = confmat(1,1)/sum(confmat(:,1))*100;
precisionCHF = confmat(2,2)/sum(confmat(:,2))*100 ;
precisionNSR = confmat(3,3)/sum(confmat(:,3))*100 ;
recallARR = confmat(1,1)/sum(confmat(1,:))*100;
recallCHF = confmat(2,2)/sum(confmat(2,:))*100;
recallNSR = confmat(3,3)/sum(confmat(3,:))*100;
F1ARR = 2*precisionARR*recallARR/(precisionARR+recallARR);
F1CHF = 2*precisionCHF*recallCHF/(precisionCHF+recallCHF);
F1NSR = 2*precisionNSR*recallNSR/(precisionNSR+recallNSR);
% Construct a MATLAB Table to display the results.
PRTable = array2table([precisionARR recallARR F1ARR;...
    precisionCHF recallCHF F1CHF; precisionNSR recallNSR...
    F1NSR],'VariableNames',{'Precision','Recall','F1_Score'},'RowNames',...
    {'ARR','CHF','NSR'});

end