Simulation of a binomial random variable

Function code:

 
 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");