CODEWITHSUNDEEP
matlabrandomsampling
user131983
user131983

Reputation: 3927

Matlab: Rejection sampling

I want to sample from only the tails ([-5sigma,-3sigma] and [3sigma,5sigma]) of a Normal Distribution when I run a Monte-Carlo Simulation and therefore Rejection Sampling comes to mind. I am however, struggling to implement this in Matlab. Up till now I have been using something similar to the code below (which I know isn't rejection sampling), but would Rejection Sampling be a better way to solve this issue?

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

Thanks

Edit: Using code in the Answer enter image description here

Upvotes: 3

Views: 1729

Answers (1)

Daniel
Daniel

Reputation: 36710

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);

The accuracy is not what I expected, but I'll leave it here. Maybe someone is able to improve the code.

Upvotes: 1

Related Questions