Normal Distribution generation using acceptance rejection

6 次查看（过去 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
Ran in:
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 Center 和 File 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