function [y]=Binomial(m,n,pb,nb) // Binomial law (p,N) // P{X=n} = C_N^n p^n (1-p)^(N-n) //---------------------------------- res=[]; // we use blocks of size 100 to avoid overflows ntir=100; ntirc=ntir; y=rand(ntir,nb,'uniform'); indy= find( y < pb); y=0*ones(y); y(indy)=1; y=sum(y,'c') res=[res;y]; while ( ntirc < m*n ) y=rand(ntir,nb,'uniform'); indy= find(y< pb); y=0*ones(y); y(indy)=1; y=sum(y,'c') res=[res;y]; ntirc=ntirc+ntir; end y=matrix(res(1:m*n),m,n);Simulation code
prb=0.5;N=10; y=Binomial(1,n,prb,N); i=0:10; z=[]; for i1=i, z=[z,prod(size(find(y==i1)))],end plot2d3("onn",i',z'/n,[1,3],"161","Simulation"); //Theoritical curve i=0:N; zt=[]; deff('[y]=fact(n)','y=prod(1:n)'); deff('[z]=C(N,n)','z= fact(N)/(fact(n)*fact(N-n))'); for j=i, zt=[zt, C(N,j)*prb^j*(1-prb)^(N-j)];end plot2d1("onn",i',zt',[-2,6],"100","Theory");