Simulation of a geometric random variable

Function code:

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