Filter löschen
Filter löschen

Normal Distribution generation using acceptance rejection

15 Ansichten (letzte 30 Tage)
Tarek Mahjoub
Tarek Mahjoub am 30 Jul. 2021
Bearbeitet: Yazan am 30 Jul. 2021
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 Kommentar
Tarek Mahjoub
Tarek Mahjoub am 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;

Melden Sie sich an, um zu kommentieren.

Antworten (1)

Yazan
Yazan am 30 Jul. 2021
Bearbeitet: Yazan am 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';

Kategorien

Mehr zu Particle & Nuclear Physics finden Sie in Help Center und File Exchange

Community Treasure Hunt

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

Start Hunting!

Translated by