function [czx, cz, cx, N] = marginalize(dataset, xs, for_cluster)


if isstruct(dataset)

    if ~exist('for_cluster','var')
        for_cluster=-1;
    end

    if nargin == 1
        xs = zeros(size(dataset.measurements));
    end

        xEntries = unique(xs(~isnan(xs)));
        zEntries = unique(dataset.measurements(~isnan(dataset.measurements)));

        czx= zeros(numel(xEntries), numel(zEntries));


if ~isfield(dataset,'flags')
    if for_cluster==-1
        for i=1:numel(zEntries)
            for j=1:numel(xEntries)
                czx(i,j) = sum((xs == xEntries(j)) & (dataset.measurements == zEntries(i)));
            end
        end
    else
        for i=1:numel(zEntries)
            for j=1:numel(xEntries)
                czx(i,j) = sum((dataset.clusters==for_cluster) & (xs == xEntries(j)) & (dataset.measurements == zEntries(i)));
            end
        end

    end

else        
    if for_cluster==-1
        for i=1:numel(zEntries)
            for j=1:numel(xEntries)
                czx(i,j) = sum((dataset.flags==0) & (xs == xEntries(j)) & (dataset.measurements == zEntries(i)));
            end
        end
    else
        for i=1:numel(zEntries)
            for j=1:numel(xEntries)
                czx(i,j) = sum( (dataset.flags==0) & (dataset.clusters==for_cluster) & (xs == xEntries(j)) & (dataset.measurements == zEntries(i)));
            end
        end

    end
end


else
    if nargin == 1
        xs = zeros(size(dataset));
    end

        xEntries = unique(xs(~isnan(xs)));
        zEntries = unique(dataset(~isnan(dataset)));
        czx= zeros(numel(xEntries), numel(zEntries));
        for i=1:numel(zEntries)
            for j=1:numel(xEntries)
                czx(i,j) = sum((xs == xEntries(j)) & (dataset == zEntries(i)));
            end
        end
    
end

    cz = sum(czx')';
    cx = sum(czx)';
    N = sum(czx(:));
end