Double Integral of Bivariate joint normal distribution

I want to apply "integral2/dblquad" on bivariate normal joint pdf. The formula for multivariate joint pdf is:

mu = [1 -1]; SIGMA = [.9 .4; .4 .3]; X = mvnrnd(mu,SIGMA,10); p = mvnpdf(X,mu,SIGMA);

This formula is working, however, i don't want to generate random numbers by using "X = mvnrnd(mu,SIGMA,10);" instead I want to double integrate the function " mvnpdf(X,mu,SIGMA)" by keeping "X" as variable consisting of two variables x & y. Can you please help me in finding double integral by using "X" as vector consisting of random variables x & y. The double integral that I want to calculate is:

integral2(@(x,y)mvnpdf(X,mu,SIGMA), 0,400,0,500) in which x & y are variables.

Thanks in advance.

Answer 1

[Roger Stafford]

Using 'mvnrnd' to do double integration is an extremely poor method of integration. Just use 'integral2' on the integrand 'mvnpdf(X,mu,SIGMA)' as it was intended to be used to do your double integration.

[Muhammad]

Dear Roger Stafford

Thanks.

I tried but I need to define "X" as else it can't be solved. Can you please help me in defining "X"?

[Roger Stafford]

Just create a subfunction:

   mu = [1 -1]; sigma = [.9 .4; .4 .3];
   function P = mypdf(x,y)
     P = mvnpdf([x(:),y(:)],mu,sigma);
   return

and call it with

   integral2(@mypdf,...)

where you set the appropriate limits for the region you are integrating over. See

   http://www.mathworks.com/help/matlab/ref/integral2.html

[Muhammad]

Dear Roger,

I tried it already but its not working unfortunately. Is it working for you, if yes then please send me the detailed code.

[Roger Stafford]

In what way does what you tried not work? Please show the code you used.

[Muhammad]

First of all, I developed this subfunction:

function z18=depd(x,y)

global mu SIGMA

z18=mvnpdf([x(:),y(:)],mu,SIGMA);

end

and then called this in integral2 in command window as follow:

mu = [700 600]; SIGMA = [50 .4; .4 40];

integral2(@(x,y)depd(x,y),0,400,50,300)

the error it shows is: "error(message('MATLAB:integral2:funSizeMismatch'))"

Can you please check if any mistake. waiting for your prompt response.

[Mike Hosea]

As the error message says, depd(x,y) does not return an array the same size as the input arrays x and y. That's because you used x(:) and y(:) to conform to the requirements of mvnpdf. So just reshape the output:

z18 = reshape(mvnpdf([x(:),y(:)],mu,SIGMA),size(x));

and, while you don't have to change this, since depd is already a function of only x and y, you can refer to your m-file function like so

integral2(@depd,0,400,50,300)

Alternatively, using global variables seems heavy-handed here. You could snapshot the current values of mu and SIGMA when you create your function handle:

integral2(@(x,y)reshape(mvnpdf([x(:),y(:)],mu,SIGMA),size(x)),0,400,50,300)

[Muhammad]

Thanks Roger Stafford and Mike Hosea. It is really helpful for me. Though its working but I think, I should cross check. What do you suggest. Thanks again. Regards