function [I, dI, invfs, meanIs, stdIs] = miFromCounts_StrongEtAl_fast(czx)
    % Perform subsampling and extrapolation to get o(1/N) correction
    invfs = [1.1 1.2 1.4 1.6 1.8 2];
    nz = size(czx,1);
    nx = size(czx,2);
    for i=1:numel(invfs)
        f_sub = 1/invfs(i);
        I_sub = [];
        M = 400;
        for n=1:M
            czx_sub = floor(f_sub*czx + sqrt(f_sub*czx).*randn(nz,nx));
            I_sub(n) = miFromCounts_naive(czx_sub);
        end
        meanIs(i) = mean(I_sub);
        stdIs(i) = std(I_sub);
    end
    p = polyfit(invfs, meanIs, 1);
    I = p(end);
    
    % Determine variance using Bill's ghetto method
    i = find(invfs == 2);
    dI = stdIs(i)/sqrt(2);
    
end