概念学习之1——重点抽样法

最新推荐文章于 2024-01-05 16:37:47 发布

行僧行僧

最新推荐文章于 2024-01-05 16:37:47 发布

阅读量3.7k

点赞数 2

分类专栏：概念学习 Matlab

本文链接：https://blog.csdn.net/qq_42318112/article/details/109077726

版权

Matlab 同时被 2 个专栏收录

8 篇文章 1 订阅

订阅专栏

概念学习

1 篇文章 0 订阅

订阅专栏

文章目录

说明
1、什么是重点抽样法
- 1.1 随机抽样法
1.2 重点抽样法
2、重点抽样法求积分编程（Matlab）

说明

在学习过程中参考了以下文章或书籍：
《统计计算》

1、什么是重点抽样法

要理解重点抽样法得首先了解随机抽样法，因为重点抽样法就是在随机抽样法的基础上优化得到的。

1.1 随机抽样法

随机抽样法也叫蒙特卡罗方法，简单理解就是采用模拟的方法来逼近真实问题的理论答案。对于求积分的问题而言，随机抽样法又可以分为两种类型：

随机投点法
平均值法

以下两个例子，分别对应两种随机抽样法：
例子：

已知一个函数 $f (x)$ ，现在对它求积分，可以把它积分的区域框起来，往其中随机打点，最后统计在函数图形下方的点数 $n$ 与总点数 $M$ 的比，再根据矩形框的面积 $S$ 就可以近似求得积分： $\int ^4_0f(x)\approx \frac{nS}{M}$ 。
对于高维的情况则有类似的： $\int_{a_{d}}^{b_{d}} \cdots \int_{a_{2}}^{b_{2}} \int_{a_{1}}^{b_{1}} h\left(x_{1}, x_{2}, \ldots, x_{d}\right) d x_{1} d x_{2} \cdots d x_{d}= \frac{nS}{M}$

2.依旧是投点，并且假设 $X$ 坐标符合均匀分布，根据期望的概念: $\frac{\int_{0}^{4} f(x)d x}{b-a}$ 以及大数定律，可以推导得到这样的积分近似计算方法： $\int ^4_0f(x)\approx\frac{4-0}{M} \sum_{i=1}^{M} f\left(X_{i}\right)$
同样对于高维情况也有类似的：
$\int_{a_{d}}^{b_{d}} \cdots \int_{a_{2}}^{b_{2}} \int_{a_{1}}^{b_{1}} h\left(x_{1}, x_{2}, \ldots, x_{d}\right) d x_{1} d x_{2} \cdots d x_{d}=\prod_{j=1}^{d}\left(b_{j}-a_{j}\right) \cdot \frac{1}{M} \sum_{i=1}^{M} h\left(X_{i}\right)$

1.2 重点抽样法

所谓的重点抽样法，拿这个积分的例子来说，从上图可以看出，在0附近，函数值非常接近0，因此在0附近点落在函数以下的概率很小，大部分都落在了函数上面，从积分的公式可以看出 $n$ 越具体，求得的结果越精确，0附近显然不利于让 $n$ 更具体，于是可以让随机投点的重点放在4附近，这样结果就会更精确。

基于平均值法，要想达到这种效果，可以构造一个投点的密度函数 $g (x)$ ，使得在被积分函数在0附近时投点的概率也比较小，这样就做到了有重点的投点模拟的效果，然后构造一个新的函数: $\eta_{i}=\frac{f\left(\boldsymbol{X}_{i}\right)}{g\left(\boldsymbol{X}_{i}\right)}, i=1,2, \ldots, N$

再利用期望的概念就可以得到： $\eta_{1}=\int_{C} \frac{f(\boldsymbol{x})}{g(\boldsymbol{x})} g(\boldsymbol{x}) d \boldsymbol{x}=\int_{C} f(\boldsymbol{x}) d \boldsymbol{x}$

这样一来只需要对 $\eta$ 求期望，就相当于在求积分了，而在模拟的过程中也有了重点。

2、重点抽样法求积分编程（Matlab）

被积分函数为： $f(x)=x^3$
辅助函数为（重点抽样法中投点的密度函数）： $f(x)=\frac{3x^2}{64}$

% 定义被积函数
f =@(x) x.^3;
% 定义辅助函数（重点抽样法中投点的密度函数）
g =@(x) 3.*x.*x./64.*(x>=0&x<=4);
% 按照辅助函数生成随机数
xxx = slicesample(0.5,1000,'pdf',g)';
% 生成图形中方框的坐标
xs = [repelem(0,100)' linspace(0,4,100)' repelem(4,100)' linspace(4,0,100)'];
ys = [linspace(0,64,100)' repelem(64,100)' linspace(64,0,100)' repelem(0,100)'];
% 随机生成横纵坐标
xx = rand(1,1000)*4;
yy = rand(1,1000)*64;
% 求位于曲线下发的点的个数
zz = yy<f(xx);
times = sum(zz);
% 直接积分得到的结果
resFun1 = integral(f,0,4)
% 随机投点法积分结果
s = times/1000*4*64
% 平均值法积分结果
sAve = (4-0)/1000*sum(f(xx))
% 重点抽样法积分结果
sPUD = mean(f(xxx)./g(xxx))
% 接下来是绘图
myFigure = figure(1);
clf(myFigure)
myAxes = axes("Parent",myFigure);
myLine1 = line(myAxes,0:0.01:4,f(0:0.01:4));
myLine1.Color = 'r';
myLine1.LineStyle = '-';
myLine1.LineWidth = 0.5;
myLine2 = line(myAxes,xs,ys, ...
    'color','b','linestyle','-','linewidth',0.5);
myDots = line(myAxes,xx,yy);
myDots.Color = 'k';
myDots.LineStyle = "none";
myDots.Marker = ".";
myAxes.FontName = "宋体";
myAxes.FontSize = 12;
myAxes.XLim =  [-0.2 4.2];
myAxes.YLim =  [0 70];
myAxes.XLabel.String = "x";
myAxes.YLabel.String = "f(x)";
myAxes.Box = "on";
myAxes.LineWidth = 0.5;

% 接下来是计算结果
resFun1 = 64.0000
s = 62.4640
sAve = 63.3940
sPUD = 64.0965