
N_rep=10000;
N_p=1000;

tab0=[3 6 0
    6 20 2
    1 24 1
    1 8 5];
N=sum(sum(tab0));
[R,C]=size(tab0);

R_prior=sum(tab0,2);
R_dist=[];
for i=1:R
    R_dist=[R_dist i*ones(1,R_prior(i))];
end

C_prior=sum(tab0,1);
C_dist=[];
for i=1:C
    C_dist=[C_dist i*ones(1,C_prior(i))];
end

p_dist=zeros(N_p,R,C);
tic
for t=1:N_p
       Ntab=zeros(N_rep,R,C);
    for j=1:N_rep
        tab=zeros(R,C);
        for i = 1:N
            C_class=C_dist(ceil(N*rand));
            R_class=R_dist(ceil(N*rand));
            tab(R_class,C_class)=tab(R_class,C_class)+1;
        end
        Ntab(j,:,:)=tab;
    end
    for k=1:R
        for l=1:C
            p_dist(t,k,l)=length(find(Ntab(:,k,l)>=tab0(k,l)))/N_rep;
        end
    end
   
end
 toc
for k=1:R
    for l=1:C
        subplot(R,C,(k-1)*C+l)
        x=p_dist(:,k,l);
        hist(x,30);
        title(prctile(x,[2.5 50 97.5])')
    end
end
