test_ksubspaces.m

一个用EM算法的源程序

% Complete rewrite of the test_gpca function to use the generalized data
% generation function instead.
clear;
close all;
addpath helper_functions
clc;

roundNumber = 1;

seedTimeStamp = sum(100*clock); %% Generate a new random seed.
rand('state',seedTimeStamp);
randn('state', seedTimeStamp);


% Generate the data samples
ambientDimension = 3;
% Please put the subspace dimensions in descending order to get correct error estimation for this test
% This restriction is not required by the segmentation algorithm.
subspaceDimensions = [2 2 2];
groupCount = length(subspaceDimensions);
noiseLevel = 0.04;

basisErrorStat = 0;
segmentationErrorStat = 0;
for roundIndex=1:roundNumber
    [X, aprioriSampleLabels, aprioriGroupBases] = generate_samples('ambientSpaceDimension', ambientDimension,...
        'groupSizes', 100*subspaceDimensions,...
        'basisDimensions', subspaceDimensions,...
        'noiseLevel', noiseLevel,...
        'outlierNumber',0,...
        'minimumSubspaceAngle', pi/4,...
        'scrambleorder',true);
    sampleNumber = size(X,2);
    
    % Perform Ksubspaces on the data
    [sampleLabels, groupBases] = Ksubspaces(X, subspaceDimensions);


    % Reorder the bases and sample labels so that they are likely to agree with
    % the original ones.
    [bestMapping, errorProbability] = relabel_samples(aprioriSampleLabels, sampleLabels, subspaceDimensions);
    segmentationErrorStat = segmentationErrorStat + errorProbability;
    inverseMapping(bestMapping) = 1:length(bestMapping);
    
    fixedSampleLabels = -ones(1,sampleNumber);
    inlierIndex=find(sampleLabels~=-1);
    fixedSampleLabels(inlierIndex) = bestMapping(sampleLabels(inlierIndex));
    for subspaceIndex=1:groupCount
        fixedBases{subspaceIndex} = groupBases{inverseMapping(subspaceIndex)};
    end
    
    % Calculate the underlying structure error of subspaces by the
    % difference of their space angles.
    basisErrorStat = basisErrorStat + average_basis_error(aprioriGroupBases, fixedBases);
    
    % Display up-to-date error statistics averaged over roundIndex.
    % The plot in 2- or 3-D below can shows that the algorithm may be
    % trapped into local solutions that gives large basis errors.
    disp(['Try #' num2str(roundIndex) '. Segmentation Error: ' num2str(segmentationErrorStat/roundIndex*100) '%.  Average Basis Error: ' num2str(basisErrorStat/roundIndex)]);
end

% Plot the gpca results
if ambientDimension<=3
    figureNumber = figure();
    subplot(1,2,1);
    plot_data(X, aprioriSampleLabels);
    hold on;
    plot_subspaces(X, aprioriSampleLabels, aprioriGroupBases, subspaceDimensions);
    hold off;
    title('Apriori Data');
    set(0,'CurrentFigure',figureNumber); % Change back to the correct figure without changing its visibility.
    subplot(1,2,2);
    plot_data(X, fixedSampleLabels);
    hold on;
    plot_subspaces(X, fixedSampleLabels, fixedBases, subspaceDimensions);
    hold off;
    title('Final Segmentation')
end