模擬分段機率密度函數

 

clc
close
clear
format long
syms x a
p1=int(exp(2*x),-0.5,0);
p2=int(exp(-2*x),0,0.5)

double(p1+p2)

total=1000000;
c=0;
for i=1:total
    rnd=rand();
    if(rnd<=0.5)
        rnd1=0.5*rand();
         v=0.5 *log(2*rnd1); 
    else
         rnd2=0.5*rand();
         v=-0.5 *log(1-2*rnd2); 
    end

         if( v>=-0.5 && v<=0.5)
           c=c+1; 
         end
end
(c/total)

View Code

說明

其中的rnd1能夠用rnd代替，可是rnd2不能用rnd*0.5，這樣的話結果對不上

clc
close
clear
format long
syms t a  b  y c x1 x2 u1 u2 uy E_YX1  E_YX2 E_X1X2 E_X12 E_X22 x


total=100000;
c=0;
for i=1:total
  rnd=rand();
  if(rnd<=1/4)
     v= sqrt(rand()*(1/4) * 4)+2;
  else
     u=rand()*(3/4);
     v=-sqrt(-12*u +9)+6;
  
  end
  if( v>2.5 && v<3.5)
     c=c+1; 
  end
end
c/total

p1=int(0.5*(x-2),2.5,3);
p2=int(0.5*(2-(x/3)),3,3.5);
double(p1+p2)

View Code

這個在開根號時須要考慮3<x<=6時， 0<u<3/4 ,因此開根號後加負號

p=int(4+0*x,0.1,0.2);
double(p)
total=1000000;
c=0;
for i=1:total
    rnd=rand()*(1/4);
    v=rnd;

     if( v>=0.1 && v<=0.2)
       c=c+1; 
     end
end
(c/total)

View Code

X時0-1上的均勻分佈隨機變量兩那麼  Y=(1/4)X就是 0-1/4上的均勻分佈隨機變量，機率密度是 fy(y)=4

若是要6.0-6.25上的均勻分佈隨機變量那麼就令 Y=(1/4)X+6

clc
close
clear
format long
syms t a  b  y c x1 x2 u1 u2 uy E_YX1  E_YX2 E_X1X2 E_X12 E_X22 x


total=1000000;
c=0;

for i=1:total
  rnd=rand();
  if(rnd<=1/2)
     v= 6*(rand()*0.5) -3 ;
  else
     u=rand()*(1/2) +0.5;
     v=4*sqrt(2*u -1);
  
  end
  if( v>-2 && v<2)
     c=c+1; 
  end
end
c/total

p1=(0.5+2^2/32)-(0.5+0^2/32);
p2=(0.5+-0/6) -(0.5+-2/6);
double(p1+p2)

View Code

 

取捨法模擬變量

clc;clear;
syms x

num=1000000
total=0
c=0;
for i=1 : num
    rnd=rand();
    Y=rnd^(1/3);

    if ((1-2*Y+Y^2)>=rand())
       X= Y;
       total=total+1;
       if(X>=0.5 && X<=0.8) 
           c=c+1; 
        end
    end

end
% for i=1 : num
%     U1=rand();
% 
%     Y=U1;
% 
%     if (16*(Y^2-2*Y^3+Y^4)>=rand())
%        X= Y;
%        total=total+1;
%        if(X>=0.5 && X<=0.8) 
%            c=c+1; 
%         end
%     end
% 
% end

c/total
p=int(30*(x^2-2*x^3+x^4),0.5,0.8);
%p=int(3*x^2,0.5,0.8);
double(p)

View Code

採用g(x)=3x^2,c=10