Marginal Density from a joint DIstribution

Mo
2012 8 7
1 个回答
回答已采纳
更新时间：2020 11 19
11 次查看（30 天）

0 个投票

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

请先登录，再回答此问题。

 采纳的回答

Mike Hosea
Mike Hosea 2012-8-7

0 个投票

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

5 个评论
显示 3更早的评论 隐藏 3更早的评论

Mo
Mo 2012-8-8
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
Mike Hosea 2012-8-8
编辑：Mike Hosea 2012-8-8
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
Mike Hosea 2012-8-9
编辑：Mike Hosea 2012-8-9
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
Mo 2012-8-9
Everything works now ! Thanks a lot for your help Mike.
Antonio Marino
Antonio Marino 2020-11-19
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.

请先登录，再进行评论。

更多回答（0 个）

类别

在 帮助中心 和 File Exchange 中查找有关 Numerical Integration and Differentiation 的更多信息

标签

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by