首页/文章/ 详情

重要性抽样方法实例分享

过冷水

3年前浏览3013

过冷水诚挚邀请你加入Matlab仿真秀官方交流群进行Matlab学习、问题咨询、 Matlab相关资料下载，群号：927550334

经过连续不断的推送Monte Carlo方法，所以我们对其了解透彻了吗？NO！当然还得日日精进,大家经常使用的Monte Carlo方法并不完美，我估计大多数人也听不懂我在说什么，是因为你不知道错在哪了。

对该图像用Monte Carlo求积分

clear all

warning off

feature jit off

n=1000;

x=1 9*rand(1,n);

y=50*( 0.9664-1.22*exp(-0.1811*x.^1.5).*sin(-2.805*x 4.49)) 2.25;

I=9*sum(y)/n

I =

449.9730

y2=50*( 0.9664-1.22*exp(-0.1811*x.^1.5).*sin(-2.805*x 4.49)) 2.25;

N=vpa(int(y2,x,1,10))

N =

451.21022742407789201299004991588

比较可知Monte Carlo在抽样只有1000次的时候有明显的误差，我们不能总归结于通过抽样次数的增加来增加计算精度，应该考虑一下其他可能性。过冷水以前关于Monte Carlo方法求定积分问题没有在随机数的抽样上下功夫，之前都是在积分域内均匀随机抽样，称为直接抽样法。直接抽样法完全不考虑被积函数的特点。所以，当被积函数f(x)在积分区域内起伏很大的话，直接抽样法在函数峰值左右取到的样本数目相对偏少，于是求积分的误差就很大，反之，如果所有抽样点的函数值都很接近，直接抽样法的精度就很高。所以对于提高抽样效率来说对于积分值贡献很大的区域抽样要多取些，被积函数值接近于0的区域可以少取些点。这就是重要抽样法，也就是要对函数的分布情况改变抽样的分布。

随机抽样法：若随机变量X的累计分布函数G(x)连续，则随机数r=G(x)在区间[0,1]内均匀分布将等式两边均被G^-1(.)，作用得到G^-1(r)=x,连续，可见以上定理提供了连续型随机数的生成办法。

步骤1：由分布的概率密度分布函数g(x)的积分

得到累计密度分布函数。

步骤2：令G(x)=r,然后从反函数求得G^-1(r)=x，该x的取值就能够符合概率密度分布函数g(x)

通俗的讲就是对原分布的y值进行均匀抽样，则x就是非均匀抽样。大致思路为：

x取均匀值，y为非均匀分布。

对y进行均匀抽样，这x就是非均匀抽样点。该思路就是这样简单，我们具体怎样在实际中使用呢？

假设积分函数为：

把积分函数的密度函数看做，则:

原来是对x进行均匀抽样，现在是对G(x)进行均匀抽样，反推x的抽样

代码为：

syms x x3

n=500;

X1=rand(1,n);

y=x.^2.*sin(x)/2;

y1=int(sin(x)/2,0,x);

g=finverse(y1);

x1=arrayfun(@(x)(acos(1 - 2*x)),X1);

y1=x1.^2.*sin(x1)/2;

I1=(pi)*sum(y1)/n

I1 =

3.5202

n=500;

x2=pi*rand(1,n);

y2=x2.^2.*sin(x2)/2;

I2=(pi)*sum(y2)/n

I2 =

2.7474

I3=vpa(int(x3.^2.*sin(x3)/2,0,pi))

I3 =

由程序可知当直接对x均匀抽样的结果比对sin(x)/2均匀抽样反推x的结果要好！这是因为过冷水只是为了给大家展示案例，实际选取的抽样函数并不合理，如何选取均匀抽样的函数这个要根据实际函数进行分析，在此过冷水了解的不多，就不做拿来主义了，有需要了解的读者可以查看相关资料。

期望估计值法：期望值估计法的原理就是数学中的变量替换。

变换后的等式右边g(x)dx是一个新的分布，g(x)为分布的密度函数，原来要求对dx均匀抽样，然后对抽样点出的f(x)求值，现在我们变成了对g(x)dx均匀抽样，对抽样点处的f(x)/g(x)求值。此时就相当于对dx不均匀抽样，即对这些不均匀分布的抽样点上的f(x)求值。

通俗的讲就是对图像上累计概率密度进行均匀抽样，然后求对应的x值，再用x进行大数定理的计算。

计算步骤为：

（1）根据分布密度函数g(x)产生随机点；

（2）求出抽样点x处的f(x)/g(x)函数值，积分

源代码：

clear all

warning off

feature jit off

syms x

n=900;

x1=rand(1,n);

y1= 0.1128*x-0.1269;%拟合积分图像

g=finverse(y1);

x2=arrayfun(@(x)((1250*x)/141 9/8),x1);

y2=50*( 0.9664-1.22*exp(-0.1811*x2.^1.5).*sin(-2.805*x2 4.49)) 2.25;

I=9*sum(y2)/n

I =

449.9730

本期关于重要性抽样方法的分享就这么多。知识是逐渐积累的过程，过冷水最初只知道的用int（）函数求积分，接触到用Monte Carlo求积分，然后又看到用大数定理求积分，最后抽样方法的改变对大数定理、Monte Carlo都有影响，学问做细后发现好多有趣的点，希望过冷水的抛砖引玉能够激起大家探索的兴趣。