function [I, dI, invfs, meanIs, stdIs] = miFromCounts_StrongEtAl(czx)
    % Create list from which to subsample
    N = sum(czx(:));
    nz = size(czx,1);
    nx = size(czx,2);
    xs = [];
    zs = [];
    for i=1:nz
        for j=1:nx
            xs = [xs; j*ones(round(czx(i,j)),1)];
            zs = [zs; i*ones(round(czx(i,j)),1)];
        end
    end

    % Perform subsampling and extrapolation to get o(1/N) correction
    invfs = [1.1 1.2 1.4 1.6 1.8 2 2.5 3 3.5 4];
    meanIs = zeros(1,numel(invfs));
    stdIs = zeros(1,numel(invfs));
    for i=1:numel(invfs)
        f_sub = 1/invfs(i);
        M = 100;
        I_sub = zeros(M,1);
        for n=1:M
            N_sub = floor(f_sub*N);
            indices = ceil(N*rand(N_sub,1));
            czx_sub = marginalize(zs(indices), xs(indices));
            I_sub(n) = miFromCounts_naive(czx_sub)/N_sub;
        end
        meanIs(i) = mean(I_sub);
        stdIs(i) = std(I_sub);
    end
    p = polyfit(invfs, meanIs, 1);
    I = p(end)*N;
    
    % Plot extrapolation
    %if false
        figure;
        errorbar(invfs, meanIs, stdIs)
        xl = [0, 4];
        hold on
        plot(xl, xl*p(1)+p(2), 'r');
        set(gca, 'xlim', xl);
    %end
    
    % Determine variance using Bill's ghetto method
    i = find(invfs == 2);
    dI = N*stdIs(i)/sqrt(2);
    
end