Integrating to get volume under bivariate normal distribution

Question

supernoob 2018-7-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/409280-integrating-to-get-volume-under-bivariate-normal-distribution

编辑： supernoob 2018-7-8

data_from_mc.mat

Hi there, I have 2 gaussian random variables which together form a bivariate normal distribution. Let's call them A and B. Say I did this:

AB = [A,B]; %A and B are gaussian distributed RV's
AB_mean = mean(AB); 
AB_std = std(AB); 
AB_var = var(AB); 
AB_covar = cov(AB);

I can create and plot the multivariate normal distribution as follows:

xsteps = linspace(-5e-3,5e-3, 1000);
ysteps = xsteps;
[X,Y] = meshgrid(xsteps,ysteps);
F = mvnpdf([X(:) Y(:)],AB_mean,AB_covar);
F = reshape(F,length(ysteps),length(xsteps));
figure(53); s = surf(xsteps,ysteps,Fnorm, 'FaceAlpha',0.5); s.EdgeColor = 'none';
xlabel('x'); ylabel('y'); zlabel('Probability Density');

How do I integrate this distribution in 2d to get the volume under some portion of the surface? I haven't been able to find a way which works. Note what I really want is to integrate over the pdf in polar coordinates. I tried doing this (using the cartesian definition of the pdf from http://mathworld.wolfram.com/BivariateNormalDistribution.html):

AB_mean = mean(AB); mu1 = AB_mean(1); mu2 = AB_mean(2);
AB_std = std(AB); sigma1 = AB_std(1); sigma2 = AB_std(2);
AB_var = var(AB); var1 = AB_var(1); var2 = AB_var(2);
AB_covar = cov(AB);
rmax = 2e-3;
fun = @(x,y) (1/(2*pi*sigma1*sigma2*sqrt(1-rho^2)))...
    .*exp((-1/(2*(1-rho^2)))*(((x-mu1)/var1).^2+((y-mu2)/var2).^2-((2*rho*(x-mu1).*(y-mu2))/(sigma1*sigma2))));
polarfun = @(theta,r) fun(r.*cos(theta),r.*sin(theta)).*r;
q = integral2(polarfun,0,2*pi,0,rmax);

But this always gives me zero, even if I make r very large, and I can't figure out what is wrong. If things were working, I would expect to get 1 for very large r. Please help! I have attached my 2 RV's as a matrix AB.

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

supernoob 2018-7-8

编辑：supernoob 2018-7-8

Update! I found my mistake: I was squaring the variance in the denominators of the exponential but they are already squared. Now everything works as expected. I'll post this as the answer and leave this question up in case it is useful to somebody. If you think I should take it down, let me know.

请先登录，再进行评论。

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

Answer 1

supernoob 2018-7-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/409280-integrating-to-get-volume-under-bivariate-normal-distribution#answer_327983

在 MATLAB Online 中打开

The mistake is squaring the variance in the denominator of the exponential terms. The variance is already squared. The correct expression for the probability density function is the following:

fun = @(x,y) (1/(2*pi*sigma1*sigma2*sqrt(1-rho^2)))...
    .*exp((-1/(2*(1-rho^2)))*(((x-mu1).^2/var1)+((y-mu2).^2/var2)-((2*rho*(x-mu1).*(y-mu2))/(sigma1*sigma2))));