diff options
Diffstat (limited to 'sourcecodes/parameter_learning/kfoldCrossValid.m')
| -rw-r--r-- | sourcecodes/parameter_learning/kfoldCrossValid.m | 322 |
1 files changed, 322 insertions, 0 deletions
diff --git a/sourcecodes/parameter_learning/kfoldCrossValid.m b/sourcecodes/parameter_learning/kfoldCrossValid.m new file mode 100644 index 00000000..6306c511 --- /dev/null +++ b/sourcecodes/parameter_learning/kfoldCrossValid.m @@ -0,0 +1,322 @@ +function kfoldCrossValid(pre,predict_label,nfolds) +% This function will peform k-fold cross-validation. +% This requires the specification of the name of the variable +% that you want to predict and the number of folds that the +% data should be divided into. + +nfolds = uint8(str2num(nfolds)); + +sfile=strcat(pre,'structure_input.txt'); +dfile=strcat(pre,'continuous_input.txt'); +nnodefile=strcat(pre,'nnode.txt'); +fnnode = fopen(nnodefile,'r'); +nnodes = fscanf(fnnode,'%d'); + +Std_flag=true; +[labels,cases,bnet]=readInput(dfile,sfile,nnodes,Std_flag); + +for i=1:nnodes + if strcmp(labels(i),predict_label) + predict_node = i; + end +end + +predict_cases = bnet.node_sizes(predict_node); + +ncases = size(cases,2); + +kfold_index = kfold_bin(ncases,nfolds); + +if predict_cases == 1 + kfoldCV_continuous(pre,predict_label,nnodes,labels,cases,bnet,predict_node,predict_cases,nfolds,kfold_index); +else + kfoldCV_discrete(pre,predict_label,nnodes,labels,cases,bnet,predict_node,predict_cases,nfolds,kfold_index); +endif + + +end + + +function kfoldCV_continuous(pre,predict_label,nnodes,labels,cases,bnet,predict_node,predict_cases,nfolds,kfold_index) + +ncases = size(cases,2); + +%%Read in original means and standard deviations to report output as +%% untransformed values. +means_orig=cell(1,nnodes); +stdevs_orig=cell(1,nnodes); +labels_orig=cell(1,nnodes); +%Read in original means and standard deviations +mapfile = strcat(pre,'map.txt'); +fmap = fopen(mapfile,'r'); +for i=1:nnodes + buffer = fgetl(mapfile); + temp = cell(1,3); + for j=1:3 + [next,buffer] = strtok(buffer); + temp{j} = next; + end + labels_orig{i} = temp{1}; + means_orig{i} = str2num(temp{3}); + stdevs_orig{i} = str2num(temp{2}); +end +fclose(fmap); +%Need to map the means and stdevs to the correct labels +means = cell(1,nnodes); +stdevs = cell(1,nnodes); +for i = 1:nnodes + for j = 1:nnodes + if strcmp(labels{i},labels_orig{j}) + means{i} = means_orig{j}; + stdevs{i} = stdevs_orig{j}; + break + end + end +end + +kfoldPredictions=zeros(size(cases,2),2); + + +%t=cputime; +%First get predictions +for i =1:nfolds +% i + test_flag = false(ncases,1); + for j=1:ncases + if kfold_index(j) == i + test_flag(j) = true; + end + end + trainData = cases(:,~test_flag); + testData = cases(:,test_flag); + temp_index = zeros(size(testData,2),1); + temp = 1; + for j=1:ncases + if kfold_index(j) == i + temp_index(temp) = j; + temp = temp + 1; + end + end + [bnet] = parameterLearning(bnet,trainData); + for j=1:size(testData,2) + evidence = testData(:,j); + evidence{predict_node} = {}; + engine = jtree_inf_engine(bnet); + [engine,loglik] = enter_evidence(engine,evidence); + predict = marginal_nodes(engine,predict_node); + adj_mu = predict.mu*stdevs{predict_node}+means{predict_node}; + adj_sigma = stdevs{predict_node}*predict.Sigma; + kfoldPredictions(temp_index(j),1) = adj_mu; + kfoldPredictions(temp_index(j),2) = adj_sigma; + end +end +%e=cputime-t; + + +%Open output file. +filename = strcat(pre,'kfoldCV.txt'); +fileID = fopen(filename,'w'); + +fprintf(fileID,'Variable that was predicted: %s\n\n',predict_label); + +%Calculate RMSEP (root mean square error of prediction) and q^2 +%First, calculate TSS (total sum of squares) and +% PRESS (sum of squares of prediction errors) +average = mean(cell2mat(cases'))(predict_node); +%undo standardization +average = average*stdevs{predict_node}+means{predict_node}; +tss = 0; +press = 0; +case_adj=zeros(size(cases,2),1); +for i =1:ncases + case_adj(i) = cases{predict_node,i}*stdevs{predict_node}+means{predict_node}; + tss = (case_adj(i)-average)^2 + tss; + press = (kfoldPredictions(i,1)-case_adj(i))^2 + press; +end +rmsep = sqrt(press/ncases); +q_squared = 1 - press/tss; + +%% Print rmseq and q^2 +fprintf(fileID,'RMS error of predictions: %6.4f\n',rmsep); +fprintf(fileID,'Q^2 of predictions: %6.4f\n\n',q_squared); + + +%%Print the predictions +fprintf(fileID,'Predicted mean and standard deviation for each case:\n'); +fprintf(fileID,'CaseRow\tFoldNumber\tActualValue\tPredictionMean\tPredictionStDev\n'); +for i = 1:ncases + fprintf(fileID,'%i\t%i\t%6.4f\t',i,kfold_index(i),case_adj(i)); + fprintf(fileID,'%6.4f\t%6.4f\n',kfoldPredictions(i,:)); +end + +end + + +function kfoldCV_discrete(pre,predict_label,nnodes,labels,cases,bnet,predict_node,predict_cases,nfolds,kfold_index) + +ncases = size(cases,2); + +%This next section just gets the original names of the levels. +% so they can be written to the output file. +%%Get the maximum_number of states so array will be big enough +%%Add 1 because the input includes the node name +max_states = max(bnet.node_sizes) + 1; +disc_nodes = size(bnet.dnodes,2); + +%%Get mapping of discrete levels. +levelfile = strcat(pre,'nlevels.txt'); +flevels = fopen(levelfile,'r'); +levels = cell(disc_nodes,max_states); +ndisc_nodes = 0; +for i=1:disc_nodes + ndisc_nodes = ndisc_nodes + 1; + buffer = fgetl(flevels); + for j = 1:max_states + [next,buffer] = strtok(buffer); + if j == 1 + levels{i,j} = next; + else + levels{i,j} = next; + end + if length(buffer) < 1 + break + end + end +end + +pred_levels = cell(1,predict_cases); +for i = 1:disc_nodes + if strcmp(levels{i,1},predict_label); + for j = 1:predict_cases + pred_levels{j} = levels{i,j+1}; + end + break + end +end + +kfoldPredictions=zeros(size(cases,2),predict_cases); + +%t=cputime; +%First get kfold CV predictions +for i =1:nfolds + test_flag = false(ncases,1); + for j=1:ncases + if kfold_index(j) == i + test_flag(j) = true; + end + end + trainData = cases(:,~test_flag); + testData = cases(:,test_flag); + temp_index = zeros(size(testData,2),1); + temp = 1; + for j =1:ncases + if kfold_index(j) == i + temp_index(temp) = j; + temp = temp + 1; + end + end + [bnet] = parameterLearning(bnet,trainData); + for j = 1:size(testData,2) + evidence = testData(:,j); + evidence{predict_node} = {}; + engine = jtree_inf_engine(bnet); + [engine,loglik] = enter_evidence(engine,evidence); + predict = marginal_nodes(engine,predict_node); + for k = 1:predict_cases + kfoldPredictions(temp_index(j),k) = predict.T(k); + end + end +end +%e=cputime-t; + +%Now compare with actual outcomes +%actual_states = zeros(1,predict_cases); +%for i=1:ncases +% for j = 1:predict_cases +% if cell2mat(cases(predict_node,i)) == j +% actual_states(j) = actual_states(j) + 1; +% end +% end +%end + +pred_states = zeros(1,ncases); + +for i=1:ncases + max_state = 1; + for j = 2:predict_cases + if kfoldPredictions(i,j) > kfoldPredictions(i,max_state) + max_state = j; + end + end + pred_states(i) = max_state; +end + +correct = 0; +for i=1:ncases + if pred_states(i) == cell2mat(cases(predict_node,i)) + correct = correct + 1; + end +end + +accuracy = correct/ncases; + +%Open output file. +filename = strcat(pre,'kfoldCV.txt'); +fileID = fopen(filename,'w'); + +fprintf(fileID,'Variable that was predicted: %s\n\n',predict_label); + + +%%Print the accuracy +fprintf(fileID,'Fraction of accurate predictions: %6.4f\n\n',accuracy); + +%%Print the predictions +fprintf(fileID,'Predicted likelihood of each state for each case:\n'); +fprintf(fileID,'%s\t%s\t%s\t','CaseRow','FoldNumber','ActualState'); +fprintf(fileID,'%s\t',pred_levels{1:end-1}); +fprintf(fileID,'%s\n',pred_levels{end}); +for i = 1:ncases + fprintf(fileID,'%i\t%i\t',i,kfold_index(i)); + case_level = pred_levels{cases{predict_node,i}}; + fprintf(fileID,'%s\t',case_level); + fprintf(fileID,'%6.4f\t',kfoldPredictions(i,1:end-1)); + fprintf(fileID,'%6.4f\n',kfoldPredictions(i,end)); +end + +end + + +function [kfold_index] = kfold_bin(ncases,nfolds); +%This returns indexes for the different cases to separate data into folds. + +%Get array with random permutation of the number of cases +p = randperm(ncases); + +base_size = idivide(ncases,nfolds); +remainder = rem(ncases,nfolds); + +kfold_index = zeros(ncases,1); +group = 1; +count = 0; +for i=1:ncases + count = count + 1; + kfold_index(p(i)) = group; +%Need to do some checks if groups cannot be exactly equally sized +%If you have already added an extra member to the group, go to next group + if count > base_size + count = 0; + group = group + 1; +%If you have filled the group, check if an extra is needed + elseif count == base_size + if remainder > 0 + remainder = remainder - 1; + else + count = 0; + group = group + 1; + end + end +end + + +end + |