Simulation of a Poisson random variable

Function code:

 
function [y]=Poisson(m,n,pmean) 
// P{n} = exp(-lambda)lambda^n/n!
// pmean =lambda 
//----------------------------
y=0*ones(m,n)
bound= exp(-pmean);
for i=1:m*n, 
	count=0
	lprod=1
	while( lprod >= bound), lprod=lprod*rand(1,1,'uniform');
	count=count+1;end
	y(i)=count-1;
end
y=matrix(y,m,n)

Simulation code:

n=1000;
pmean=3;
y=Poisson(1,n,pmean);
N=20;
i=0:N;
z=[];
for i1=i, z=[z,prod(size(find(y==i1)))],end
plot2d3("onn",i',z'/n,1,"161");
deff('[y]=fact(n)','if n==0 then y=1;else y=n*fact(n-1);end');
zt=[];for i1=0:N; zt=[zt, exp(-pmean) *pmean^i1/fact(i1)];end
plot2d1("onn",i',zt',[-2,6],"100","Theory");