function newData = reduceSampleData(data, dataStatistic, nSample, nBatch)
% reduce given sample data by removing points proportional to statistic
% data

% get data to work on
reducedHcMatrix = data.hcMatrix;
reducedTimeVector = data.timeVector;
assignedBinMatrix = dataStatistic.assignedBinMatrix;
binValueMatrix = dataStatistic.binValueMatrix;
rng(0);
while size(reducedHcMatrix, 2) - nSample > 0
    % compute "probability"
    probabilityVector = ones(1, size(reducedHcMatrix, 2));
    for iDim = 1:size(data.hcMatrix, 1)
        probabilityVector = probabilityVector .* binValueMatrix(iDim, ...
            assignedBinMatrix(iDim, :));
    end
    cumProbabilityVector = [0 cumsum(probabilityVector)];
    totalProbability = cumProbabilityVector(end);
    % select values (this will reduce the total number by 1 to nBatch)
    selectVector = zeros(1,nBatch);
    for iSelect = 1:nBatch
        wheelResult = rand*totalProbability;
        thisFind = find(...
            (cumProbabilityVector(1:(end-1)) <= wheelResult) & ...
            (wheelResult <= cumProbabilityVector(2:end)), 1, 'first');
        selectVector(iSelect) = thisFind;
    end
    selectVector = unique(selectVector);
    
    % remove points
    for iDim = 1:size(data.hcMatrix, 1)
    for iSelect = 1:length(selectVector)
        thisAssignedBin = assignedBinMatrix(iDim, selectVector(iSelect));
        binValueMatrix(iDim, thisAssignedBin) = ...
            binValueMatrix(iDim, thisAssignedBin) - 1;
    end
    end
    reducedHcMatrix(:, selectVector) = [];
    reducedTimeVector(selectVector) = [];
    assignedBinMatrix(:, selectVector) = [];
    
    if (size(reducedHcMatrix, 2) - nSample) < nBatch
        nBatch = size(reducedHcMatrix, 2) - nSample;
    end
end

newData.hcMatrix = reducedHcMatrix;
newData.timeVector = reducedTimeVector;