Draw random numbers from a custom probability density function in Matlab

Question

I want to sample R random numbers from a custom probability density function in Matlab.

This is the expression of the probability density function evaluated at x.

I thought about using slicesample

R=10^6;
f =  @(x) 1/(2*pi^(1/2))*(1/(x^(3/2)))*exp(-1/(4*x));
epsilon= slicesample(0.3,R,'pdf',f,'thin',1,'burnin',1000);

However, it does not work because I get the error

Error using slicesample (line 175)
The step-out procedure failed.

I tried to change starting value and values of thin and burning parameters but it does not seem to work. Could you advise, either on how to make slicesample work or on alternative solutions to sample random numbers from a custom probability density function in Matlab?

Answer 1

Let X be a random variable distributed according to your target pdf. Applying the change of variable y = 1/x and using the well-known theorem for a function of a random variable, the distribution of Y = 1/X is recognized to be a Gamma distribution with parameters α = 1/2, β = 1/4.

Therefore, it suffices to generate a Gamma random variable (using gamrnd) with those parameters and take the inverse. Note that Matlab's definition of the Gamma distribution uses parameters A = α, B = 1/β.

R = 1e5; % desired sample size
x = 1./gamrnd(1/2, 4, [1 R]); % result

Check:

histogram(x, 'Normalization', 'pdf', 'BinEdges', 0:.1:10)
hold on
f = @(x) 1/2/sqrt(pi)./x.^(3/2).*exp(-1/4./x); % target pdf
fplot(f, 'linewidth', .75)