Tags: cdf, chi-square, complex, comprised, gaussian, identically, independently, matlab, ofi, programming, random, rightimagine, simulation, vector

CDF of chi-square

On Programmer » Matlab

6,695 words with 3 Comments; publish: Wed, 07 May 2008 00:54:00 GMT; (20062.50, « »)

Hi there,

I need some help in getting a simulation right!

Imagine I have a complex Gaussian random vector comprised of

i.i.d(identically & independently distributed) random samples each

having mean 0 and variance \sigma^2.

The probability of the sum of the squares of these Gaussian rv's is a

Chi-square random variable with the degree of freedom equal to the

number of Gaussian rv's .

I am having problems computing the P ( this chi-square random

variable < some constant value). The answer I get is lesser than

the corresponding Monte-Carlo simulation. This is because I think the

function (ncx2cdf which I use to compte the probability with the cdf

of the generalized Chi-square distribution) assumes a normal Gaussian

distribution of the random vector, which is not the case in general.

Say, I assume a 32 - dimensional Gaussian rv, where each sample of

the random vector has a zero mean Gaussian distribution with variance

0.5 and mean 0, what is the probability that the euclidean norm (sum

of the squares) of this rv is less than 15

%%%%%%%%%%%%% MATLAB COD OF MS SIMULATION%%%%%%%%%%%%%%%%%%

% Generating the X sequence

C = 15

X = sqrt(.5) * complex(randn(16,1000000),randn(16,10000

00));

X = [real(X);imag(X)];

X = abs(X).^2;

% Computing the sum of the squares

X = sum(X,1);

XL = length(find(X<C))/length(X)

% Theoretical estimate

XT = ncx2cdf(C,32,0)

Monte-Carlo simulation gives the probability that the sum of the

squares of the samples of the random vector is less than 15 is .4320.

Using ncx2cdf(15,32,0) I get .0046. Could someone help me figure out

as to what I am doing wrong! I have got some normalization wrong! Any

answer along with the correct code would be appreciated!

Thanks for taking thetime to help me out!

Best wishes,

Bikral

All Comments

Leave a comment...

  • 3 Comments
    • Bikral Saluja wrote:

      > %%%%%%%%%%%%% MATLAB COD OF MS SIMULATION%%%%%%%%%%%%%%%%%%

      > % Generating the X sequence

      > C = 15

      > X = sqrt(.5) * complex(randn(16,1000000),randn(16,10000

      00));

      > X = [real(X);imag(X)];

      > X = abs(X).^2;

      > % Computing the sum of the squares

      > X = sum(X,1);

      > XL = length(find(X<C))/length(X)

      Hi Bikral -

      Your problem is the sqrt(.5) -- a chi-square r.v. is defined as sum of squar

      es

      of _standard_ normals, and you're generating scaled normals.

      > % Theoretical estimate

      > XT = ncx2cdf(C,32,0)

      Try ncx2cdf(C/.5,32,0) instead, and I think you'll find the values agree nic

      ely.

      By the way, there's no reason to use the NC chi-square, unless you need a

      non-zero mean for the underlying normals. chi2cdf would be faster.

      Hope this helps.

      - Peter Perkins

      The MathWorks, Inc.

      #1; Wed, 07 May 2008 00:56:00 GMT
    • Peter Perkins wrote:

      >

      > Bikral Saluja wrote:

      >

      > Hi Bikral -

      > Your problem is the sqrt(.5) -- a chi-square r.v. is defined as sum

      > of squares

      > of _standard_ normals, and you're generating scaled normals.

      >

      > Try ncx2cdf(C/.5,32,0) instead, and I think you'll find the values

      > agree nicely.

      > By the way, there's no reason to use the NC chi-square, unless

      > you need a

      > non-zero mean for the underlying normals. chi2cdf would be faster.

      > Hope this helps.

      > - Peter Perkins

      > The MathWorks, Inc.

      >

      Hi Dr. Perkins,

      Thank you for your reply! Well I was trying a similar estimation of

      the probability with the non-central chi-square distribution. Imagine

      I have a random vector as before having a non - zero mean (all are

      non-zero) and each having variance sigma^2. How can I use the

      function ncx2cdf to procure the answer. I am just messing up on the

      normalization. (In the sense that for the Gaussian case, if I want to

      find the cdf of say x, concerning a non-zero mean Gaussian random

      variable, all one does is (x-mu)/sqrt(sigma^2)))

      How do I normalize for the non-central chi-square distribution

      %%%%%%%%%%%%%%%%% MATLAB CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      M = randn(4,1)

      % Generating the X sequence

      X = diag(M)*ones(4,1000000) +

      sqrt(. 5)*complex(randn(4,1000000),randn(4,1000

      000));

      X1 = [real(X);imag(X)];clear X

      X2 = abs(X1).^2;

      % Computing the sum of the squares

      X2 = sum(X2,1);

      XL = length(find(X2<C))/length(X2)

      % Theoretical estimate

      NCP1 = mean(X1,2);

      NCP1 = sum(NCP1.^2,1)/.5;

      XT = ncx2cdf(C,8,NCP1)

      1. A similar example as before

      2. I create M, and then add it to the zero mean Gaussian random

      vector.

      3. Then do as before, computing the Monte-Carlo estimate (XL)

      4. Since the non-centrality parameter of a non-centreal chi-square

      distribution is defined as summation (mu_i/sigma)^2, I compute NCP1

      as above.

      5. I want to find that the sum of the squares is less than some

      constant C.

      I have computed the non-centrality parameter exactly as it states in

      the textbooks, in the sense that it is the sum of the squares of the

      means divided by the variance, which is .5 for each sample. But the

      answers, once again do not tally! I tried computing C as [ C(Actual

      value) - NCP1 / .5 ] but this does'nt give the answer.

      An example : C = 7 , M =[ -0.7944 ; -1.3584; -0.7425; -0.7556]

      gives XL = 0.4812 and XT = 0.0857. My normalization does not give

      the eigfht answer ! Could you kindly tell me as to what I am doing

      wrong!

      Best wishes,

      Bijral

      #2; Wed, 07 May 2008 00:57:00 GMT
    • Bijral wrote:

      > %%%%%%%%%%%%%%%%% MATLAB CODE %%%%%%%%%%%%%%%%%%%%%%%%%%%%%

      > M = randn(4,1)

      > % Generating the X sequence

      > X = diag(M)*ones(4,1000000) +

      > sqrt(. 5)*complex(randn(4,1000000),randn(4,1000

      000));

      > X1 = [real(X);imag(X)];clear X

      > X2 = abs(X1).^2;

      > % Computing the sum of the squares

      > X2 = sum(X2,1);

      > XL = length(find(X2<C))/length(X2)

      > % Theoretical estimate

      > NCP1 = mean(X1,2);

      > NCP1 = sum(NCP1.^2,1)/.5;

      > XT = ncx2cdf(C,8,NCP1)

      > An example : C = 7 , M =[ -0.7944 ; -1.3584; -0.7425; -0.7556]

      > gives XL = 0.4812 and XT = 0.0857. My normalization does not give

      > the eigfht answer ! Could you kindly tell me as to what I am doing

      > wrong!

      Again: your MC simulation is not generating NC chi-square random values, so

      there is no reason to expect thatXL and XT should be the same.

      #3; Wed, 07 May 2008 00:58:00 GMT