Oct 192010
 

有了前面的一维积分的基础,我们可以开始研究蒙特卡洛模拟方法计算二维积分。其原理也和一维积分差不多,具体步骤是生成分别服从g1(), g2() 概率密度函数的随机数x1, x2, 然后计算 f(x1,x2)/(g1(x1)*g2(x2)) ,不断重复此过程,最后将得到的这些值加起来取平均值就行了。

我们看一个简单的有限积分上下限的积分例子,计算下列函数的积分:

\int_1^2\int_1^3 exp(-x^2-y^2) dx dy

这种例子的积分区域是一个有限的矩形区域,所以直接生成[1,2], [1,3] 区域内的均匀分布随机数,这里需要注意这两个随机数的概率密度函数分别是 g1=1, g2=1/2; 于是可以用下列的代码计算:

>> x1=rand(1e7,1)+1; %[1,2]上的均匀分布随机数
>> x2=rand(1e7,1)*2+1; %[1,3]上均匀分布随机数
>> temp=exp(-x1.^2-x2.^2)/(1*1/2); %计算f(x1,x2)/(g1(x1)*g2(x2))
%在Matlab代码中,我们用x1 代表x, 用x2代表y
>> mean(temp)
ans =
        0.0189168353178035

然后看无穷限的二维积分,其实也不难,就和一维的情况同样处理:分别生成x1和x2都服从标准正态分布,用N() 来表示标准正态分布的概率密度函数 (即g1(), g2() 都是N()) :

N=@(x)(exp(-x.^2/2)/sqrt(2*pi));

最后计算 f(x1,x2)/(N(x1)*N(x2)) 并求平均值就可以了。例如下面这个例子:

\int_{-\infty}^{\infty}\int_{-\infty}^{\infty} exp(-x^2-y^2) dx dy

注意:这个积分的值就是\pi

>> x1=randn(1e6,1); %标准正态分布
>> x2=randn(1e6,1); %标准正态分布
>> temp=exp(-x1.^2-x2.^2)./(N(x1).*N(x2)); 
%计算f(x1,x2)/(g1(x1)*g2(x2)),
%此处g1, g2都是标准正态分布的概率密度函数,
%而且由于x1和x2独立,所以直接相乘即可。
>> mean(temp)
ans =
          3.13962663007921

由上可见计算上下限都是无穷的积分操作简单,接下来我们要解决上下限并非无穷限的二维积分——此时又分两种情况,一种是积分区域还是矩形形状,另一种是积分区域不是矩形形状(即积分上下限中含有积分变量)。这两种情况的解决方法差不多,只是后一种麻烦一些。具体思路是先生成无穷区域内的标准随机变量,然后使得位于积分区域之外随机数对 (x1,x2) 的 f(x1,x2)=0, 区域内的 f(x1,x2) 照常计算,最后将所有的 f(x1,x2)/(g1(x1)*g2(x2)) 相加取平均值即可。 此处的判断 (x1, x2) 是否在积分区域内是技术难度最高的地方, 如果对Matlab很熟练,推荐用基础班讲述的“逻辑寻址”方法做,那样代码少,计算迅速;如果不熟练,也可以用循环+判断语句做,但是那样速度慢。本文的例子全部是使用逻辑寻址的办法处理。

我们先看第一种的例子:

计算:\int_{-\infty}^{0}\int_{0}^{\infty} exp(-x^2-y^2) dx dy

注意:此处选的积分限都是0的好处是,这个函数在四个象限的积分是pi,而此时的积分区域是四象限中的一个象限,又由于函数对称性,可知这个积分的结果是pi/4。

>> x1=randn(1e6,1); %标准正态分布
>> x2=randn(1e6,1); %标准正态分布
>> idx= (x1< =0 & x2>=0); %做一个逻辑筛子,
%筛子中为1的位置是位于积分区域内的 (x1,x2) 对。
>> temp=zeros(size(x1)); 	%预先定义用于存储f(x1,x2)/(g1(x1)*g2(x2))的份额变量 
>> temp(idx)=exp(-x1(idx).^2-x2(idx).^2)./(N(x1(idx)).*N(x2(idx))); 
%计算f(x1,x2)/(g1(x1)*g2(x2)),
%此处g1, g2都是标准正态分布的概率密度函数,
%而且由于x1和x2独立,所以直接相乘即可。
>> mean(temp)
ans =
         0.785398163397448

