Noise samples of gaussian mixture distribution

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.

[José-Luis]

The function you describe might yield a bunch if imaginary numbers. How do you intend to deal with that?

[Harsha Vats]

'u' lies between 0 to 0.5. for that range of u I don't think the function will yield any imaginary values.

Answer 1

[John D'Errico]

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.

[Harsha Vats]

Thank you so much. This was really helpful. This distribution was given in a paper "Design and performance analysis of a signal detector based on suprathreshold stochastic resonance" by V.N Hari. This was a reference paper.

[John D'Errico]

It is an interesting variation of a Gaussian that I've never seen before this.

[Harsha Vats]

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]

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? :)

Answer 2

[José-Luis]

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

[Ankita kumari]

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]

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.

Answer 3

[Bjorn Gustavsson]

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

[John D'Errico]

Yes. You are correct that it is not properly scaled for a pair of Gaussian modes as a mixture. As such, the PDF is incorrect.

[Ankita kumari]

I have checked the pdf, the integration from minus infinity to infinty gives 1.

The method you suggested

x_samples = sigma1*randn(sz_samples) + sigma2*randn(sz_samples)

what does sz_samples contains?

[Bjorn Gustavsson]

OK, in that case sz_samples should be your noise-sample-size i.e. [1,200].

[John D'Errico]

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.

[Ankita kumari]

Yes, you're right. This gives wrong results.