Generate random numbers from a mixture distribution
28 views (last 30 days)
Show older comments
Thomas PLOCOSTE
on 2 Aug 2019
Commented: cynthia thing
on 26 Jan 2021
Hi, I would like to generate random numbers from a mixture distribution defined as:
pdf_mix = p*burrpdf(x,alpha,c,k) + (1-p)*wblpdf(x,a,b);
The aim is to compute a kstest2 between the observed data and the mixture distribution.
Is this possible? Can I get some help on this?
T PLOCOSTE
2 Comments
Accepted Answer
John D'Errico
on 3 Aug 2019
Edited: John D'Errico
on 3 Aug 2019
Read David's answer. But let me explain it like this, because I think many people do not understand what a mixture distribution means or what it is.
A mixture distribution means you have two (well, or more) different distributions. I'll call them A and B. With probability p, a sample arises from distribution A, then with probability 1-p, the sample arises from distribution B.
The solution is simple. For each sample, first, choose a uniform random number. If the number is less than or equal to p, then you will sample from distribution A. If the number exceeds p, then you choose from distribution B. Easy, peasy.
Lets see how it might work. I'll pick two distributions that will be clearly different.
Distribution A: Uniform between limits [2,3]
Distribution B: Normally distributed, with mean 0, standard deviation 1.
Lets say 1000000 samples overall, enough to get a nice clean looking histogram. I'll set p as 0.4, arbitrarily. So roughly 40% of the samples will be uniform, 60% normal.
N = 1000000;
p = 0.4;
% U is a binary indicator variable, that tells how to
% sample for each sample. Here, 1 --> A, 0 --> B.
U = rand(N,1) <= p;
nA = sum(U);
nB = N - nA;
% preallocate X
X = zeros(N,1);
X(U) = rand(nA,1)*1 + 2;
X(~U) = randn(nB,1);
% finally, a histogram, just to show we did what was expected.
histogram(X,1000)
![](https://www.mathworks.com/matlabcentral/answers/uploaded_files/232507/image.jpeg)
So we have some overlap between the distributions, but clearly a mixture of a Gaussian and a uniform, as we expected.
sum(U)
ans =
399632
We would have expected roughly 400000 of the events to be from the uniform case. Pretty close.
0 Comments
More Answers (3)
David Goodmanson
on 3 Aug 2019
Edited: David Goodmanson
on 3 Aug 2019
Hi Thomas.
I presume you mean something like the following. A uniform random variable is used to create an index that picks from the first distribution with probability p, and from the second one with probability 1-p. Here I just arbitrarily made two different distributions from the standard normal distribution. You don't say if p = 2253/3849 in your case, but that's not necessary and p could be anything.
% normal distribution parameters
s1 = 1; mu1 = 1;
s2 = 2; mu2 = 5;
p = .2;
n = 1000;
R = zeros(1,n); % will be the final result
r = rand(1,n);
ind = r<p; % index for first distribution
np = sum(ind);
R( ind) = s1*randn(1,np ) + mu1; % np instances
R(~ind) = s2*randn(1,n-np) + mu2; % n-np instances
0 Comments
Jeff Miller
on 3 Aug 2019
> With this formula I obtain a concatenation from 2 distribtutions but not a mixture.
Well, a mixture is just a concatenation from 2 distributions in random order. So you could just randomize the order of R and have your mixture.
As far as kstest2 goes, the order is not important, so you could in fact just use R and get the same kstest2 output.
But if your goal is to test whether the data come from a known distribution, you would really be better to treat that distribution as known and use a one-sample kstest. (I am not sure whether you can build the mixture distribution that you want as a MATLAB distribution object, however, so you might have to code up the kstest yourself if you want to go that route.)
0 Comments
cynthia thing
on 24 Jan 2021
Hi Thomas, thank you for this post.
I am currently finding out about this too.
Hopefully you could help me if you are able to obtain.
Thank you
2 Comments
cynthia thing
on 26 Jan 2021
HI Thomas,
You were so helpful.
However I am wondering if I have a random data, which I would like to try to fit the dato into the function, will there be a generalized formula to it?(for example if I do not have a particular number of data that i know)
Thank you
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!