%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % Matlab script to perform cross validation on a model generated for % numeric data from agriculture % % purpose: data mining with SVMs/MLPs/RBFs/RegTree % % requires: % - statistics toolbox for cvpartition (from Matlab2008a) % - neural network toolbox % - readColData script (see link below) % - SVMTorch (compiled version, of course) % - data structure used for the results (cv_results) % % new: % - comparison of SVM result vs. NNet result vs. RBF result vs. regtree % result % - fixed SVM and fixed NNET parameters % - included normalization % ----- % Georg Ru{\ss} % russ@iws.cs.uni-magdeburg.de % 2008-10-17 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% %% Preparation steps, workspace % clean workspace clear all; % set clock clock_start = clock; % seed random for reproducible results %rand('seed',1); % change paths to wherever your data is located and readable to matlab % uses script readColData from % http://web.cecs.pdx.edu/~gerry/MATLAB/plotting/loadingPlotData.html#colHeadings [label,id_column,data]=readColData('sorted_all',10,10,1); % convert labels to strings instead of character arrays; transpose label = transpose(cellstr(label)) %% data-specific stuff % generate three data sets for the nnet to play with % one: N1, Yield2003, EM38 -- target: Yield2004 % two: N1, Yield2003, EM38, N2, REIP32 -- target: Yield2004 % three: N1, Yield2003, EM38, N2, REIP32, N3, REIP49 -- target: Yield2004 % works by eliminating the respective columns from the 'data' matrix above set_1 = data; set_1(:,7) = []; set_1(:,5) = []; set_1(:,4) = []; set_1(:,3) = []; set_1(:,2) = []; set_2 = data; set_2(:,7) = []; set_2(:,5) = []; set_2(:,3) = []; % only set_3 is actually used here set_3 = data; set_3(:,7) = []; label(:,1) = [] % delete "ID" label label(:,7) = [] % labels of input columns label_inputs = label(:,1:7); % target label (probably YIELD'year') label_target = label(:,8); Size_set_3 = size(set_3); %% modeling stage % partitioning parameters for cross validation k=10; % k-fold holdout cross validation p = 1/k; % normalize data (requires double transposition of data matrix) [set_3_norm, PS] = mapstd(transpose(set_3)); set_3_norm = transpose(set_3_norm); % where to store the results modelresults = struct('j',[],'mae_svm',[],'rmse_svm',[], 'mae_mlp',[], 'rmse_mlp', []); modelresults.mae_rbf = []; modelresults.rmse_rbf = []; modelresults.mae_regtree = []; modelresults.rmse_regtree = []; j_res=[]; for j = 1:50 % generate data partition (cvpartition from matlab's statistics toolbox) % size of holdout (test) set: 1/k % (yields indices of data set) cvdata = cvpartition(Size_set_3(1,1),'Holdout',p) % structured object to store actual, predicted and error values (for each % model) svmresults = struct('actual',[], 'prediction', [], 'abserr', [], 'squerr',[]); mlpresults = struct('actual',[], 'prediction', [], 'abserr', [], 'squerr',[]); rbfresults = struct('actual',[], 'prediction', [], 'abserr', [], 'squerr',[]); regtreeresults = struct('actual',[], 'prediction', [], 'abserr', [], 'squerr',[]); j_res = vertcat(j_res,j); TrainSet = struct(); TrainSet.all = set_3_norm(cvdata.training,:); TrainSet.size = size(TrainSet.all); TrainSet.examples = TrainSet.all(:,1:TrainSet.size(1,2)-1); TrainSet.targets = TrainSet.all(:,TrainSet.size(1,2)); TestSet = struct(); TestSet.all = set_3_norm(cvdata.test,:); TestSet.size = size(TestSet.all); TestSet.examples = TestSet.all(:,1:TestSet.size(1,2)-1); TestSet.targets = TestSet.all(:,TestSet.size(1,2)); % The regression tree requires (or should have) un-normalized % values for better understanding of the tree. % Hence, the same data set is used, but with % un-normalized values. TrainSet.all_unnorm = set_3(cvdata.training,:); TrainSet.examples_unnorm = TrainSet.all_unnorm(:,1:TrainSet.size(1,2)-1); TrainSet.targets_unnorm = TrainSet.all_unnorm(:,TrainSet.size(1,2)); TestSet.all_unnorm = set_3(cvdata.test,:); TestSet.size = size(TestSet.all); TestSet.examples_unnorm = TestSet.all_unnorm(:,1:TestSet.size(1,2)-1); TestSet.targets_unnorm = TestSet.all_unnorm(:,TestSet.size(1,2)); % -------------------------------------- % begin SVMTorch stuff % write training data to file dlmwrite('svm-train',TrainSet.size,' '); dlmwrite('svm-train',TrainSet.all,'-append','delimiter',' '); % write test data to file dlmwrite('svm-test',TestSet.size,' '); dlmwrite('svm-test',TestSet.all,'-append','delimiter',' '); %run svmtorch as the model % requires writing a script file first ... fid = fopen('svmtorch.sh','w'); fprintf(fid,'%s\n%s%u%s\n','#! /bin/bash','SVMTorch -rm -t 2 -eps 0.2 -std 4 -c ',60,' svm-train svm-model'); fclose(fid); % ... and making it executable ... ! chmod 700 svmtorch.sh % ... before running it: % generate model ! ./svmtorch.sh % test model and store results in file 'svm-results' ! SVMTest -oa svm-results svm-model svm-test % remove model ! rm svm-model % read in model output on test values (from file) svmres_read_norm = dlmread('svm-results'); % un-normalize the SVM results and TestSet.targets (see help mapstd for calculation) svmres_read = (svmres_read_norm - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); test_targets = (TestSet.targets - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); % append values svmresults.prediction = vertcat(svmresults.prediction,svmres_read); %svmresults.actual = vertcat(svmresults.actual,TestSet.targets); svmresults.actual = vertcat(svmresults.actual,test_targets); % end SVMTorch stuff % ---------------------------------------- % ---------------------------------------- % begin MLP stuff mlp = newff(transpose(TrainSet.examples),transpose(TrainSet.targets),[10]); mlp.trainParam.min_grad = 0.001; mlp.trainParam.epochs = 10000; mlp.trainParam.lr = 0.25; mlp.trainParam.showWindow = false; mlp.trainParam.showCommandLine = true; mlp.trainParam.show = 50; [mlptrained, mlptrainrecord] = train(mlp,transpose(TrainSet.examples),transpose(TrainSet.targets)); mlpresults.prediction = transpose(sim(mlptrained,transpose(TestSet.examples))); mlpresults.actual = TestSet.targets; % un-normalize mlp results to original ranges mlpresults.prediction = (mlpresults.prediction - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); mlpresults.actual = (mlpresults.actual - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); % end MLP stuff % ---------------------------------------- % ---------------------------------------- % begin RBF stuff target_error = 0; spread = 1; max_neurons = 70; df = 100; [rbf, tr] = newrb(transpose(TrainSet.examples),transpose(TrainSet.targets), target_error, spread, max_neurons, df); rbfresults.prediction = transpose(sim(rbf,transpose(TestSet.examples))); rbfresults.actual = TestSet.targets; % un-normalize rbf results to original ranges rbfresults.prediction = (rbfresults.prediction - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); rbfresults.actual = (rbfresults.actual - PS.ymean) * PS.xstd(8)/PS.ystd + PS.xmean(8); % end RBF stuff % ---------------------------------------- % ---------------------------------------- % begin regression tree stuff % construct tree tree = classregtree(TrainSet.examples_unnorm, TrainSet.targets_unnorm, 'names', label_inputs); %,y,param1,val1,param2,val2) % calculate prediction regtreeresults.prediction = eval(tree,TestSet.examples_unnorm); regtreeresults.actual = TestSet.targets_unnorm; % end regression tree stuff % ---------------------------------------- % ---------------------------------------- % collect error values % generate error measures and store them inside the respective results struct svmresults.abserr = abs(svmresults.actual - svmresults.prediction); svmresults.squerr = (svmresults.abserr).^2; mlpresults.abserr = abs(mlpresults.actual - mlpresults.prediction); mlpresults.squerr = (mlpresults.abserr).^2; rbfresults.abserr = abs(rbfresults.actual - rbfresults.prediction); rbfresults.squerr = (rbfresults.abserr).^2; regtreeresults.abserr = abs(regtreeresults.actual - regtreeresults.prediction); regtreeresults.squerr = (regtreeresults.abserr).^2; % calculate mae and rmse svmmae = mean(svmresults.abserr); svmrmse = sqrt(mean(svmresults.squerr)); mlpmae = mean(mlpresults.abserr); mlprmse = sqrt(mean(mlpresults.squerr)); rbfmae = mean(rbfresults.abserr); rbfrmse = sqrt(mean(rbfresults.squerr)); regtreemae = mean(regtreeresults.abserr); regtreermse = sqrt(mean(regtreeresults.squerr)); % store mae and rmse modelresults.mae_mlp = vertcat(modelresults.mae_mlp, mlpmae); modelresults.rmse_mlp = vertcat(modelresults.rmse_mlp, mlprmse); modelresults.mae_svm = vertcat(modelresults.mae_svm, svmmae); modelresults.rmse_svm = vertcat(modelresults.rmse_svm, svmrmse); modelresults.mae_rbf = vertcat(modelresults.mae_rbf, rbfmae); modelresults.rmse_rbf = vertcat(modelresults.rmse_rbf, rbfrmse); modelresults.mae_regtree = vertcat(modelresults.mae_regtree, regtreemae); modelresults.rmse_regtree = vertcat(modelresults.rmse_regtree, regtreermse); end %% plot mae and rmse against j (trial run number, simply collected) e = figure; plot(j_res,modelresults.mae_mlp,'-.b'); hold on; plot(j_res,modelresults.mae_svm,'-r'); plot(j_res,modelresults.mae_rbf,'kx-'); plot(j_res,modelresults.mae_regtree,'go--'); legend('mae mlp','mae svm', 'mae rbf', 'mae regtree'); title('mean absolute error of MLP vs SVM vs RBF vs RegTree'); xlabel('trial #'); ylabel('error'); hold on; h = figure; plot(j_res,modelresults.rmse_mlp,'-.b'); hold on; plot(j_res,modelresults.rmse_svm,'-r') plot(j_res,modelresults.rmse_rbf,'kx-') plot(j_res,modelresults.rmse_regtree,'go--'); legend('rmse mlp','rmse svm', 'rmse rbf', 'rmse regtree'); title('root mean squared error of MLP vs SVM vs RBF vs RegTree'); xlabel('trial #'); ylabel('error'); %% get script stats duration = etime(clock, clock_start)