generate random number from a function that serves as PDF ?

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

[Ameer Hamza]

Can you share your pdf function?

[Hakim Jemaa]

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

Answer 1

[Ameer Hamza]

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

[Hakim Jemaa]

Thanks. If I understood correctly, I have to get the CDF out of my PDF function, inverse it, generate a random number and compute its value ??

how to get the cdf ? and how to shrink the results to a specific segment ?

Thank you for your time

[Ameer Hamza]

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)

Answer 2

[Torsten]

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