pkg load statistics
N = 10000;
x0= 0;sig = 0.1; 

s = normrnd (x0, sig, [N,1]);

bins =linspace(min(s), max(s), 100);
[h, c] =hist(s,bins);
bar (c, h/trapz(bins,h));
hold on;
plot(bins, 1/(2*pi*sig.^2).^0.5*exp(-(bins-x0).^2/(2*sig.^2)));
trapz(bins, 1/(2*pi*sig.^2).^0.5*exp(-(bins-x0).^2/(2*sig.^2)))
xlim([-5, 5]);
hold off;

ans =  0.99985

N = 10000;
xmin=0; xmax =6;
s = random("uniform",xmin, xmax, [N,1]);
bins =linspace(min(s), max(s), 10);
[h, c] =hist(s,bins);
bar (c, h/trapz(bins,h));
hold on;
plot(bins, ones(size(bins))/(xmax-xmin));
hold off;
xlim([xmin-0.1, xmax+0.1]);

N = 10000;
a=2; m= 1;
s = (1+burrrnd(1,1,a,[N,1]))*m;
[h, c]= hist(s,1000);
bar(c, h/trapz(c, h));
xlim([0, 10]);
hold on;
x = linspace(m, 10);
f = a*m^a./(x.^(a+1));
plot(x, f, '--m')
xlim([0, 10])

function y=mypdf(a,m,x)
   y= a*m.^a./(x.^(a+1));
end

a=2; m = 1;

x = linspace(m, 6); p = mypdf(a,m,x); p = p./trapz(x,p); plot(x,p,"-*r");
hold on;
c = cumtrapz(x,p); plot(x, c,'-*m',"linewidth",2); line(xlim,[1 1], "color",'k');

N = 10000;
s = unifrnd (0, 1,[N,1]); 
xi = interp1(c,x,s); plot(xi, s, 'ok',"markersize",2);
[counts,bins] = hist(xi, 100);
bar(bins, counts/trapz(bins, counts));