睡眠阶段分期——SVM和ELM分别与粒子群算法结合（function）
生活随笔
收集整理的這篇文章主要介紹了
睡眠阶段分期——SVM和ELM分别与粒子群算法结合（function）
小編覺得挺不錯的,現(xiàn)在分享給大家,幫大家做個參考.
目錄
importslpdb函數(shù)
loadmatobject函數(shù)
extractfeatures函數(shù)
getindexrange函數(shù)
PSOforELM函數(shù)
?PSOforSVM函數(shù)
?extractresults函數(shù)
importslpdb函數(shù)
——將頭文件、RR文件、注釋文件構(gòu)成結(jié)構(gòu)體：OutputData = struct('filename', fileName, 'time', anTimeGeneratedCell, ...'rr', rrCollection, 'annotation', anClassGeneratedCell, 'age', age, ...'gender', gender, 'weight', weight);
function OutputData = importslpdb(fileName) %Import and synchronize a slpdb recording % Syntax: % OutputData = importslpdb(fileName) % % Input: % *) fileName - slpdb file name to be imported. Example: 'slp01a'. % file must be located in 'slpdb' folder, % three file formats needed: .hea, .rr, and .an % % Output: % *) OutputData - struct contains synchronized RR and annotationfileName = strcat('slpdb/', cell2mat(fileName)); SEC_PER_EPOCH = 30; % amount of seconds in one epoch (value for slpdb is 30) OutputData = []; fprintf('\n%s DATA IMPORT...\n', fileName); %% IMPORT HEADER DATA 導(dǎo)入頭文件數(shù)據(jù)fprintf('Importing header file...\n'); fid = fopen(strcat(fileName, '.hea'), 'r'); headerFile = textscan(fid, '%s%s%s%s%s%s%s%*[^\n]'); fclose(fid); heaSamplingFreq = strsplit(char(headerFile{3}(1)), '/'); heaSamplingFreq = str2double(cell2mat(heaSamplingFreq(1))); heaTotalSamples = str2double(cell2mat(headerFile{4}(1))); heaRecLengthInSec = ceil(heaTotalSamples/heaSamplingFreq); heaTotalEpoch = ceil(heaRecLengthInSec/SEC_PER_EPOCH); heaIdx = size(headerFile{1}, 1); % get last line index of fileif cell2mat(headerFile{1}(end-1)) == '#'% decease index by 1 for 'slp37.hea',% because the last line is not age, gender, and weight information heaIdx = heaIdx - 1; endage = headerFile{2}(heaIdx); gender = headerFile{3}(heaIdx); weight = headerFile{4}(heaIdx);% output of "IMPORT HEADER DATA" section: % *) heaTotalEpoch - total epoch according to header data % *) age - age of the subject % *) gender - gender of the subject % *) weight - weight of the subject % END OF IMPORT HEADER DATA %% IMPORT ANNOTATION DATA 導(dǎo)入注釋數(shù)據(jù)fprintf('Importing annotation file...\n'); fid = fopen(strcat(fileName, '.an'), 'r'); anFile = textscan(fid, '%s%s%s%s%s%s%s%*[^\n]'); fclose(fid);% remove header of annotation data (first line) % remove this string: 'Elapsed time Sample # Type Sub Chan Num Aux'for i=1:size(anFile, 2)anFile{i}(1) = []; end% change first epoch's 'start time' into 0:00.000 anTemp = cell2mat(anFile{1}(1)); anTemp(end-2:end) = 48; anFile{1}(1) = cellstr(anTemp); anTime = anFile{1}; anClass = anFile{7};% output of "IMPORT ANNOTATION DATA" section: % *) anTime - time from annotation file (cell array) % *) anClass - annotation (cell array) % END OF IMPORT ANNOTATION DATA %% IMPORT RR DATA 導(dǎo)入RR數(shù)據(jù)fprintf('Importing RR file...\n'); fid = fopen(strcat(fileName, '.rr'), 'r'); rrFile = textscan(fid, '%s%s%s%s%s%*[^\n]'); fclose(fid);rrConvertedTime = rrFile{1};for i=1:size(rrConvertedTime, 1);rrStartTimeChar = cell2mat(rrConvertedTime(i)); % convert cell into charrrStartTimeChar(end-2:end) = 48; % xx:xx:xx.aaa -> change 'aaa' part to '000'% split start time by ":" into matrix 求時長rrStartTimeMat = strsplit(char(rrConvertedTime(i)), ':')'; % get seconds from the last element rrSecond = str2double(cell2mat(rrStartTimeMat(end)));% epoch groupingrrWhichGroup = floor(rrSecond/SEC_PER_EPOCH)*SEC_PER_EPOCH;% set epoch groupingif rrWhichGroup == 0rrStartTimeChar(end-5) = 48;elseif rrWhichGroup == 30rrStartTimeChar(end-5) = 51;endrrStartTimeChar(end-4) = 48;rrConvertedTime(i) = mat2cell(rrStartTimeChar, 1);end% change RR value from 'array of cell' into 'array of double' 改變RR值屬性rrNum = zeros(size(rrFile{3}, 1), 1);for i=1:size(rrFile{3}, 1)rrNum(i) = str2double(cell2mat(rrFile{3}(i))); endrrTime = rrFile{1};% output of "IMPORT RR DATA" section: % *) rrConvertedTime - rounded RR start time according to the epoch % example: 1:34:31.328 -> 1:34:30.000 % 1:50:13.616 -> 1:50:00.000 % *) rrNum - RR value (array of double) % *) rrTime - start time of un-rounded RR % END OF IMPORT RR DATA %% VALIDITY CHECK 數(shù)據(jù)有效性檢驗（A.注釋、B.RR數(shù)據(jù)）fprintf('Data Validity Check:\n');% A. Annotation File Check % *) generate annotation time acocrding to total epoch from header file % *) result: anTimeGeneratedMat -> Matrix size: number of epoch X 3 (h,m,s)anTimeGeneratedMat = zeros(heaTotalEpoch, 3);for i=2:heaTotalEpochanTimeGeneratedMat(i, 3) = anTimeGeneratedMat(i-1, 3) + SEC_PER_EPOCH;anTimeGeneratedMat(i, 2) = anTimeGeneratedMat(i-1, 2);anTimeGeneratedMat(i, 1) = anTimeGeneratedMat(i-1, 1);if anTimeGeneratedMat(i, 3) >= 60anTimeGeneratedMat(i, 3) = 0;anTimeGeneratedMat(i, 2) = anTimeGeneratedMat(i-1, 2) + 1;if anTimeGeneratedMat(i, 2) >= 60anTimeGeneratedMat(i, 2) = 0;anTimeGeneratedMat(i, 1) = anTimeGeneratedMat(i-1, 1) + 1;endendend% convert anTimeGeneratedMat into anTimeGeneratedCell for easier comparisonanTimeGeneratedCell = cell(size(anTimeGeneratedMat, 1), 1);for i = 1:size(anTimeGeneratedMat, 1)if anTimeGeneratedMat(i, 1) == 0 % when the 'hour' is 0anTimeGeneratedCell(i) = ...cellstr(strcat(sprintf('%d',anTimeGeneratedMat(i, 2)), ...sprintf(':%02d.000',anTimeGeneratedMat(i, 3))));elsetemp = strcat(sprintf('%d', anTimeGeneratedMat(i, 1)), ...sprintf(':%02d',anTimeGeneratedMat(i, 2)));anTimeGeneratedCell(i) = ...cellstr(strcat(temp, sprintf(':%02d.000',anTimeGeneratedMat(i, 3))));end end% *) ANNOTATION FILE CHECK 1 (Total epoch of each data): fprintf(' CHECK 1: ');if heaTotalEpoch == size(anTime, 1) && ...size(unique(rrConvertedTime), 1) == size(anTimeGeneratedMat, 1) && ...heaTotalEpoch == size(unique(rrConvertedTime), 1)fprintf('[SUCCESS] heaTotalEpoch (%d) == size(anTime, 1) (%d) == ', ...'size(unique(rrConvertedTime), 1) (%d) == ', ...'size(anTimeGeneratedMat, 1) (%d)\n', heaTotalEpoch, size(anTime, 1), ...size(unique(rrConvertedTime), 1), size(anTimeGeneratedMat, 1)); elsefprintf('[WARNING] heaTotalEpoch (%d) != size(anTime, 1) (%d) != ', ...'size(unique(rrConvertedTime), 1) (%d) != ', ...'size(anTimeGeneratedMat, 1) (%d)\n', heaTotalEpoch, size(anTime, 1), ...size(unique(rrConvertedTime), 1), size(anTimeGeneratedMat, 1)); end% *) ANNOTATION FILE CHECK 2 (Check equality of anTimeGeneratedCell and anTime): fprintf(' CHECK 2: '); if size(anTime, 1) == heaTotalEpochfor i=1:heaTotalEpochif ~strcmp(anTimeGeneratedCell{i}, anTime{i})fprintf('[FAILED ] anTimeGeneratedCell is NOT EQUAL to anTime\n');returnendendfprintf('[SUCCESS] anTimeGeneratedCell is EQUAL to anTime\n'); elsefprintf('[WARNING] size(anTime, 1) (%d) != heaTotalEpoch (%d), ', ...'anTimeGeneratedCell will be used\n', size(anTime, 1), heaTotalEpoch);end% *) ANNOTATION FILE CHECK 3 (Check annotation value must be '1', '2', '3','4', 'W', 'R', or {'MT', 'M' -> these two will be removed later}): fprintf(' CHECK 3: '); distinctClass = char(unique(anClass));for i=1:size(distinctClass, 1)if distinctClass(i) ~= '1' && distinctClass(i) ~= '2' ...&& distinctClass(i) ~= '3' && distinctClass(i) ~= '4' ...&& distinctClass(i) ~= 'W' && distinctClass(i) ~= 'R' ...&& distinctClass(i) ~= 'M'fprintf('[WARNING ] Annotation values is NOT OK\n');returnendendfprintf('[SUCCESS] Annotation values is OK\n');% B. RR File Check % *) RR FILE CHECK 1 (Check equality of size(unique(rrConvertedTime), 1) % and heaTotalEpoch):fprintf(' CHECK 4: ');if size(unique(rrConvertedTime), 1) ~= heaTotalEpochfprintf('[WARNING] size(unique(rrConvertedTime), 1) (%d) != ', ...'heaTotalEpoch (%d)\n', size(unique(rrConvertedTime), 1), heaTotalEpoch); elsefprintf('[SUCCESS] size(unique(rrConvertedTime), 1) (%d) == ', ...'heaTotalEpoch (%d)\n', size(unique(rrConvertedTime), 1), heaTotalEpoch); end% END OF VALIDITY CHECK %% SYNCHRONIZE RR AND ANNOTATION DATA 同步RR和注釋數(shù)據(jù)epochCounter = 1; rrCounter = 1;% rrCollection = each row contains RRs of associated epoch rrCollection = cell(heaTotalEpoch, 1);% rrTimeCollection = each row contains RR time of associated epoch rrTimeCollection = cell(heaTotalEpoch, 1);for i=1:size(rrConvertedTime, 1) % looping for each rrConvertedTime in that fileif strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter))% when i-th RR time is equal to annotation time of current epochrrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;elseif ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter)) ...&& ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))% when i-th RR time is not equal to annotation time of current epoch% and i-th RR time is not equal to annotation time of the next epochwhile ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))epochCounter = epochCounter + 1;endrrCounter=1;epochCounter=epochCounter+1;rrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;elseif ~strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter)) ...&& strcmp(rrConvertedTime(i), anTimeGeneratedCell(epochCounter+1))% when i-th RR time is not equal to annotation time of current epoch% and i-th RR time is equal to annotation time of the next epochrrCounter=1;epochCounter=epochCounter+1;rrCollection{epochCounter}(rrCounter) = rrNum(i);rrTimeCollection{epochCounter}(rrCounter) = rrTime(i);rrCounter=rrCounter+1;endend% END OF SYNCHRONIZE RR AND ANNOTATION DATA %% SYNCHRONIZED DATA VALIDITY CHECK 同步數(shù)據(jù)有效性檢驗% generate new annotation matrix, fill the time without annotation 生成新的注釋矩陣,在沒有注釋的情況下填充時間 anClassGeneratedCell = cell(heaTotalEpoch, 1);je = 1;for i=1:heaTotalEpochif strcmp(anTimeGeneratedCell(i), anTime(je))% if the time is the same, copy the annotationanClassGeneratedCell(i) = anClass(je);if je < size(anClass, 1)je = je + 1;endelse% if the time is different, fill with 'none'anClassGeneratedCell(i) = {'none'};endendfprintf('Removing invalid epoch:\n'); isExists = 0;for i=heaTotalEpoch:-1:1flag = 0;if sum(rrCollection{i}) < 28 || sum(rrCollection{i}) > 32% set flag to remove incomplete RR data of that epoch by:% check the sum of RR interval from each epoch,% can't be below 28 or higher than 32% (according to slp04 data, min sum is 29 and max is 30)flag = 1;fprintf(' Epoch %d (time: %s) of %s data is removed because ', ...'incomplete RR data\n', i, anTimeGeneratedCell{i}, fileName);elseif strcmp(anClassGeneratedCell(i), {'none'})% set flag to remove no annotation epochflag = 1;fprintf(' Epoch %d (time: %s) of %s data is removed because ', ...'no annotation\n', i, anTimeGeneratedCell{i}, fileName);elseif strcmp(anClassGeneratedCell(i), {'MT'}) || ...strcmp(anClassGeneratedCell(i), {'M'})% set flag to remove 'MT' or 'M' annotation epochflag = 1;fprintf(' Epoch %d (time: %s) of %s data is removed because ', ...'the annotation is %s\n', i, anTimeGeneratedCell{i}, fileName, ...anClassGeneratedCell{i});end% when the flag is 1, remove the dataif flag == 1anTimeGeneratedCell{i} = [];rrCollection{i} = [];anClassGeneratedCell{i} = [];isExists = 1;endend% print message if no invalid epoch if ~isExistsfprintf('No invalid epoch\n'); end% delete empty row anTimeGeneratedCell = ...anTimeGeneratedCell(~cellfun(@isempty, anTimeGeneratedCell)); rrCollection = rrCollection(~cellfun(@isempty, rrCollection)); anClassGeneratedCell = ...anClassGeneratedCell(~cellfun(@isempty, anClassGeneratedCell));%END OF SYNCHRONIZED DATA VALIDITY CHECK %% PREPARE THE OUTPUT 輸出信息OutputData = struct('filename', fileName, 'time', anTimeGeneratedCell, ...'rr', rrCollection, 'annotation', anClassGeneratedCell, 'age', age, ...'gender', gender, 'weight', weight);% END OF PREPARE THE OUTPUT endloadmatobject函數(shù)
——導(dǎo)入第index例樣本的字段及對應(yīng)數(shù)據(jù)?
function Data = loadmatobject(fileName, index)%Load .mat object by index % Syntax: % Data = loadmatobject(dir, index) % % Input: % *) fileName - file name % *) index - index of .mat's variable to be returned % % Output: % *) Data - index-th variable returnedData = load(fileName);fieldName = fieldnames(Data); %返回結(jié)構(gòu)體的字段名稱Data = Data.(fieldName{index}); %返回數(shù)據(jù)endextractfeatures函數(shù)
——沒有輸出變量，輸出非歸一化、歸一化的特征文件（'xlsx', 'mat'）
function extractfeatures(SlpdbData, destination, outputFormat)%Extract HRV Features % Syntax: % extractfeatures(SlpdbData, destination, outputFormat) % % Input: % *) SlpdbData - struct generated from importslpdb() function % 從importslpdb()函數(shù)生成的結(jié)構(gòu) % *) destination - directory of the result % *) outputFormat - output format: 'xlsx', 'mat', 'all' % % Output: % No output variables, but there are two files output: % hrv_features_unorm - unnormalized features % hrv_features_norm - normalized features % target - matrix total samples X 6 (1 - 6 classes target)nSamples = size(SlpdbData, 1);nClasses = length(unique([SlpdbData.annotation])); % unique：找出數(shù)據(jù)矩陣中所有不重復(fù)數(shù)，確定分類類數(shù)hrv = zeros(nSamples, 25);target = zeros(nSamples, nClasses);target(:, [1 5]) = NaN;for i=1:nSamplesrr_diff = diff(SlpdbData(i).rr); % diff函數(shù)式用于求導(dǎo)數(shù)和差分hrv(i, 1) = HRVFeature.AVNN(SlpdbData(i).rr);hrv(i, 2) = HRVFeature.SDNN(SlpdbData(i).rr);hrv(i, 3) = HRVFeature.RMSSD(rr_diff);hrv(i, 4) = HRVFeature.SDSD(rr_diff);hrv(i, 5) = HRVFeature.NNx(50, rr_diff);hrv(i, 6) = HRVFeature.PNNx(hrv(i, 5), size(SlpdbData(i).rr, 2));hrv(i, 7) = HRVFeature.HRV_TRIANGULAR_IDX(SlpdbData(i).rr);hrv(i, 8) = HRVFeature.SD1(hrv(i, 4));hrv(i, 9) = HRVFeature.SD2(hrv(i, 2), hrv(i, 4));hrv(i, 10) = HRVFeature.SD1_SD2_RATIO(hrv(i, 8), hrv(i, 9));hrv(i, 11) = HRVFeature.S(hrv(i, 8), hrv(i, 9));[TP,pLF,pHF,LFHFratio,VLF,LF,HF,f,Y,NFFT] = ...HRVFeature.fft_val_fun(SlpdbData(i).rr,2);hrv(i, 12) = TP;hrv(i, 13) = pLF;hrv(i, 14) = pHF;hrv(i, 15) = LFHFratio;hrv(i, 16) = VLF;hrv(i, 17) = LF;hrv(i, 18) = HF;% set class annotationswitch SlpdbData(i).annotationcase '1'target(i,6) = 1; % 標記類別target(i,4) = 1;target(i,3) = 1;target(i,2) = 1;case '2'target(i,6) = 2;target(i,4) = 1;target(i,3) = 1;target(i,2) = 1;case '3'target(i,6) = 3;target(i,4) = 2;target(i,3) = 1;target(i,2) = 1;case '4'target(i,6) = 4;target(i,4) = 2;target(i,3) = 1;target(i,2) = 1;case 'R'target(i,6) = 5;target(i,4) = 3;target(i,3) = 2;target(i,2) = 1;case 'W'target(i,6) = 6;target(i,4) = 4;target(i,3) = 3;target(i,2) = 2;otherwisefprintf('Invalid Annotation');returnendendhrv( :, ~any(hrv,1) ) = [];% create a new dir if not existsdirList = dir;isDirExists = 0;for i=1:length(dir)if dirList(i).isdir && strcmp(dirList(i).name, destination)isDirExists = 1;endendif ~isDirExistsmkdir(destination);end% save the data into destinationhrv_features_unorm = hrv;hrv_features_norm = normalizedata(hrv, -1, 1);if strcmp(outputFormat, 'xlsx') || strcmp(outputFormat, 'all') % strcmp是用于做字符串比較，保存成相應(yīng)的特征文件xlswrite(strcat(destination, 'hrv_features_unorm.xlsx'), ...hrv_features_unorm);xlswrite(strcat(destination, 'hrv_features_norm.xlsx'), ...hrv_features_norm);xlswrite(strcat(destination, 'target.xlsx'), target);endif strcmp(outputFormat, 'mat') || strcmp(outputFormat, 'all')save(strcat(destination, 'hrv_features_unorm.mat'), 'hrv_features_unorm');save(strcat(destination, 'hrv_features_norm.mat'), 'hrv_features_norm');save(strcat(destination, 'target.mat'), 'target');endendgetindexrange函數(shù)
——通過index得到輸入向量的范圍
function range = getindexrange(nSamplesEachData, index)%Get range of inputted vector by index. For example [2 3 4 5] is the %inputted nSamplesEachData and index is 2. Then, the output is [3 4 5]. %Explanation: %The sum of [2 3 4 5] is 14 (there are 14 items). %If index = 1, so the output is [1 2] -> total elements are 2 %If index = 2, so the output is [3 4 5] -> total elements are 3 %If index = 3, so the output is [6 7 8 9] -> total elements are 4 %If index = 4, so the output is [10 11 12 13 14] -> total elements are 5 % Syntax: % range = getindexrange(nSamplesEachData, index) % % Input: % *) nSamplesEachData - total number of data in each index % 每個index中的數(shù)據(jù)總數(shù) % *) index - index to be retrieved % % Output: % *) range - a vector contains ordered number of associated index % 包含相關(guān)index的有序數(shù)量向量if sum(index > length(nSamplesEachData)) >= 1disp('Index limit exceeded');returnendrange = [];for i=1:length(index)if index(i) == 1startNum = 1;endNum = nSamplesEachData(index(i));elsestartNum = sum(nSamplesEachData(1:index(i)-1)) + 1;endNum = startNum + nSamplesEachData(index(i)) - 1;endrange = [range startNum:endNum];end endPSOforELM函數(shù)
——PSO的初始化、結(jié)合ELM的PSO迭代
function [result, startTime, endTime] = PSOforELM(nFeatures, trainingData, ...testingData, PSOSettings)%Running PSO with ELM for feature selection and number of hidden nodes % 使用ELM運行PSO來選擇特征和隱藏節(jié)點的數(shù)量 %optimization % Syntax: % [result, startTime, endTime] = PSOforELM(nFeatures, trainingData, ... % testingData, PSOSettings) % % Input: % *) nFeatures - total number of features to be selected % *) trainingData - training data (Matrix size: total samples X nFeatures) % *) testingData - testing data (Matrix size: total samples X nFeatures) % *) PSOSettings - struct contains PSO parameters, examples: % PSOSettings.MAX_ITERATION = MAX_ITERATION; % PSOSettings.nParticles = 20; % PSOSettings.W = 0.6; % PSOSettings.c1 = 1.2; % PSOSettings.c2 = 1.2; % PSOSettings.Wa = 0.95; % PSOSettings.Wf = 0.05; % % Output: % *) result - struct contains records of PSO ELM result % 此結(jié)構(gòu)包含PSO ELM結(jié)果的記錄 % *) startTime - time when the experiment starts % *) endTime - time when the experiment endsstartTime = clock;%% PSO PARAMETER PREPARATION PSO各參數(shù)定義及初始化% max total bits for hidden nodes nHiddenBits = length(dectobin(size(trainingData, 1))); % dectobin十-二進制轉(zhuǎn)換 populationPosition = rand(PSOSettings.nParticles, nFeatures+nHiddenBits) > 0.5; % 隨機矩陣中的每一個數(shù)與0.5比較，若值小于0.5，populationPosition矩陣中相對應(yīng)的值返回1，否則返回0% 不符合條件的重新更新 for i=1:PSOSettings.nParticleswhile bintodec(populationPosition(i, nFeatures+1:end)) < nFeatures || ...bintodec(populationPosition(i, nFeatures+1:end)) > size(trainingData, 1) || ...sum(populationPosition(i, 1:nFeatures)) == 0populationPosition(i, :) = rand(1, nFeatures+nHiddenBits) > 0.5; end endpopulationVelocity = int64(zeros(PSOSettings.nParticles, 1)); % 初始化為0（十進制）% 定義pBest結(jié)構(gòu)體的字段：position、fitness、trainingAccuracy、testingAccuracy pBest(PSOSettings.nParticles).position = []; pBest(PSOSettings.nParticles).fitness = []; pBest(PSOSettings.nParticles).trainingAccuracy = []; pBest(PSOSettings.nParticles).testingAccuracy = [];% 各字段初始化 for i=1:PSOSettings.nParticlespBest(i).position = false(1, nFeatures+nHiddenBits); % max fitness valuepBest(i).fitness = repmat(-1000000, PSOSettings.nParticles, 1);pBest(i).trainingAccuracy = 0;pBest(i).testingAccuracy = 0; end% 定義gBest結(jié)構(gòu)體的字段：position、fitness、trainingAccuracy、testingAccuracy、fromIteration、fromParticle gBest.position = false(1, nFeatures+nHiddenBits); gBest.fitness = -1000000; % max fitness value all particle all iteration gBest.trainingAccuracy = []; gBest.testingAccuracy = []; gBest.fromIteration = []; gBest.fromParticle = [];% 定義result結(jié)構(gòu)體的字段 result(PSOSettings.MAX_ITERATION+1).iteration = []; result(PSOSettings.MAX_ITERATION+1).populationPosition = []; result(PSOSettings.MAX_ITERATION+1).pBest = []; result(PSOSettings.MAX_ITERATION+1).time = []; result(PSOSettings.MAX_ITERATION+1).trainingAccuracy = []; result(PSOSettings.MAX_ITERATION+1).testingAccuracy = []; result(PSOSettings.MAX_ITERATION+1).model = []; result(PSOSettings.MAX_ITERATION+1).gBest = [];% END OF PSO PARAMETER PREPARATION %% INITIALIZATION STEP 初始化步驟 %fitness function evaluation 適應(yīng)度函數(shù)的評價[modelArr, trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest); % evaluatefitness函數(shù)具體見本節(jié)最后gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, populationFitness, ...populationPosition, gBest, 0); % gbestupdate函數(shù)具體見本節(jié)最后% save initial data 保存原始數(shù)據(jù) result(1).iteration = 0; result(1).populationPosition = populationPosition; result(1).pBest = pBest; result(1).time = timeArr; result(1).trainingAccuracy = trainAccArr; result(1).testingAccuracy = testAccArr;%result(1).model = modelArr; result(1).gBest = gBest;% END OF INITIALIZATION STEP %% PSO ITERATION PSO的迭代for iteration=1:PSOSettings.MAX_ITERATION% Update Velocityr1 = rand();r2 = rand();for i=1:PSOSettings.nParticles% calculate velocity value 計算速度值positionDec = int64(bintodec(populationPosition(i, :)));populationVelocity(i, 1) = PSOSettings.W * populationVelocity(i, 1) + ...PSOSettings.c1 * r1 * (bintodec(pBest(i).position) - positionDec) ...+ PSOSettings.c2 * r2 * (bintodec(gBest.position) - positionDec);% update particle positionnewPosDec = abs(int64(positionDec + populationVelocity(i, 1)));newPosBin = dectobin(newPosDec);% if the total bits is lower than nFeatures + nHiddenBits,add zeros in frontif size(newPosBin, 2) < (nFeatures + nHiddenBits)newPosBin = ...[zeros(1, (nFeatures + nHiddenBits) - size(newPosBin, 2)) ...newPosBin];end% if the number of hidden node is more than the number of samplesif bintodec(newPosBin(1, nFeatures+1:end)) > size(trainingData, 1) ...|| size(newPosBin(1, nFeatures+1:end), 2) > nHiddenBitsnewPosBin = ...[newPosBin(1, 1:nFeatures) dectobin(size(trainingData, 1))];end% if the number of selected features is 0 如果選擇的特征數(shù)量為0while sum(newPosBin(1, 1:nFeatures)) == 0newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;end% set the new value of positionpopulationPosition(i, :) = newPosBin;end% fitness function evaluation[modelArr, trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration+1);% save dataresult(iteration+1).iteration = iteration;result(iteration+1).populationPosition = populationPosition;result(iteration+1).pBest = pBest;result(iteration+1).time = timeArr;result(iteration+1).trainingAccuracy = trainAccArr;result(iteration+1).testingAccuracy = testAccArr;%result(iteration+1).model = modelArr;result(iteration+1).gBest = gBest; end% END OF PSO ITERATION endTime = clock; end?PSOforSVM函數(shù)
function [result, startTime, endTime] = PSOforSVM(nFeatures, trainingData, ...testingData, PSOSettings)%Running PSO with SVM for feature selection % Syntax: % [result, startTime, endTime] = PSOforSVM(nFeatures, trainingData, ... % testingData, PSOSettings) % % Input: % *) nFeatures - total number of features to be selected % *) trainingData - training data (Matrix size: total samples X nFeatures) % *) testingData - testing data (Matrix size: total samples X nFeatures) % *) PSOSettings - struct contains PSO parameters, examples: % PSOSettings.MAX_ITERATION = MAX_ITERATION; % PSOSettings.nParticles = 20; % PSOSettings.W = 0.6; % PSOSettings.c1 = 1.2; % PSOSettings.c2 = 1.2; % PSOSettings.Wa = 0.95; % PSOSettings.Wf = 0.05; % % Output: % *) result - struct contains records of PSO SVM result % *) startTime - time when the experiment starts % *) endTime - time when the experiment endsstartTime = clock;%% PSO PARAMETER PREPARATION populationPosition = rand(PSOSettings.nParticles, nFeatures) > 0.5;for i=1:PSOSettings.nParticleswhile sum(populationPosition(i, 1:nFeatures)) == 0populationPosition(i, :) = rand(1, nFeatures) > 0.5;end endpopulationVelocity = int64(zeros(PSOSettings.nParticles, 1)); % in decimal value% pBest pBest(PSOSettings.nParticles).position = []; pBest(PSOSettings.nParticles).fitness = []; pBest(PSOSettings.nParticles).trainingAccuracy = []; pBest(PSOSettings.nParticles).testingAccuracy = [];for i=1:PSOSettings.nParticlespBest(i).position = false(1, nFeatures);% max fitness valuepBest(i).fitness = repmat(-1000000, PSOSettings.nParticles, 1);pBest(i).trainingAccuracy = 0;pBest(i).testingAccuracy = 0; end% gBest gBest.position = false(1, nFeatures); gBest.fitness = -1000000; % max fitness value all particle all iteration gBest.trainingAccuracy = []; gBest.testingAccuracy = []; gBest.fromIteration = []; gBest.fromParticle = [];% initialize struct data result(PSOSettings.MAX_ITERATION+1).iteration = []; result(PSOSettings.MAX_ITERATION+1).populationPosition = []; result(PSOSettings.MAX_ITERATION+1).pBest = []; result(PSOSettings.MAX_ITERATION+1).time = []; result(PSOSettings.MAX_ITERATION+1).trainingAccuracy = []; result(PSOSettings.MAX_ITERATION+1).testingAccuracy = []; result(PSOSettings.MAX_ITERATION+1).gBest = [];% END OF PSO PARAMETER PREPARATION %% INITIALIZATION STEP%fitness function evaluation [trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, populationFitness, ...populationPosition, gBest, 0);% save initial data result(1).iteration = 0; result(1).populationPosition = populationPosition; result(1).pBest = pBest; result(1).time = timeArr; result(1).trainingAccuracy = trainAccArr; result(1).testingAccuracy = testAccArr; result(1).gBest = gBest;% END OF INITIALIZATION STEP %% PSO ITERATIONfor iteration=1:PSOSettings.MAX_ITERATION% Update Velocityr1 = rand();r2 = rand();for i=1:PSOSettings.nParticles% calculate velocity valuepositionDec = int64(bintodec(populationPosition(i, :)));populationVelocity(i, 1) = PSOSettings.W * populationVelocity(i, 1) + ...PSOSettings.c1 * r1 * (bintodec(pBest(i).position) - positionDec) ...+ PSOSettings.c2 * r2 * (bintodec(gBest.position) - positionDec);% update particle positionnewPosDec = abs(int64(positionDec + populationVelocity(i, 1)));newPosBin = dectobin(newPosDec);% if the total bits is lower than nFeatures, add zeros in frontif size(newPosBin, 2) < nFeaturesnewPosBin = [zeros(1, (nFeatures)-size(newPosBin, 2)) newPosBin];end% if the total bits is higher than nFeatures, get first nFeaturesif size(newPosBin, 2) > nFeaturesnewPosBin = newPosBin(1, 1:nFeatures);end% if the number of selected features is 0while sum(newPosBin(1, 1:nFeatures)) == 0newPosBin(1, 1:nFeatures) = rand(1, nFeatures) > 0.5;end% set the new value of positionpopulationPosition(i, :) = newPosBin;end% fitness function evaluation[trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest);gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration+1);% save dataresult(iteration+1).iteration = iteration;result(iteration+1).populationPosition = populationPosition;result(iteration+1).pBest = pBest;result(iteration+1).time = timeArr;result(iteration+1).trainingAccuracy = trainAccArr;result(iteration+1).testingAccuracy = testAccArr;result(iteration+1).gBest = gBest; end% END OF PSO ITERATION endTime = clock;end?
?
?extractresults函數(shù)
——沒有輸出變量，有Excel輸出文件（每次實驗結(jié)果、所有實驗的最好結(jié)果）、保存gBest隨迭代變化圖片
function extractresults(resultRootFolder, nFeatures, classNum, nExperiments, ...nIterations)%Extract raw results of experiment using PSOELM or PSOSVM method % Syntax: % extractresults(resultRootFolder, nFeatures, classNum, nExperiments, ... % nIterations) % % Input: % *) resultRootFolder - root directory of the results % *) nFeatures - total features % *) classNum - number of class in vector -> [2 3 4 6] % *) nExperiments - total experiments of PSO % *) nIterations - total iterations of each PSO % % Output: % No output variables, there are excel output files: % [method]_[filename]_extracted_result - results of each experiment % [method]_result - best of all experiments % % Parameter Example: % resultRootFolder = 'PSOELM_raw_result'; % nFeatures = 18; % classNum = [2 3 4 6]; % nExperiments = 25; % nIterations = 100;fileNames = {'slp01a' 'slp01b' 'slp02a' 'slp02b' 'slp03' 'slp04' ...'slp14' 'slp16' 'slp32' 'slp37' 'slp41' 'slp45' 'slp48' ...'slp59' 'slp60' 'slp61' 'slp66' 'slp67x'}; method = strsplit(resultRootFolder, '_');method = method{1}; % PSOSVM | PSOELMnClassClassifiers = length(classNum); % 4headerEachExp = [];switch methodcase 'PSOELM'headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...'TestAcc', 'ProcessTime(Sec)', 'HiddenNodes', 'SelectedFeatures'};case 'PSOSVM'headerEachExp = {'Experiment', 'gBestFitness', 'TrainAcc', ...'TestAcc', 'ProcessTime(Sec)', 'SelectedFeatures'};endheaderBestExp = headerEachExp;headerBestExp{1} = 'RecordingName';% write header for main excel result (only 1 excel)for i=1:length(classNum)xlswrite(sprintf('%s/%s_result.xlsx', resultRootFolder, method), ...headerBestExp, sprintf('%d classes', classNum(i)));endfor iFile=1:length(fileNames) % loop for each file% eachFileFolder example: PSOELM_raw_result/PSOELM_slp01a_raw_resulteachFileFolder = sprintf('%s/%s_%s_raw_result', resultRootFolder, ...method, fileNames{iFile});for iClass=1:nClassClassifiers % loop for each class numbermatFileName = sprintf('%s/%s_%s_%dclasses_raw_result.mat', ...eachFileFolder, method, fileNames{iFile}, classNum(iClass));ExperimentResult = loadmatobject(matFileName, 1);temp = zeros(nExperiments, length(headerEachExp)-1);nBits = ...length(ExperimentResult(4). ...iterationResult(end).gBest.position) - nFeatures;gBestParticles = false(nExperiments, nFeatures+nBits);tempCell = cell(nExperiments, 1);% get the last gBest result of each experimentfor iExp=1:nExperimentslastResult = ExperimentResult(iExp).iterationResult(end);temp(iExp, 1) = iExp;temp(iExp, 2) = lastResult.gBest.fitness;temp(iExp, 3) = lastResult.gBest.trainingAccuracy;temp(iExp, 4) = lastResult.gBest.testingAccuracy;temp(iExp, 5) = ...etime(ExperimentResult(iExp).endTime, ...ExperimentResult(iExp).startTime);if strcmp(method, 'PSOELM')temp(iExp, 6) = ...bintodec(lastResult.gBest.position(nFeatures+1:end));endtempCell(iExp, 1) = ...{bintostringorder(lastResult.gBest.position(1, 1:nFeatures))};gBestParticles(iExp, :) = lastResult.gBest.position;endeachFileExcelPath = sprintf('%s/%s_%s_extracted_result.xlsx', ...eachFileFolder, method, fileNames{iFile});xlswrite(eachFileExcelPath, headerEachExp, ...sprintf('%d classes', classNum(iClass)), 'A1');xlswrite(eachFileExcelPath, temp, ...sprintf('%d classes', classNum(iClass)), 'A2');xlswrite(eachFileExcelPath, tempCell, ...sprintf('%d classes', classNum(iClass)), ...sprintf('%s2', getexcelcolumncode(length(headerEachExp))));% get the best experiment of each classificationbestExpIdx = find(temp(:, 2) == max(temp(:, 2)));if length(bestExpIdx) > 1% if have the same gBest fitness value, get the max of testAccbestExpIdx = ...bestExpIdx(temp(bestExpIdx, 4) == max(temp(bestExpIdx, 4)));if length(bestExpIdx) > 1% if have the same testAcc, get the max of trainAccbestExpIdx = bestExpIdx(temp(bestExpIdx, 3) == ...max(temp(bestExpIdx, 3)));if length(bestExpIdx) > 1% if have the same trainAcc,% get the min of selected featuresbestExpIdx = bestExpIdx( ...sum(gBestParticles(bestExpIdx, 1:nFeatures), 2) == ...min(sum(gBestParticles(bestExpIdx, 1:nFeatures), 2)));if length(bestExpIdx) > 1% if have the same selected feature,% check the method usedswitch methodcase 'PSOELM'bestExpIdx = ...bestExpIdx(temp(bestExpIdx, 6) == ...min(temp(bestExpIdx, 6)));if length(bestExpIdx) > 1% if have the same selected feature,% get the firstbestExpIdx = bestExpIdx(1);endcase 'PSOSVM'bestExpIdx = bestExpIdx(1);endendendendend% mark the best indexxlswrite(eachFileExcelPath, {'BEST EXPERIMENT'}, ...sprintf('%d classes', classNum(iClass)), ...sprintf('%s%d', ...getexcelcolumncode(length(headerEachExp)+1), bestExpIdx+1));%{% gather gBest fitness of the best experimentgBest = zeros(nIterations, 1);for iItr=1:nIterationsgBest(iItr) = ...ExperimentResult(bestExpIdx). ...iterationResult(iItr+1).gBest.fitness;end% save graphicsf = figure;plot(1:nIterations, gBest);ylabel('gBest Fitness'); xlabel('Iteration');title(sprintf('[%s] Best Experiment of %s (%d classes)', ...method, fileNames{iFile}, classNum(iClass)));saveas(f, ...sprintf('%s/[%s] Best Experiment of %s (%d classes).png', ...eachFileFolder, method, fileNames{iFile}, classNum(iClass)));close all;%}% save result to main excelswitch methodcase 'PSOELM'xlswrite( ...sprintf('%s/%s_result.xlsx', resultRootFolder, ...method), ...[fileNames(iFile) temp(bestExpIdx, 2) ...temp(bestExpIdx, 3) temp(bestExpIdx, 4) ...temp(bestExpIdx, 5) temp(bestExpIdx, 6) ...tempCell(bestExpIdx)], ...sprintf('%d classes', classNum(iClass)), ...sprintf('A%d', iFile+1));case 'PSOSVM'xlswrite( ...sprintf('%s/%s_result.xlsx', resultRootFolder, ...method), ...[fileNames(iFile) temp(bestExpIdx, 2) ...temp(bestExpIdx, 3) temp(bestExpIdx, 4) ...temp(bestExpIdx, 5) tempCell(bestExpIdx)], ...sprintf('%d classes', classNum(iClass)), ...sprintf('A%d', iFile+1));endendendend?evaluatefitness、gbestupdate函數(shù)
function [modelArr, trainAccArr, testAccArr, timeArr, populationFitness, ...pBest] = evaluatefitness(PSOSettings, nFeatures, trainingData, ...testingData, populationPosition, pBest)modelArr(PSOSettings.nParticles).inputWeight = []; % 定義modelArr結(jié)構(gòu)體字段inputWeight、outputWeightmodelArr(PSOSettings.nParticles).outputWeight = []; trainAccArr = zeros(PSOSettings.nParticles, 1);testAccArr = zeros(PSOSettings.nParticles, 1);timeArr = zeros(PSOSettings.nParticles, 1);populationFitness = zeros(PSOSettings.nParticles, 1);for i=1:PSOSettings.nParticlestic;% TRAININGmaskedTrainingFeature = featuremasking(trainingData, ... % featuremasking函數(shù)具體見最后populationPosition(i, 1:nFeatures)); % remove unselected features% prepare the target data% (example: transformation from 4 into [0 0 0 1 0 0])trainingTarget = full(ind2vec(trainingData(:,end)'))'; [Model, trainAcc] = trainELM(maskedTrainingFeature, trainingTarget, ... % trainELM函數(shù)具體見最后bintodec(populationPosition(i, nFeatures+1:end)));% TESTINGmaskedTestingFeature = featuremasking(testingData, ... % featuremasking函數(shù)具體見最后populationPosition(i, 1:nFeatures)); % remove unselected features% prepare the target data% (example: transformation from 4 into [0 0 0 1 0 0])testingTarget = full(ind2vec(testingData(:,end)'))';testAcc = testELM(maskedTestingFeature, testingTarget, Model); % testELM函數(shù)具體見最后populationFitness(i, 1) = fitness(PSOSettings.Wa, PSOSettings.Wf, ... testAcc, populationPosition(i, 1:nFeatures)); % function fitnessValue = fitness(Wa, Wf, acc, featureMask)% fitnessValue = Wa * acc + Wf * (1 - (sum(featureMask)/length(featureMask)));% end% pBest updateischanged = 0;% 滿足以下任意一條件，更新變量if populationFitness(i, 1) > pBest(i).fitnessischanged = 1;elseif populationFitness(i, 1) == pBest(i).fitnessif pBest(i).testingAccuracy < testAccischanged = 1;elseif pBest(i).trainingAccuracy < trainAccischanged = 1;elseif sum(pBest(i).position(1, 1:nFeatures)) > ...sum(populationPosition(i, 1:nFeatures))ischanged = 1;elseif bintodec(pBest(i).position(1, nFeatures+1:end)) > ...bintodec(populationPosition(i, nFeatures+1:end))ischanged = 1;endendif ischangedpBest(i).fitness = populationFitness(i, 1);pBest(i).position = populationPosition(i, :);pBest(i).trainingAccuracy = trainAcc;pBest(i).testingAccuracy = testAcc;end% end of pBest updatemodelArr(i) = Model;timeArr(i) = toc;trainAccArr(i) = trainAcc;testAccArr(i) = testAcc;end end function gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration)if max(populationFitness) >= gBest.fitnessfound = find(populationFitness == max(populationFitness));if length(found) > 1% if have the same gBest fitness value, get the max of testAcc 如果具有相同的gBest適應(yīng)度值，則獲取testAcc的最大值found = found(testAccArr(found) == max(testAccArr(found)));if length(found) > 1% if have the same testAcc, get the max of trainAcc 如果有相同的testAcc，得到trainAcc的最大值found = found(trainAccArr(found) == max(trainAccArr(found)));if length(found) > 1% if have the same trainAcc, get the min of selected features 如果具有相同的trainAcc，則獲取所選特征的最小值found = ...found(sum(populationPosition(found, 1:nFeatures), 2) ...== min(sum(populationPosition(found, 1:nFeatures), 2)));if length(found) > 1% if have the same selected feature,get the min of hidden nodehn = zeros(length(found), 1);for i=1:length(found)hn(i, 1) = bintodec(populationPosition(found(i), ...nFeatures+1:end));endfound = found(hn == min(hn));if length(found) > 1found = found(1);endendendendendgBest.fitness = populationFitness(found);gBest.position = populationPosition(found, :);gBest.trainingAccuracy = trainAccArr(found);gBest.testingAccuracy = testAccArr(found);gBest.fromIteration = iteration;gBest.fromParticle = found;end end?
function [trainAccArr, testAccArr, timeArr, populationFitness, pBest] = ...evaluatefitness(PSOSettings, nFeatures, trainingData, testingData, ...populationPosition, pBest)trainAccArr = zeros(PSOSettings.nParticles, 1);testAccArr = zeros(PSOSettings.nParticles, 1);timeArr = zeros(PSOSettings.nParticles, 1);populationFitness = zeros(PSOSettings.nParticles, 1);for i=1:PSOSettings.nParticlestic;% TRAININGmaskedTrainingFeature = featuremasking(trainingData, ...populationPosition(i, 1:nFeatures)); % remove unselected featuresModel = trainSVM(maskedTrainingFeature, trainingData(:,end), 'RBF');trainAcc = testSVM(maskedTrainingFeature, trainingData(:,end), Model);% TESTINGmaskedTestingFeature = featuremasking(testingData, ...populationPosition(i, 1:nFeatures)); % remove unselected featurestestAcc = testSVM(maskedTestingFeature, testingData(:,end), Model);populationFitness(i, 1) = fitness(PSOSettings.Wa, PSOSettings.Wf, ...testAcc, populationPosition(i, 1:nFeatures));% pBest updateischanged = 0;if populationFitness(i, 1) > pBest(i).fitnessischanged = 1;elseif populationFitness(i, 1) == pBest(i).fitnessif pBest(i).testingAccuracy < testAccischanged = 1;elseif pBest(i).trainingAccuracy < trainAccischanged = 1;elseif sum(pBest(i).position(1, 1:nFeatures)) > ...sum(populationPosition(i, 1:nFeatures))ischanged = 1;endendif ischangedpBest(i).fitness = populationFitness(i, 1);pBest(i).position = populationPosition(i, :);pBest(i).trainingAccuracy = trainAcc;pBest(i).testingAccuracy = testAcc;end% end of pBest updatetimeArr(i) = toc;trainAccArr(i) = trainAcc;testAccArr(i) = testAcc;end endfunction gBest = gbestupdate(nFeatures, trainAccArr, testAccArr, ...populationFitness, populationPosition, gBest, iteration)if max(populationFitness) >= gBest.fitnessfound = find(populationFitness == max(populationFitness));if length(found) > 1% if have the same gBest fitness value, get the max of testAccfound = found(testAccArr(found) == max(testAccArr(found)));if length(found) > 1% if have the same testAcc, get the max of trainAccfound = found(trainAccArr(found) == max(trainAccArr(found)));if length(found) > 1% if have the same trainAcc, get the min of selected featuresfound = ...found(sum(populationPosition(found, 1:nFeatures), 2) ...== min(sum(populationPosition(found, 1:nFeatures), 2)));if length(found) > 1% if have the same selected feature, get the firstfound = found(1);endendendendgBest.fitness = populationFitness(found);gBest.position = populationPosition(found, :);gBest.trainingAccuracy = trainAccArr(found);gBest.testingAccuracy = testAccArr(found);gBest.fromIteration = iteration;gBest.fromParticle = found;end endfeaturemasking、trainELM、testELM函數(shù)
function maskedFeature = featuremasking(feature, mask) %Retrieve masked features % Syntax: % maskedFeature = featuremasking(feature, mask) % % Input: % *) feature - feature collection % (Matrix size: total samples X total features) % *) mask - logical matrix, 1 means selected, 0 is not selected % (Matrix size: 1 X total features) % Output: % *) maskedFeature - matrix with selected features only % (Matrix size: total samples X total selected featuresmaskedFeature = zeros(size(feature, 1), sum(mask));j = 1;for i=1:sum(mask)if mask(1, i) == 1maskedFeature(:, j) = feature(:, i);j = j + 1;endend endfunction [ELMModel, trainAcc] = trainELM(feature, target, nHiddenNode) %Train Extreme Learning Machine (ELM) model % Syntax: % [ELMModel, trainAcc] = trainELM(feature, target, nHiddenNode) % % Input: % *) feature - feature collection % (Matrix size: total samples X total features) % *) target - target of each sample % (Matrix size: total samples X total classes) % Example: class 4 -> target is [0 0 0 1] % *) nHiddenNode - total hidden nodes of ELM % % Output: % *) ELMModel.inputWeight - input weight of ELM % (Matrix size: nHiddenNode (+1 for bias) X total features) % *) ELMModel.outputWeight - output weight of ELM % (Matrix size: total classes X nHiddenNode) % *) trainAcc - training accuracyif size(feature, 2) == 0fprintf('Someting went wrong, no feature selected.');returnend% STEP 1: RANDOMLY ASSIGN INPUT WEIGHT AND BIAS 隨機分配輸入權(quán)重和偏差minWeight = -1;maxWeight = 1;inputWeight = (maxWeight-minWeight) .* ...rand(nHiddenNode, size(feature, 2)+1) + minWeight;% STEP 2: CALCULATE THE HIDDEN LAYER OUTPUT MATRIX H 計算隱含層輸出矩陣H% linear combination of hidden output 線性組合的隱藏輸出hiddenOutput = (inputWeight(:, 1:end-1) * feature')+ ...repmat(inputWeight(:, end), 1, size(feature, 1));% apply activation function on hidden output 對隱藏輸出應(yīng)用激活函數(shù)hiddenOutput = sigmoid(hiddenOutput); % function y = sigmoid(x)% y = 1./(1 + exp(-1.*x));% end% STEP 3: CALCULATE THE OUTPUT WEIGHT B 計算輸出權(quán)值B% estimate output weightoutputWeight = target' * pinv(hiddenOutput); % pinv函數(shù)-求矩陣的偽逆矩陣% STEP 4: APPLY MODEL TO TRAINING DATA 將模型應(yīng)用訓練數(shù)據(jù)% linear combination of predicted output 預(yù)測輸出的線性組合predictedOutput = outputWeight * hiddenOutput;% apply activation function on predicted output 對預(yù)測輸出應(yīng)用激活函數(shù)predictedOutput = sigmoid(predictedOutput);maxPred = max(predictedOutput); predictedClass = zeros(size(predictedOutput, 2), 1);for i=1:size(predictedOutput, 2)class = find(predictedOutput(:, i) == maxPred(i));predictedClass(i) = class(1, 1);endtrainAcc = sum(predictedClass == vec2ind(target')')/ ... % vec2ind：向量到索引，向量中每一列中的所有元素有且只能有一個為1，向量的列數(shù)即為索引矩陣的列數(shù)size(predictedOutput, 2) * 100;ELMModel.inputWeight = inputWeight;ELMModel.outputWeight = outputWeight; endfunction testAcc = testELM(feature, target, ELMModel) %Test Extreme Learning Machine (ELM) model % Syntax: % testAcc = testELM(feature, target, ELMModel) % % Input: % *) feature - feature collection % (Matrix size: total samples X total features) % *) target - target of each sample % (Matrix size: total samples X total classes) % Example: class 4 -> target is [0 0 0 1] % *) ELMModel - ELMModel generated from trainELM() function % % Output: % *) testAcc - testing accuracy% linear combination of hidden outputhiddenOutput = (ELMModel.inputWeight(:, 1:end-1) * feature')+ ...repmat(ELMModel.inputWeight(:, end), 1, size(feature, 1)); % apply activation function on hidden outputhiddenOutput = sigmoid(hiddenOutput);% linear combination of predicted outputpredictedOutput = ELMModel.outputWeight * hiddenOutput;% apply activation function on predicted outputpredictedOutput = sigmoid(predictedOutput);maxPred = max(predictedOutput);predictedClass = zeros(size(predictedOutput, 2), 1);for i=1:size(predictedOutput, 2)class = find(predictedOutput(:, i) == maxPred(i));predictedClass(i) = class(1, 1);endtestAcc = sum(predictedClass == vec2ind(target')')/ ...size(predictedOutput, 2) * 100; end?
?
總結(jié)
以上是生活随笔為你收集整理的睡眠阶段分期——SVM和ELM分别与粒子群算法结合（function）的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。
- 上一篇： windows修改mac地址
- 下一篇： 参数整定类毕业论文文献有哪些？