% Generating NonUniform distribution
% using 3 methods:
% 1. The rejection
% 2. The transform/inverse
% 3. The metropolis
clear  % Clears command window
clc    % Clears command history 
clf    % Removes anything in the figure window before simulation. 

rng('default')
rng('shuffle')      % use seed as current time

% input 1: Max number of points and array of uniform random numbers
Nmax = 100000;
xu = rand(1,Nmax);

% input 2: the area
xmin = 0.0;
xmax = 6.0;
ymin = 0.0;
ymax = 1.0;

% part 1: plot the original function
figure (1)
fplot(@(x) fw(x),[xmin xmax],'b')

% part 2: preparing the histogram - bins
bin = 0.20;     % bin size
nbmax = 100;    % max number of bins
% check for the bin size (change if needed)
nbins = int64((xmax-xmin)/bin);
if nbins > nbmax
    bin = (xmax - xmin)/nbmax;
end
for j = 1:nbmax
    edges(j) =xmin + bin*(j-1);
    if edges(j) >= xmax
        break
    end
end

%Method 1: the rejection method
ki = 0;
for k = 1:Nmax
    x = xmin + (xmax-xmin)*rand;
    y = ymin + (ymax-ymin)*rand;
    if y > fw(x)
        continue
    end
    ki = ki +1;
    xr(ki) = x;
end
fprintf('Succes rate for rejection  %5.3f \n',double(ki)/Nmax)
figure(2)
histogram(xr,edges);
str = sprintf('REJECTION method: N = %7i,  Success rate R = %5.3f ',Nmax,double(ki)/Nmax);
title(str);

% inverse method
for k=1:Nmax
    xi(k) = -log(rand);             % w = exp(-x)
    %xi(k) = acos(2.0*rand-1.0);
end
figure(3)
%histogram(xi)
histogram(xi,edges);
str = sprintf('TRANSFORM method: N = %7i  ',Nmax);
title(str);

% Metropolis method
h = (xmax - xmin)/10.0;     % max step delta in the algorithm
km = 1;
% Choose the first point where w(x) is a maximum
% xm(1) = xmin + (xmax-xmin)*rand;  % middle point
xm(1) = xmin; % left point
for k = 1:Nmax
    xtest = xm(km) + h*(2*rand - 1.0);
    xtest = abs(xtest);
    r = fw(xtest)/fw(xm(km));
    if r < rand
        continue
    else
        km = km + 1;
        xm(km) = xtest;
    end
end
fprintf('Succes rate for Metropolis %5.3f \n',double(km)/Nmax)
figure(4)
%histogram(xm)
histogram(xm,edges);
str = sprintf('METROPOLIS method: N = %7i, success rate R = %5.3f ',Nmax,double(km)/Nmax);
title(str);    


% Probability distribution w(x)
function w = fw(x)
    w = exp(0.0-x);
    %w = 0.5*sin(x);
end