# generate random number from a function that serves as PDF ?

58 views (last 30 days)
Hakim Jemaa on 8 Apr 2020
Answered: Torsten on 8 Apr 2020
Hello! I have mean zero second order gaussian distribution model and I want to generate random variables from this distribution within a specific interval [ -pi/2, pi/2 ]. the PDF depends on 4 parameters that are 2 variances and 2 other coefficients for the seocnd order. I built the PDF within a mtlab function and I would like to know how can I generate random numbers from this PDF ??
Thanks
##### 2 CommentsShowHide 1 older comment
Hakim Jemaa on 8 Apr 2020
Of course
y=(w1/(sqrt(2*pi)*sigma1))*exp(-(x^2/(2*sigma1^2)))+(w2/(sqrt(2*pi)*sigma2))*exp(-(x^2/(2*sigma2^2)));
sorry for the bad equation writing. my parameters are w1, w2, sigma1 and sigma2 which are the variances

Ameer Hamza on 8 Apr 2020
You can use the inverse transform sampling method: https://en.wikipedia.org/wiki/Inverse_transform_sampling to generate a random number from your specified distribution. See my answer on this question for details: https://www.mathworks.com/matlabcentral/answers/514287-how-to-generate-random-variable-x-for-specific-cdf-fx-1-exp-x-2-2-sigma-2-for-sigma-2
##### 2 CommentsShowHide 1 older comment
Ameer Hamza on 8 Apr 2020
Hakim, the pdf function you shared cannot be inverted in closed-form, so the method I mentioned cannot be directly applied to it. There are ways to apply this method, but it will be a bit more complicated. The easiest way is to use this package from FEX: https://www.mathworks.com/matlabcentral/fileexchange/26003-random-numbers-from-a-user-defined-distribution. This package will also restrict the random number in the specified range. Download this package and place it in MATLAB's path. Then run the following code. The variable rand_num contains 100000 random samples from the pdf you gave.
w1 = 0.5;
w2 = 0.5;
sigma1 = 1;
sigma2 = 0.1;
pdf = @(x) (w1/(sqrt(2*pi)*sigma1))*exp(-(x.^2/(2*sigma1^2)))+(w2/(sqrt(2*pi)*sigma2))*exp(-(x.^2/(2*sigma2^2)));
sample = linspace(-pi/2, pi/2, 1000);
pdf_sample = pdf(sample);
rand_num = randpdf(pdf_sample, sample, [100000 1]);
histogram(rand_num,200)

Torsten on 8 Apr 2020
N = 1000;
U = rand(N,1);
n1 = numel(U(U<=w1));
n2 = N - n1;
X1 = sigma1*randn(n1,1);
X2 = sigma2*randn(n2,1);
X = [X1,X2];
X = X(abs(X)<=pi/2);