%
% Generating non-uniform distribution using the Rejection method
%
% by AG: January 2022

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: set N number of points to generate
N = 1000000;
fp = @f;        
%input 2: set parameters for the area
xmin = -5.0;
xmax =  5.0;
ymin =  0.0;
ymax = f(0.0);

% part 1: 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

%Preparation (reserve array w)
%w = zeros(N,1);

% part 2: call the rejection method
[w,ki] = Reject(fp,xmin,xmax,ymin,ymax,N);

% part 3: print and plot
fprintf('Success rate for rejection  %5.3f \n',double(ki)/N)
figure(1)
%histogram(w)
histogram(w,edges);
str = sprintf('N = %7i,  Success rate R = %5.3f ',N,double(ki)/N);
title(str);


% the main function - generates a nonuniform distribution
function [wr,ki] = Reject(fp, xmin, xmax, ymin, ymax, N)
ki = 0;
for k = 1:N
    x = xmin + (xmax-xmin)*rand;
    y = ymin + (ymax-ymin)*rand;
    if y > fp(x)
        continue
    end
    ki = ki +1;
    wr(ki) = x;
end
end

% function 1: normal distribution
function y = f(x)
si = 1.0;
mu = 0.0;
y = (1.0/(si*sqrt(2.0*pi)))*exp(1.0-0.5*power((x-mu),2)/power(si,2));
end