function [y]=Geom(m,n,p) // P(0)= 0 P(i) = p*(1-p)^{n-1} P(inf)=0 // E = 1/p ; sig2= (1-p)/p^2 //-------------------------------------- if p >= 1 then write(%io(2),'p must be < 1');end y=0*ones(m,n) for i=1:m*n, samples=1 z=rand(1,1,'uniform'); while( z < 1-p) ,z=rand(1,1,'uniform'); samples=samples+1;end y(i)= samples; end y=matrix(y,m,n)
Simulation code:
n=1000; pr=0.2 y=Geom(1,n,pr) N=20 i=0:N; z=[]; for i1=i, z=[z,prod(size(find(y==i1)))],end plot2d3("onn",i',z'/n,[1,3],"161","Simulation"); zt=[0];for i1=1:N; zt=[zt,pr*(1-pr)^(i1-1)];end plot2d1("onn",i',zt',[-2,6],"100","Theory");