n=50000;
R=1;
B=1;
A=2;
d=[1 3 5];
e=ones(n,1)*d;
y=R*(B*(A*rand(n,3)-1).^e)./e;
[h,x]=hist(y,100);
figure(1); subplot(2,2,[1 2]);
plot(0,0,'w',x,h/n,'LineWidth',2); legend('A=2, B=1','\delta=1','\delta=3','\delta=5')
xlabel('\Delta x'); ylabel('p(x)');
A=1.5;
B=-1;
d=3;
y=R*(B*(A*rand(n,1)-1).^d)/d;
[h,x]=hist(y,100);
figure(1); subplot(223);
plot(x,h/n,'LineWidth',2); legend('A = 1.5, B = -1');
xlabel('\Delta x'); ylabel('p(x)');
B=1;
y=R*(B*(A*rand(n,1)-1).^d)/d;
[h,x]=hist(y,100);
figure(1); subplot(224);
plot(x,h/n,'LineWidth',2); legend('A = 1.5, B = +1')
xlabel('\Delta x'); ylabel('p(x)');
