N=78;
N_rep=100000;

tab0=[3 6 0
      6 20 2
      1 24 1
      1 8 5];
tab_p=zeros(4,3);
  
Ntab=zeros(N_rep,4,3);
for j=1:N_rep
    k1=0;
    k2=0;
    tab=zeros(4,3);
    for i = 1:N
        
        EME = N*rand(1,1);
        if EME <12,
            k1=1;
        else if EME < 12+58,
                k1=2;
            else
                k1=3;
            end
        end
        SME = N*rand(1,1);
        if SME < 10
            k2=1;
        elseif SME < 38
            k2=2;
        elseif SME < 64
            k2=3;
        else
            k2=4;
        end
    tab(k2,k1)=tab(k2,k1)+1;
    end
    Ntab(j,:,:)=tab;
end

for k=1:4
    for l=1:3
        p_tab(k,l)=length(find(Ntab(:,k,l)>=tab0(k,l)))/N_rep;
    end
end

disp(p_tab)