Marginal Density from a joint DIstribution

Hey, I have a really simple question. How can I obtain a marginal density fx(x) from a joint distribution (x,y) ? In my case the joint distribution follows a log-normal distribution. I cannot use Quad since it requires both integrals (x and y). Thanks a lot for your help. Mo

Answer 1

[Mike Hosea]

You might try to do it symbolically with INT. Numerically, you could do this:

fx = @(t)arrayfun(@(x)integral(@(y)f(x,y),-inf,inf),t)

Naturally you would use whatever the correct range is on y if it's not -inf to inf. I just wrapped it with arrayfun so you could easily integrate it or plot it. I also used arrayfun because you can substitute quadgk for integral if you don't have R2012a. With the new integral function in particular you also have the option of using 'ArrayValued' option. Here's an example

BivariateNormalPDF = @(x,y,mux,sigmax,muy,sigmay,rho) ...
  exp(-(((x-mux)/sigmax).^2 ...
       + ((y-muy)/sigmay).^2 ...
       - 2*rho*((x-mux).*(y-muy)/(sigmax*sigmay)) ...
       )/(2*(1-rho*rho)) ...
     )/(2*pi*sigmax*sigmay*sqrt(1-rho*rho));
f = @(x,y)BivariateNormalPDF(x,y,3,1,1,2,0.5);
fx = @(x)integral(@(y)f(x,y),-inf,inf,'ArrayValued',true);

You can now plot or integrate fx.

>> x = -2:0.1:8;
>> plot(x,fx(x));
>> integral(fx,-inf,inf)
ans =
            1

-- Mike

[Mo]

Thanks a lot Mike ! It's exactly what I was looking for. Still, I have two problems.

1. When I re-do your example using the built-in function 'MVNPDF', I don't get the same figure as you :

x = 0:0.1:8;
  y = 0:0.1:8;
BivariateNormalPDF = @(x,y,mux,sigmax,muy,sigmay,rho) ...
  exp(-(((x-mux)/sigmax).^2 ...
       + ((y-muy)/sigmay).^2 ...
       - 2*rho*((x-mux).*(y-muy)/(sigmax*sigmay)) ...
       )/(2*(1-rho*rho)) ...
     )/(2*pi*sigmax*sigmay*sqrt(1-rho*rho));
f = @(x,y)BivariateNormalPDF(x,y,3,1,1,2,0.5);
  plot(x,f(x,y),'-r');  
  hold on;    

    % Using MVNPDF
   mu = [3 1];
  Sigma = [1 .5; .5 2];
  phi =  @(x,y) mvnpdf([x(:) y(:)],mu,Sigma) ;
  plot(x,phi(x,y),'-k')
  hold on;

2. I use the function 'MVNPDF' since I need to get the log normal distribution :

    logphi = @(x,y)  reshape(mvnpdf([log(x(:)) log(y(:))],mu,Sigma) ,size(x))./x./y

The problem is, when I use 'MVNPDF', I cannot compute/plot the marginal distributions : fx = @(x)integral(@(y)f(x,y),-inf,inf,'ArrayValued',true); phix = @(x)integral(@(y)phi(x,y),-inf,inf,'ArrayValued',true); logphix = @(x)integral(@(y)logphi(x,y),-inf,inf,'ArrayValued',true); plot(x,fx(x),'-g'); plot(x,phix(x),'-o'); plot(x,logphix(x),'-o');

Thanks a lot, Mo

[Mike Hosea]

The equivalent Sigma matrix is

v12 = sigmax*sigmay*rho; Sigma = [sigmax.^2 v12; v12 sigmay.^2];

I'll figure out the remaining vectorization issue and get back to you.

[Mike Hosea]

You have to keep in mind that INTEGRAL is passing in a scalar for y. The integrand functions will not be convenient for plotting because they have to be coded for an array input for x and a scalar input for y.

phi = @(x,y) reshape(mvnpdf([x(:),y*ones(numel(x),1)],mu,Sigma),size(x)) ; logphi = @(x,y)reshape(mvnpdf([log(x(:)),log(y)*ones(numel(x),1)],mu,Sigma),size(x))./x./y;

Also, you will probably need to re-think the range on the integral computing logphix. Maybe you want to integrate from 0 to inf? There is still a problem because the integrand returns NaN when either x or y is zero (it's doing 0/0), so you might want to plug that hole somehow. You can, of course, define logphi as a .m function and avoid that with an extra line of code. A quick and dirty way out is to integrate from realmin (and also define x and y at top to start at realmin rather than zero).

[Mo]

Everything works now ! Thanks a lot for your help Mike.

[Antonio Marino]

Hi Mike, I would like to ask you if the same procedure is suitable to calculate the conditional distribution dividing the joint and marginal.

If it is possible Can you explain how?

Would be very useful to me for my thesis.

thank you in advance.