Matlab：拒绝抽样

Question

当我运行蒙特卡罗模拟时，我希望只从正态分布的尾部（[-5sigma，-3sigma]和[3sigma，5sigma]）中抽样，因此想到了拒绝抽样。然而，我在Matlab中实现这一点时遇到了困难。到目前为止，我一直在使用类似下面代码的东西（我知道这不是拒绝抽样），但是拒绝抽样是否是解决这个问题的更好方法？

function [new_E11] = elasticmodulusrng()

new_E11 = normrnd(136e9,9.067e9,[1 1]);

while new_E11>=136e9-3*9.067e9 && new_E11<=136e9+3*9.067e9
        new_E11 = normrnd(136e9,9.067e9,[1 1]);
end

感谢您的提问。以下是关于使用答案中的代码的说明：

编辑：使用答案中的代码

- user131983

---

以下代码是拒绝抽样，只是针对不同的条件。顺便说一句：避免使用魔法数字。将您的常量（136e9和9.067e9）分配给变量，这使得代码更易于阅读和维护。 - Daniel

@Daniel 谢谢。有没有一种方法可以增加从尾部选择数字的概率，而不仅仅是拒绝不来自尾部的样本？ - user131983

@user131983 我认为没有一种简单的方法可以做到这一点。 - Luis Mendo

你的分布的累积分布函数已知，normcdf 在 [-5sigma,-3sigma] 和 [3sigma,5sigma] 之间的部分被归一化为1的积分。符号工具箱可用吗？ - Daniel

@Daniel 是的，我认为这是个好主意。将尾部归一化，使它们的面积等于1。不幸的是，我没有符号工具箱。 - user131983

Answer 1

mu=0
sigma=1
%get scaling factor
scale=(normcdf(5*sigma+mu,mu,sigma)-normcdf(3*sigma+mu,mu,sigma))*2
%define pdf
cpdf=@(x)((1/scale)*normpdf(x,mu,sigma).*((abs(x-mu)<5.*sigma)&(abs(x-mu.*sigma)>3)))
%get cdf via integral
ccdf=@(x)integral(cpdf,mu-5.*sigma,x)
%allow vector inputs
ccdf=@(x)arrayfun(ccdf,x)

%inverse cdf
icdf=@(y)(fzero(@(x)(ccdf(x)-y),.5));
%allow vector inputs    
icdf=@(x)arrayfun(icdf,x);
%icdf is very slow, thus evaluate some numbers and use the cached and interpolated version:
cachedicdf=nan(n+1,1);
x=0:0.01:1;
y=icdf(x);
icdf=@(uni)interp1(x,y,uni);
%plot some example data
hist(icdf(rand(10000000,1)),1000);

准确度不如我预期，但我会把它留在这里。也许有人能够改进代码。

- Daniel

---

如果您有反向累積分布函數（稱為icdf），則icdf(rand(n,m))是您分佈的隨機樣本。 - Daniel

再次感谢。只是确认一下，这对正态分布的另一个区间（-3sigma, 3sigma）没有影响吗？仅仅是尾部的面积现在为1，而其他区间的面积仍然为0.682左右？ - user131983

再次感谢。我只是想知道您是否认为我可以在从（-3 * sigma，3 * sigma）选择数据样本的同时使用它？这样，选择尾部的概率将更高，这也是最初的目的。 - user131983

对于(-3*sigma，3*sigma)的情况，我会使用细化方法，因为只有很少的样本被拒绝。 - Daniel

谢谢。我很抱歉，我有点傻，但是你所说的“thinning”是什么意思？此外，根据你的回答，我不能只使用代码的这一部分吗：x=0:0.01:1; y=icdf(x); icdf=@(uni)interp1(x,y,uni); hist(icdf(rand(10000000,1)),1000);它似乎可以完成工作。谢谢。 - user131983