然后看另外一个例子,此处的积分限中包含了其中一个积分变量——这意味着积分区域变得不规则,至于此例中积分区域到底什么形状,感兴趣的话可以自己用Matlab画个图。

计算:\int_{-\infty}^{0}\int_{y^2}^{\infty} exp(-x^2-y^2) dx dy

这里需要特别注意搞清楚: x的积分区间是 [y^2, +无穷], y的积分区间是 [-无穷, 0],

>> x1=randn(1e6,1); %标准正态分布
>> x2=randn(1e6,1); %标准正态分布
>> idx= (x2< =0 & x1>=x2.^2); %做一个逻辑筛子,
%筛子中为1的位置是位于积分区域内的 (x1,x2) 对。
>> temp=zeros(size(x1)); 	%预先定义用于存储f(x1,x2)/(g1(x1)*g2(x2))的份额变量 
>> temp(idx)=exp(-x1(idx).^2-x2(idx).^2)./(N(x1(idx)).*N(x2(idx))); 
%计算f(x1,x2)/(g1(x1)*g2(x2)),
%此处g1, g2都是标准正态分布的概率密度函数,
%而且由于x1和x2独立,所以直接相乘即可。
>> mean(temp)
ans =
         0.497240279746054

No related posts.

Bookmark/FavoritesSina WeiboGoogle+FacebookQQTwitterYahoo BookmarksBaiduDiggEmailGoogle GmailOutlook.comEvernotePrintAIMLinkedInBlogger PostKindle ItOrkutShare

