generate_samples.m

一个用EM算法的源程序

groupCount = length(basisDimensions);

if max(basisDimensions) >= ambientSpaceDimension,
    error('for generate_samples.m, group dimensions must correspond to proper subspaces of the ambient space.')
end
if min(basisDimensions) <= 0,
    error('for generate_samples.m, the group dimensions must be greater than zero.')
end

if ~exist('groupSizes','var'),
    dimensionSampleCounts = zeros(1,ambientSpaceDimension-1);
    % Decide how many points groups of various dimensions should have.
    for dimensionIndex = 1:ambientSpaceDimension-1;
        % Try to preserve the arial density over dimensions.
        % This may be tragically incorrect, especialy since the area of the
        % distribution will depend on whether we are 'uniformSphere' or
        % 'uniformCube'
        switch groupDistributionType
            case {'normal' 'uniformCube'}
                dimensionSampleCounts(dimensionIndex) = ceil(BASE_SAMPLE_COUNT * 1);  
            case 'uniformSphere'
                % The number of points placed in the uniform
                % disc(generally, sphere) should be proportional to its
                % volume?
                dimensionSampleCounts(dimensionIndex) = ceil(BASE_SAMPLE_COUNT * hypersphere_volume(dimensionIndex));
            case 'uniformSphereSurface'
                % The number of points points placed one the surface of the
                % sphere should be proportional to its area (i.e. volume of
                dimensionSampleCounts(dimensionIndex) = max(ceil(BASE_SAMPLE_COUNT * hypersphere_area(dimensionIndex)),BASE_SAMPLE_COUNT);
        end
    end
    groupSizes = dimensionSampleCounts(basisDimensions);
end

% Initialization of the group Bases.
groupBases = cell(1,groupCount);

minimumSubspaceAngleViolated = true;
iterationIndex = 0;
newGroupOrientation = cell(1,groupCount);
sampleNumber = 0;
% Generate the data for each group
sampleLabels = zeros(1,sum(groupSizes));
X = zeros(ambientSpaceDimension,sum(groupSizes));
while(minimumSubspaceAngleViolated && iterationIndex <= maxIterations)
    iterationIndex = iterationIndex + 1;

    for groupIndex = 1:groupCount
        
        % Rotate the group to an arbitrary orientation.
        newGroupOrientation{groupIndex} = rand_special_orthogonal(ambientSpaceDimension);
        if ALIGN_FIRST_GROUP && groupIndex==1
            % Without loss of generality, leave the first group aligned with
            % the original axes.
            groupBases{groupIndex} = eye(ambientSpaceDimension, basisDimensions(groupIndex));
        else
            groupBases{groupIndex} = newGroupOrientation{groupIndex}(:,1:basisDimensions(groupIndex));
        end
        
    end %for
    
    % Use the group Bases to determine if the data that was generated
    % satisfies the minimum subspace angle criteria.
    smallestPairwiseAngle = pi/2; 
    for firstGroupIndex = 1:groupCount,
        for secondGroupIndex = firstGroupIndex+1:groupCount,
            firstGroupBases = groupBases{firstGroupIndex};
            secondGroupBases = groupBases{secondGroupIndex};
            currentPairwiseAngle = subspace_angle(firstGroupBases, secondGroupBases);
            if currentPairwiseAngle < smallestPairwiseAngle,
                smallestPairwiseAngle = currentPairwiseAngle;
            end % if     
        end % for
    end % for
    
    if smallestPairwiseAngle >= minimumSubspaceAngle
        % Generate samples according to the bases
        for groupIndex = 1:groupCount
            % Generate the data for the current Group, aligned with the first
            % coordinate axes for now.
            currentGroup = zeros(ambientSpaceDimension, groupSizes(groupIndex));
            switch groupDistributionType
                case 'uniformCube'
                    currentGroup(1:basisDimensions(groupIndex),:) = (rand(basisDimensions(groupIndex), groupSizes(groupIndex)) - .5);
                case 'uniformSphere'
                    currentGroup(1:basisDimensions(groupIndex),:) = rand_uniform_inside_hypersphere(basisDimensions(groupIndex), groupSizes(groupIndex));
                case 'uniformSphereSurface'
                    currentGroup(1:basisDimensions(groupIndex),:) = rand_uniform_on_hypersphere(basisDimensions(groupIndex), groupSizes(groupIndex));
                case 'normal'
                    currentGroup(1:basisDimensions(groupIndex),:) = randn(basisDimensions(groupIndex), groupSizes(groupIndex));
            end

            if ALIGN_FIRST_GROUP
                % Without loss of generality, leave the first group aligned with
                % the original axes.
                if groupIndex ~= 1
                    currentGroup =  newGroupOrientation{groupIndex} * currentGroup;
                end
            else
                currentGroup =  newGroupOrientation{groupIndex} * currentGroup;
            end

            % Combine the data and generate labels for the samples.
            X(:, sampleNumber+1:sampleNumber+groupSizes(groupIndex)) = currentGroup;
            sampleLabels(sampleNumber+1:sampleNumber+groupSizes(groupIndex)) = groupIndex;
            sampleNumber = sampleNumber + groupSizes(groupIndex);
        end
        
        minimumSubspaceAngleViolated = false;
    end % if
    
end %while

if minimumSubspaceAngleViolated
    error(['Unable to find a data set that enforces the minimum angle after ' num2str(maxIterations) ' iterations. Try lowering ''minimumSubspaceAngle''']);
end

% Normalize the sample magnitude
if ~strcmp(groupDistributionType,'normal')
    for sampleIndex=1:sampleNumber
        Xnorm(sampleIndex) = norm(X(:,sampleIndex));
    end
    X=X/max(Xnorm);
end

% Apply the noise function globally to the data.
switch noiseStatistic
    case 'uniform'
        noise = noiseLevel * 2*(rand(size(X)) - .5);
    case 'normal'
        noise = noiseLevel * randn(size(X));
end
switch noiseType
    case 'multiplicative'
        X = X .* noise;
    case 'additive'
        X = X + noise;
end

% Avoid Intersections
if avoidIntersection
    
    % Randomly change the center of the sample cluster on each subspace
    % away from the origin. 
    % Notice: by default the magnitudes range in [-1, 1].
    offsetMagnitude = 5;
    for groupIndex=1:groupCount
        offset = 2*offsetMagnitude*(rand(basisDimensions(groupIndex),1)-0.5*ones(basisDimensions(groupIndex),1));
        subspaceIndices = find(sampleLabels==groupIndex);
        offsetCenter = repmat(groupBases{groupIndex}*offset,1, length(subspaceIndices));
        X(:,subspaceIndices) = X(:,subspaceIndices) + offsetCenter;
    end
else
    offsetMagnitude = 1;
end

% Add outliers
if outlierPercentage>0 || outlierNumber>0
    if outlierNumber==0
        outlierNumber = ceil(sampleNumber * outlierPercentage / (1-outlierPercentage));
    end
    
    if strcmp(groupDistributionType,'normal')
        % the samples have a normal distribution
        for sampleIndex=1:sampleNumber
            Xnorm(sampleIndex) = norm(X(:,sampleIndex));
        end
        outlierMagnitude=max(Xnorm);
    else
        % otherwise, they have been normalized with magnitude max == 1;
        outlierMagnitude = offsetMagnitude;
    end
    outliers = outlierMagnitude*2*(rand(ambientSpaceDimension, outlierNumber)-0.5*ones(ambientSpaceDimension, outlierNumber));
    X = [X outliers];
    sampleLabels = [sampleLabels -ones(1,outlierNumber)];
end

% Scramble the order of the data in the array, scrambling the labels
% correspondingly.
if scrambleOrder,
    dataPermutationVector = randperm(size(X, 2));
    X = X(:,dataPermutationVector);
    sampleLabels = sampleLabels(dataPermutationVector);
end