Noise samples of gaussian mixture distribution

4 次查看（过去 30 天）
Harsha Vats
Harsha Vats 2016-6-28
编辑： John D'Errico 2016-6-30
I'm new to matlab. I want to generate noise samples of Gaussian mixture with PDF= sqrt((u.^3)./pi).*exp(-u.*(x.^2)) + sqrt((1-u)^3/pi).*exp(-(1-u).*x.^2);
The length of noise sample is (1,200) Please help me out.
  2 个评论
José-Luis
José-Luis 2016-6-28
编辑：José-Luis 2016-6-28
The function you describe might yield a bunch if imaginary numbers. How do you intend to deal with that?
Harsha Vats
Harsha Vats 2016-6-28
'u' lies between 0 to 0.5. for that range of u I don't think the function will yield any imaginary values.

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

John D'Errico
John D'Errico 2016-6-28
编辑：John D'Errico 2016-6-28
Ok, so you have formulated what seems to me to be a rather strange mixture model.
You have two modes, with variances that are related to each other, the mean of both terms is zero.
So, if the PDF is:
PDF = sqrt((u.^3)./pi).*exp(-u.*(x.^2)) + sqrt((1-u)^3/pi).*exp(-(1-u).*x.^2)
So we have two terms there.
syms u positive
syms x real
Look at the first term:
pdf1 = sqrt((u.^3)./pi).*exp(-u.*(x.^2));
int(pdf1,x,[-inf,inf])
ans =
u
So, it integrates to u. The second term is similar,
pdf2 = sqrt(((1-u).^3)./pi).*exp(-(1-u).*(x.^2));
int(pdf2,x,[-inf,inf])
ans =
piecewise([u == 1, 0], [1 < u, Inf*1i], [u < 1, 1 - u])
For u <= 1, it integrates to 1-u.
So the integral of the sum is indeed 1, therefore it is a valid PDF.
Essentially, you have two Gaussian modes, such that the variances are 2/u and 2/(1-u) respectively. You have chosen the mixture parameter as the inverse of those variances. Effectively, when u is small (or very near 1), you have a mixture distribution that rarely, you will generate large outliers. When u is exactly 1/2, this reduces to a standard normal distribution, with a unit variance.
I'm not sure why, but I suppose that is not really important. Normally, one would choose a mixture parameter that is independent of the variances.
The solution is simple.
1. Pick some value for u. u determines both variances for each Gaussian mode, as well as the mixture coefficient.
2. For a sample size of N, pick N uniformly distributed random numbers. These will determine which mode you will sample from. Then just use a classical tool like randn to do the sampling.
For a sample size of 1e7, and some arbitrary value for u, say 1/100. I've picked a tiny value for u to make the result clear. I've also chosen a large sample to make it easier to see as a histogram.
N = 10000000;
u = 0.01;
% random mixture selection
% Here s will be 1 if the element is sampled from first pdf.
% s will be 0 if the element is sampled from second pdf.
s = rand(1,N) < u;
% sample using randn, with a unit variance initially
x = randn(1,N);;
% now scale each element based on the desired variance
x(s) = x(s)*sqrt(2/u);
x(~s) = x(~s)*sqrt(2/(1-u));
hist(x,1000)
So, for this fairly tiny value of u, see that the histogram looks much like a traditional Gaussian, BUT with very wide tails.
I'll next plot the separate pieces of the PDF, split into two parts. Here I've chosen u as 0.25.
pdf1 = @(x,u) sqrt((u.^3)./pi).*exp(-u.*(x.^2));
pdf2 = @(x,u) sqrt(((1-u).^3)./pi).*exp(-(1-u));
h1 = ezplot(@(x) pdf1(x,0.25),[-10,10]);
hold on
h2 = ezplot(@(x) pdf2(x,0.25),[-10,10]);
title 'pdf1 and pdf2, split apart'
So, the sampling was rather simple. No need to do anything strenuous. I still have no idea why you want this distribution, but whatever floats your boat is a good thing.
  4 个评论
显示 2更早的评论
Harsha Vats
Harsha Vats 2016-6-29
Sir,
When I tried to convert this equation in form of sum of gaussian, I was getting u*Gaussian with variance sqrt(1/(2u)) + (1-u)*Gaussian with variance sqrt(1/(2(1-u))).
How in your case the variances are sqrt(2/u) and sqrt(2/(1-u))?
John D'Errico
John D'Errico 2016-6-30
编辑：John D'Errico 2016-6-30
Sheepish response. :)
Oh, because I did the algebra without even using pencil and paper, so I made a mistake? In fact as I was thinking about it long after I answered your question, I was thinking I had possibly put that factor of 2 in the wrong place, but then I forgot to check it.
What is a factor of 2 either way anyway? Did my mission to Mars end up in the wrong place, so a crash landing? :)

请先登录，再进行评论。

更多回答（2 个）

José-Luis
José-Luis 2016-6-28
编辑：José-Luis 2016-6-28
Proof of concept, not the smart way to do it
If this was to be done in an efficient way:
  1. Derive the cdf.
  2. Get a random number between 0 and 1: rand()
  3. Compute the inverse cdf for that number.
This is possible to do analytically.
We can do it numerically for u = 0.5:
u = 0.5;
myfun = @(x) sqrt((u.^3)./pi).*exp(-u.*(x.^2)) + sqrt((1-u)^3/pi).*exp(-(1-u).*x.^2);
Getting the random value
val = rand;
The integral() part computes the cdf corresponding to a value in the interval ]0,1[
intFun = @(x) integral(myfun,-Inf,x) - val;
The value you are looking for will be the root for the above function:
your_value = fzero(intFun,0);
And you can check that you actually obtained the right result.
integral(myfun,-Inf,your_value)
Which should be equal to val
Inefficient: loop 200 times to get your random vector. Better: solve it analytically, or see if your function is a built-in distribution in Matlab (I don't think so though).
  2 个评论
Ankita kumari
Ankita kumari 2016-6-28
I'm unable to understand why you subtracted val from integral value. The integral will give the cdf. Earlier you said that after deriving the CDF we need to compute the inverse of that.
I have computed the cdf it is:
cdf=u*(1-qfunc(x/(sqrt(1/(2*u))))) + (1-u)*(1-qfunc(x/(sqrt(1/(2*(1-u))))));
next what I need to do?
José-Luis
José-Luis 2016-6-28
Because you need to find the value for which the cdf gives val and finding the roots of the cdf minus said value is one way to go.
I you have the cdf then you can ignore the integral() part of what I gave you and use that instead. Better yet, compute the inverse cdf and obtain your value directly.

请先登录，再进行评论。

Bjorn Gustavsson
Bjorn Gustavsson 2016-6-28
For any given u between 0 and 1/2 this looks like the sum of 2 zero-centred normal-distributions (however there seems to me that the scaling is a bit off, so maybe check that the pdf actually adds up to 1?), in that case it should be as simple as
x_samples = sigma1*randn(sz_samples) + sigma2*randn(sz_samples)
right?
HTH
  5 个评论
显示 3更早的评论
John D'Errico
John D'Errico 2016-6-28
编辑：John D'Errico 2016-6-28
Actually, this procedure is not correct. It is not the same to simply add samples from two different gaussians as what the OP asked for.
The sum of x and y, where x and y are both gaussian with different variances has a gaussian distribution also. The PDF posed by the OP was a MIXTURE distribution. That is a VERY different animal.

请先登录，再进行评论。

类别

在 Help Center 和 File Exchange 中查找有关 Gaussian Mixture Distribution 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by