function [y]=Erlang(m,n,pMean,pVariance) //------------------------------- k = int( (pMean * pMean ) / pVariance + 0.5 ); if (k <= 0) then k=1;end a = k / pMean; // we use blocks of size 100 to avoid overflows res=[]; ntir=100; ntirc=ntir; y=rand(ntir,k,'uniform'); y= -log(prod(y,'r'))/a; res=[res;y]; while ( ntirc < m*n ) y=rand(ntir,k,'uniform'); y= -log(prod(y,'r'))/a; res=[res;y]; ntirc=ntirc+ntir; end y=matrix(res(1:m*n),m,n);
Simulation code:
n=10000; y=Erlang(1,n,10,1) histplot(20,y,[1,1],'061');