﻿Matlab: CDF of chi-square - matlab.itags.org

CDF of chi-square

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