接收-拒绝采样
x=0:6/9999:6;
for i=1:10000
if x(i)<=1
y(i)=0.5*betapdf(x(i),8,5);
elseif (4<=x(i))&(x(i)<5)
y(i)=0.5*(x(i)-4);
elseif (56667)&(h(j)<10000))
if g(h(j))*m*u(h(j))<=y(h(j))
accept=accept+1;
plot(x(h(j)),g(h(j))*m*u(h(j)))
hold on
end
end
end
rejectratio=1-accept/10000
Then, to get 1000 samples in two envelopes, I divided 1000 samples into 500 and 500 samples in each envelope. Finally, I get 508 samples in 1500 candidates from Beta function and 504 samples in 1000 candidates from triangle function. The rejection rate is 56.67%, which means it rejects nearly 1.76 candidates per sample.
x1=0:1/1499:1;
y1=0.5*betapdf(x1,8,5);
x2=4:1/499:5;
y2=0.5*(x2-4);
x3=5:1/499:6;
y3=-0.5*(x3-6);
m1=max(y1);
m2=max(max(y2),max(y3));
g1=ones(1,length(x1));
g2=ones(1,(length(x2)+length(x3)));
u1=rand(1,length(x1));
u2=rand(1,length(x2));
u3=rand(1,length(x3));
accept1=0;
accept2=0;
plot(x1,y1,'b',x2,y2,'b',x3,y3,'b',x1,g1*m1,'r',[x2 x3],g2*m2,'r');
hold on
plot([1,1],[0,m1],'r',[4,4],[0,m2],'r',[6,6],[0,m2],'r')
hold on
h=floor(rand(1,length(x1))*length(x1))+1;
for i=1:length(x1)
if u1(h(i))*m1*g1(h(i))<=y1(h(i))
accept1=accept1+1;
plot(x1(h(i)),u1(h(i))*m1*g1(h(i)))
hold on
end
end
k=floor(rand(1,(length(x2)))*(length(x2)))+1;
for i=1:length(x2)
if u2(k(i))*m2*g2(k(i))<=y2(k(i))
accept2=accept2+1;
plot(x2(k(i)),u2(k(i))*m2*g2(k(i)))
hold on
end
end
for i=1:length(x3)
if u3(k(i))*m2*g2(k(i))<=y3(k(i))
accept2=accept2+1;
plot(x3(k(i)),u3(k(i))*m2*g2(k(i)))
hold on
end
end
accept=accept1+accept2;
rejectionrate=1-accept/2400
本文档为【接收-拒绝采样】，请使用软件OFFICE或WPS软件打开。作品中的文字与图均可以修改和编辑，
图片更改请在作品中右键图片并更换，文字修改请直接点击文字进行修改，也可以新增和删除文档中的内容。
该文档来自用户分享，如有侵权行为请发邮件ishare@vip.sina.com联系网站客服，我们会及时删除。
[版权声明] 本站所有资料为用户分享产生，若发现您的权利被侵害，请联系客服邮件isharekefu@iask.cn，我们尽快处理。
本作品所展示的图片、画像、字体、音乐的版权可能需版权方额外授权，请谨慎使用。
网站提供的党政主题相关内容(国旗、国徽、党徽..)目的在于配合国家政策宣传，仅限个人学习分享使用，禁止用于任何广告和商用目的。
浙公网安备 33021202002483