遗忘海岸

江湖程序员 -Feiph(LM战士)

导航

模拟分段概率密度函数

 
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

posted on 2019-10-20 07:48  遗忘海岸  阅读(685)  评论(0编辑  收藏  举报