【问题标题】:Matlab: Rejection samplingMatlab:拒绝采样
【发布时间】:2014-03-16 13:05:22
【问题描述】:

当我运行蒙特卡罗模拟时,我只想从正态分布的尾部([-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

谢谢

编辑:在答案中使用代码

【问题讨论】:

  • 下面的代码是拒绝抽样,只针对不同的条件。顺便说一句:避免使用幻数。将常量(136e99.067e9)分配给变量,这使得代码更易于阅读和维护。
  • @Daniel 谢谢。有没有办法增加从尾部选择数字的概率,而不是纯粹拒绝不是来自尾部的样本?
  • @user131983 我认为没有理由这样做。至少不是一件容易的事。
  • 你的分布的 cdf 是已知的,normcdf 在 [-5sigma,-3sigma] 和 [3sigma,5sigma] 之间的部分归一化为 1 的整数。符号工具箱可用吗?跨度>
  • @Daniel 是的,我认为这是个好主意。标准化尾部,使其面积等于 1。不幸的是,我没有符号工具箱。

标签: matlab random sampling


【解决方案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);

准确性不是我所期望的,但我会留在这里。也许有人能够改进代码。

【讨论】:

  • 非常感谢。我理解代码的含义,但我只是想知道我应该如何反转 cdf?另外,如果我可能会问,为什么有必要这样做?
  • 如果你有逆 CDF(我们称之为 icdf),icdf(rand(n,m)) 是你分布的随机样本。
  • 再次感谢。只是为了检查一下,这对正态分布的另一个区间(-3*sigma,3*sigma)没有影响吗?只是尾巴的面积现在将是 1,而另一个区间的面积将保持在 0.682 左右?
  • @user131983:我不太明白你的评论。 0.682 来自哪里?我试图完成它,但意识到我不太擅长寻找数值解,精度不是我所期望的。可能是因为我从学生时代就在使用 mupad(符号数学)......
  • 再次感谢。我只是想知道您是否认为我可以在从 (-3*sigma, 3*sigma) 中选择数据样本时使用它?这样一来,从最初目标的尾部中选择的概率就会更高。
猜你喜欢
  • 2018-12-14
  • 2018-06-14
  • 1970-01-01
  • 2013-01-26
  • 2017-02-23
  • 2018-05-08
  • 1970-01-01
  • 2012-11-14
  • 1970-01-01
相关资源
最近更新 更多