Normal Distribution generation using acceptance rejection

28 views (last 30 days)
Hi guys, I want to implement the following algorithm and plot the pdf graph of X at the end.
1.Y1=−ln(U1) , U1 is a uniform random number
2.Y2=−ln(U2); U2 is a uniform random number
If Y2≥(Y1−1)^2/2, set|Z|=Y1; otherwise return to Step 1.
3. Generate U. Set Z=|Z| if U≤0.5 and Z= -|Z| if U >0.5
.4. Set X=0.5 Z -2
Does anyone have a suggestion on how to do that?
  1 Comment
Tarek Mahjoub
Tarek Mahjoub on 30 Jul 2021
Here is my attempt. Can anyone help me fix this code?
n=100;
Z = zeros(1,n);
absZ = zeros(1,n);
X = zeros(1,n);
function exp
Y1 = -log(rand(1,n));
Y2 = -log(rand(1,n));
if (Y1 -1)^2/2 <= Y2
absZ = Y1;
else
exp
end
end
U = rand(1,n);
if U <= 0.5
Z = absZ;
else
Z = -absZ;
end
X = Z * 0.5 - 2;

Sign in to comment.

Answers (1)

Yazan
Yazan on 30 Jul 2021
Edited: Yazan on 30 Jul 2021
Here is a simple script to estimqte the pdf of X.
clc, clear
% number of iterations, the higher this value is, the better is the pdf
% estimation
itt = 1000;
X = nan(itt, 1);
for n=1:itt
while 1
Y1 = -log(rand(1));
Y2 = -log(rand(1));
if Y2>=(Y1-1).^2/2
absZ = Y1;
break
end
end
U = rand(1);
if U>0.5
Z = -absZ;
else
Z = absZ;
end
X(n) = 0.5*Z -2;
end
% plot pdf
histogram(X, 'Normalization', 'probability');
ax = gca;
ax.Title.String = 'Probability distribution';
ax.XLabel.String = 'X';
ax.YLabel.String = 'Probability';

Categories

Find more on Particle & Nuclear Physics in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!