Related Posts:

  17 Responses to “二维积分-蒙特卡洛模拟法”

  1. 老师您好,我想要一份gibbs抽样程序

  2. 老师您好,我有个问题想请教,就是在利用mcmc方法计算二重积分的时候,如何实现批量计算呢?
    我遇到的问题是
    n=10
    x11=[1,2,3,4,5,6,7,8,9,10]
    x12=[11,12,13,14,15,16,17,18,19,20]

    p1=rand(1e3,1)
    p2=rand(1e3,1)
    idx=(p1+p2<1)

    temp1=zeros(size(p1))
    temp1(idx)=(p1(idx).^x11(1).*p2(idx).^x12(1).*(1-p1(idx)-p2(idx)).^(n-x11(1)-x12(1)))./(1/4)
    m1=mean(temp1)
    如果我下一步想继续算
    temp2=zeros(size(p1))
    temp2(idx)=(p1(idx).^x11(2).*p2(idx).^x12(2).*(1-p1(idx)-p2(idx)).^(n-x11(2)-x12(2)))./(1/4)
    m2=mean(temp2)
    ……
    也就是x11和x12是向量,每对应一组x11和x12,就算一次mean(temp),那该如何实现呢?希望老师帮帮我,我想了一天写了都不对,谢谢老师!等您的消息!

  3. 谢谢老师,最近也在学R,所以老是和matlab混在一起,呵呵,让您见笑了!

  4. 老师,我仿照您之前教的方法计算了一个二维的积分,属于您讲的最后一类(即其中一个积分限中包含令一个积分变量的类型),可是当算到temp(idx)时报错??? Error using ==> mrdivide
    Out of memory. Type HELP MEMORY for your options.
    这样的问题该如何解决呢?
    附我的程序:
    n=50;
    x10=33;
    x01=15;
    k1=0.5;k2=0.5;k3=0.5;k4=0.5;
    m1=x10+k2;
    m2=n-x10-x01+k1+k4;
    m3=x01+k3;
    m4=n-x01+k1+k2+k4;
    c=gamma(n+k1+k2+k3+k4)/(gamma(x10+k2)*gamma(x01+k3)*gamma(n-x10-x01+k1+k4));
    phi=-1/2;
    N=@(x)(exp(-x.^2/2)/sqrt(2*pi));
    % 用mcmc方法计算积分
    u=randn(1e6,1);
    t=rand(1e6,1);
    idx=(u=0);
    temp=zeros(size(u));
    temp(idx)=(t(idx).^(m1-1).*(1-t(idx)).^(m2-1).*(1-t(idx)/2).^(-m4).*u(idx).^(m4-1).*(1-u(idx)).^(m3-1))/((1/2).*N(u(idx)));
    J=mean(temp);
    result=c*2.^(-(x10+k2))*J;result;

    • 刚才附的程序有错,应该是
      n=50;
      x10=33;
      x01=15;
      k1=0.5;k2=0.5;k3=0.5;k4=0.5;
      m1=x10+k2;
      m2=n-x10-x01+k1+k4;
      m3=x01+k3;
      m4=n-x01+k1+k2+k4;
      c=gamma(n+k1+k2+k3+k4)/(gamma(x10+k2)*gamma(x01+k3)*gamma(n-x10-x01+k1+k4));
      phi=-1/2;
      N=@(x)(exp(-x.^2/2)/sqrt(2*pi));
      % 用mcmc方法计算积分
      u=randn(1e6,1);
      t=rand(1e6,1);
      idx=(u=0);
      temp=zeros(size(u));
      temp(idx)=(t(idx).^(m1-1).*(1-t(idx)).^(m2-1).*(1-t(idx)/2).^(-m4).*u(idx).^(m4-1).*(1-u(idx)).^(m3-1))/((1/2).*N(u(idx)));
      J=mean(temp);
      result=c*2.^(-(x10+k2))*J;result;

      • 不知为什么,这一行应该是这样的
        idx=(u=0);
        可是发表后就变成上面那样了

        • 很抱歉,页面显示有问题,我可以发您邮箱请您看一下吧?谢谢!

    • @zuoshanshan, 是你的那句话写错了,你拿出原始的表达式,拆成几个简单的部分分别计算,最后再合起来。至于回复出的问题,我先查找下原因。

      (t(idx).^(m1-1).*(1-t(idx)).^(m2-1).*(1-t(idx)/2).^(-m4).*u(idx).^(m4-1).*(1-u(idx)).^(m3-1))./((1/2).*N(u(idx)));

  5. 老师,对于之前计算复数的正态分布函数值的问题,我根据您的提示想了一下,你觉得这样对吗?
    因为Standard complex normal distribution的密度函数是g(z)=(1/pi)*exp(-|z|^2),这里用到的只是|z|,所以这是个实值,而由于普通正态分布函数的密度函数是f(x)=(1/sqrt(2*pi))*exp(-x^2/2),由于要计算的是标准正态分布函数值,所以上述两个分布的分布函数只相差1/sqrt(pi),所以,比如,要计算z=1-i,那就先去abs(z),然后用normcdf(abs(z))/sqrt(pi),就可算出z的标准正态分布函数值了
    这样到底对不对呀,望请赐教!

    • 我不需要用到复数正态分布,所以就没有兴趣花太多时间研究。大略看了一下介绍,我的想法是,你既然看不太明白Wiki上的用矩阵计算显示的表达式,那么你可以换一个思路——因为一个复数正态分布就是将两个正态分布随机数组合起来——见Wiki中的Definition内第一段。
      所以就直接看着两个正态分布随机数的联合PDF即可。Matlab中也有函数来,参考 mvnpdf() , mvncdf()

  6. 老师,我能不能再问个与这不相关的问题,在matlab中如果要计算复数的标准正态分布函数值该用什么命令啊?我看到只有一个normcdf()可是好像不能放复数进去,我在帮助里找了半天也没有找到,比如要算
    Φ(1-i)这样的东西该怎么办呢?谢谢!

  7. 老师,
    >>temp=zeros(size(x1)); %预先定义用于存储f(x1,x2)/(g1(x1)*g2(x2))的份额变量
    >> temp(idx)=exp(-x1(idx).^2-x2(idx).^2)./(N(x1(idx)).*N(x2(idx)));
    这两个语句特别不理解啊,能不能解释一下,谢谢!

    • temp用来存储模拟值,即: f(x1,x2)/(g1(x1)*g2(x2)) 其中(x1,x2) 是一组随机数, 在程序里面一共生成了 1e6 (10的6次方,具体几次方忘记了) 个这样的随机数对,因此temp一共有1e6个元素
      temp=zeros(size(x1)); 预先令temp所有的元素都是0
      idx 表示所有位于积分域内的 (x1, x2) 的下标
      temp(idx) 表示所有在积分域内的模拟值。
      exp(-x1(idx).^2-x2(idx).^2)./(N(x1(idx)).*N(x2(idx)));就是把积分域内的模拟值计算出来啊。
      原理就是这样,具体的语法在基础班里面有详细讲述,可以把这些语法拆借成基本组成部分逐个分析——如 zeros(size(x1)) 看不懂,就先看 size(x1) 是什么, 再把上一步结果代入 zeros()里面,如此类推

 Leave a Reply

(required)

(required)


*

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <strike> <strong>