Integrating multivariate Gaussian pdf in 4 dimensions

I have been trying to integrate a multivariate Gaussian pdf in 4 dimensions, but got the error message.

Any help would be appreciated.

Thanks!

f = @(x,y,z,w) mvnpdf([x y z w],zeros(1,4),eye(4));
val = integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0),x,y),-inf,0,-inf,0)

Error using mvnpdf (line 67)
X and MU must have the same number of columns.

Error in solution>@(x,y,z,w)mvnpdf([x,y,z,w],zeros(1,4),eye(4)) (line 1)
f = @(x,y,z,w) mvnpdf([x y z w],zeros(1,4),eye(4));

Error in solution>@(z,w)f(x,y,z,w) (line 2)
val = integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0),x,y),-inf,0,-inf,0)

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, ...

Error in integralCalc (line 99)
[q,errbnd] = vadapt(vfunB,interval, ...

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, ...

Error in integralCalc (line 99)
[q,errbnd] = vadapt(vfunB,interval, ...

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Error in solution>@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0) (line 2)
val = integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0),x,y),-inf,0,-inf,0)

Error in solution>@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0),x,y) (line 2)
val = integral2(@(x,y)arrayfun(@(x,y)integral2(@(z,w)f(x,y,z,w),-inf,0,-inf,0),x,y),-inf,0,-inf,0)

Error in integral2Calc>@(y)fun(xi*ones(size(y)),y) (line 18)
@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions), ...

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, ...

Error in integralCalc (line 99)
[q,errbnd] = vadapt(vfunB,interval, ...

Error in integral2Calc>@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integral2Calc>@(x)arrayfun(@(xi,y1i,y2i)integralCalc(@(y)fun(xi*ones(size(y)),y),y1i,y2i,opstruct.integralOptions),x,ymin(x),ymax(x)) (line 17)
innerintegral = @(x)arrayfun(@(xi,y1i,y2i)integralCalc( ...

Error in integralCalc>iterateScalarValued (line 334)
fx = FUN(t);

Error in integralCalc>vadapt (line 148)
[q,errbnd] = iterateScalarValued(u,tinterval,pathlen, ...

Error in integralCalc (line 99)
[q,errbnd] = vadapt(vfunB,interval, ...

Error in integral2Calc>integral2i (line 20)
[q,errbnd] = integralCalc(innerintegral,xmin,xmax,opstruct.integralOptions);

Error in integral2Calc (line 7)
[q,errbnd] = integral2i(fun,xmin,xmax,ymin,ymax,optionstruct);

Error in integral2 (line 105)
Q = integral2Calc(fun,xmin,xmax,yminfun,ymaxfun,opstruct);

Answer 1

[Torsten]

Use "mvncdf":

https://uk.mathworks.com/help/stats/mvncdf.html

If your function is more complicated than posted, try "integralN":

https://uk.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

But it will usually be less precise than "mvncdf".

f = @(x,y,z,w) mvnpdf([x y z w],zeros(1,4),eye(4));
val1 = mvncdf(zeros(1,4),zeros(1,4),eye(4))

val1 = 0.0625

val2 = integralN(@(x,y,z,w)arrayfun(f,x,y,z,w),-Inf,0,-Inf,0,-Inf,0,-Inf,0)

[LUU KHOI]

Thanks!