Simulation of a disrete random variable

Function code:

function [y]=RndDisc(m,n,x,p)
// discrete law random number
// sum p_i delta_{x_i}
//-------------------------------
p1=[0,p];p1=cumsum(p1);
y=rand(m,n,'uniform')
N=prod(size(x));
res=0*ones(m*n);
for i=1:N,z=0*ones(m*n,1),id=find( p1(i) <= y & y < p1(i+1) ),
 	z(id)=x(i)*ones(prod(size(id))),res=res+z;
end
y=matrix(res,m,n);

Simulation code:

n=10000;
x=[1,3,4,6,10,12];
pr=[0.1,0.2,0.3,0.2,0.1,0.1];
y=RndDisc(1,n,x,pr);
i=0:13
z=[];
for i1=i, z=[z,prod(size(find(y==i1)))],end
plot2d3("onn",i',z'/n,[1,3],"151","Simulation ",[0,0,14,0.5]);
plot2d1("onn",x',pr',[-2,6],"100","Theory");