How to do a double integration of a multivariate normal probability density function

Question

Trond Oesten 2015-3-12

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/182806-how-to-do-a-double-integration-of-a-multivariate-normal-probability-density-function

回答： Mike Hosea 2015-3-13

Hi,

I'm trying to do a double integral on a nultivariate normal probability density function. The two variables I want to integrate over is U and K. I created a function (fun) with the formula for multivariate normal pdf and used this to performe the double integral. When I do it like this I only get errors. Is there an easier approach?

My code:

   clc; clear all; close all;
   C = [23.875 15.75281; 15.75281 93.9842]; % C = [sigma_u^2 ro*sigma_u*sigma_k;... 
                                                   ro*sigma_u*sigma_k sigma_k] 
   mu_u = 1788.2058;
   mu_k = 70.8489;
   mu = [mu_u;mu_k];
   fun = @(U,K) ((1/sqrt(det(C)*(2*pi)^2))*exp(-0.5*transpose([U;K]-mu)*inv(C)*([U;K]-mu)));
   q = integral2(fun,1700,1900,30,120);

Thanks!

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

Greig 2015-3-13

The integral functions assume that the function being integrated is vectorized and speed the integration process by passing vectors of U and K into your function. The way your function is written, however, cannot handle vectors of U and K. Later, either I, or someone else will help you vectorize it (I have no time right now).

请先登录，再进行评论。

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

Answer 1

Mike Hosea 2015-3-13

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/182806-how-to-do-a-double-integration-of-a-multivariate-normal-probability-density-function#answer_171172

在 MATLAB Online 中打开

If you have the statistic toolbox, you'll want to use MVNCDF for a multivariate normal distribution. However, in case not everything you want to do is multivariate normal, let's fix your current approach. INTEGRAL, INTEGRAL2, and INTEGRAL3 assume that they can send in many U and K values at once and that you will return an array of corresponding p values. Yours doesn't work with non-scalar inputs. Assuming p = fun(U,K) works properly for all scalar U and K in the desired range of inputs, you can make it work with arrays using ARRAYFUN.

q = integral2(@(U,K)arrayfun(fun,U,K),1700,1900,30,120);

ARRAYFUN just takes the array inputs and loops over the entries to construct an array output.