clear trials=10000; X=randn(trials,1)+4; for i=1:trials while abs(X(i)-4)>2 X(i)=randn(1,1)+4; end end n=(1:100)/10; pX=histc(X,n)/10; % Not really a pdf, just scaled to make the plot look good hold off figure(1) bar(n,pX,'r') y=(n-4).^2+10; hold on plot(n,y,'b') yX=(X-4).^2+10; number_of_bins=100; bin_width=1/10; n=8+(1:number_of_bins)*bin_width; pY=histc(yX,n)/(trials*bin_width); figure(2) bar(n,pY,'r') sum(pY*bin_width) %make sure it's a pdf, must sum to unity