Normal Distribution generation using acceptance rejection

조회 수: 9 (최근 30일)
Tarek Mahjoub
Tarek Mahjoub 2021년 7월 30일
편집: Yazan 2021년 7월 30일
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
Tarek Mahjoub
Tarek Mahjoub 2021년 7월 30일
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;

댓글을 달려면 로그인하십시오.

답변 (1개)

Yazan
Yazan 2021년 7월 30일
편집: Yazan 2021년 7월 30일
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';

카테고리

Help CenterFile Exchange에서 Exploration and Visualization에 대해 자세히 알아보기

Community Treasure Hunt

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

Start Hunting!

Translated by