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