presto 随机抽样

文章目录

• 说明
• 1、什么是重点抽样法

• 1.1 随机抽样法

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

说明

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

1、什么是重点抽样法

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

1.1 随机抽样法

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

1. 随机投点法
2. 平均值法

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

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$坐标符合均匀分布，根据期望的概念:$E[f(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$

再利用期望的概念就可以得到：$E \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

虽然每次计算的结果可能不一样，但是总体来看，采用重点抽样法的结果要更加接近真实积分值。