接收-拒绝采样.docx

接收-拒绝采样x=0:6/9999:6;fori=1:10000ifx(i)=1y(i)=0.5*betapdf(x(i),8,5);elseif(4=x(i))(x(i)5)y(i)=0.5*(x(i)-4);elseif(5x(i))(x(i)=6)y(i)=-0.5*(x(i)-6);elsey(i)=0;endendm1=max(y(1:1667));m2=max(y(6667:10000));m=max(m1,m2);g=ones(1,10000);u=rand(1,10000);accept=0;accept1=0;accept2=0;plot(x,y,x,g*m);holdonh=floor(rand(1,length(x))*length(x))+1;forj=1:10000ifh(j)1667ifg(h(j))*m*u(h(j))=y(h(j))accept=accept+1;plot(x(h(j)),g(h(j))*m*u(h(j)))holdonendendif((h(j)6667)(h(j)10000))ifg(h(j))*m*u(h(j))=y(h(j))accept=accept+1;plot(x(h(j)),g(h(j))*m*u(h(j)))holdonendendendrejectratio=1-accept/10000Then,toget1000samplesintwoenvelopes,Idivided1000samplesinto500and500samplesineachenvelope.Finally,Iget508samplesin1500candidatesfromBetafunctionand504samplesin1000candidatesfromtrianglefunction.Therejectionrateis56.67%,whichmeansitrejectsnearly1.76candidatespersample.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,[x2x3],g2*m2,r);holdonplot([1,1],[0,m1],r,[4,4],[0,m2],r,[6,6],[0,m2],r)holdonh=floor(rand(1,length(x1))*length(x1))+1;fori=1:length(x1)ifu1(h(i))*m1*g1(h(i))=y1(h(i))accept1=accept1+1;plot(x1(h(i)),u1(h(i))*m1*g1(h(i)))holdonendendk=floor(rand(1,(length(x2)))*(length(x2)))+1;fori=1:length(x2)ifu2(k(i))*m2*g2(k(i))=y2(k(i))accept2=accept2+1;plot(x2(k(i)),u2(k(i))*m2*g2(k(i)))holdonendendfori=1:length(x3)ifu3(k(i))*m2*g2(k(i))=y3(k(i))accept2=accept2+1;plot(x3(k(i)),u3(k(i))*m2*g2(k(i)))holdonendendaccept=accept1+accept2;rejectionrate=1-accept/2400