I want a numerical or fast multidimensional integral

I have the following code:

 if true
clc
clear all
close all
tic
syms x1 x2 x3 x4 x5
warning off;
X=[x1-12.568; x2+9.568; x3-11.486; x4-4; x5-0.8];%this is [x-E[x]] matrix.
Cx=[39.4573    0.2350    1.5190   25.6028    1.6754; 0.2350    0.0089    0.0285    0.5958    0.0344;...
    1.5190    0.0285    0.2337    1.8325    0.1193; 25.6028    0.5958    1.8325  667.1131    5.5158;...
    1.6754    0.0344    0.1193    5.5158    1.3541];%this the covariance matrix.
DeterminantCx=det(Cx)
CxInverse=inv(Cx);%this is the inverse matrix of the covariance matrix (Cx).
DeterminantCxInverse=det(CxInverse)
Xtranspose=transpose(X);
Exponent=Xtranspose*CxInverse*X;
fJointDensity=(1/(sqrt(DeterminantCx)*(2*pi)^(5/2)))*exp(-(Exponent)/2);
I1=int(fJointDensity,x1,-inf,inf);
I2=int(I1,x2,-inf,inf)
I3=int(I2,x3,-inf,inf)
I4=int(I3,x4,-inf,inf)
I5=int(I4,x5,-inf,inf)
t=['Time Elapsed=',num2str(toc), ' sec'];
disp(t)
end

The expected output for I5 is 1, but the problem is that performing the integrations using "int" takes a lot of time. This program is an example of my main problem, but if I can run this program very fast then I can solve my problem.

So I am searching for a numerical integration or a fast integration method.

Thanks a lot.

Hadi.

[John D'Errico]

You and everyone else recently. Magic does not happen, except for Harry Potter.

Of course, in this problem, the answer is trivial, 1. And while you MIGht be able to reduce the complexity of this problem somewhat, it is impossible without knowing what your limits actually are, or what is the real problem.

Without that information, the simplest answer is to buy a faster computer.

[hadi]

John thank you for you reply, but I can tell you that that these limits are possible to be the limits of my actual program. My actual program is more than 600 lines so how come you want me to put it here?!....As I see this example explains my problem 100%. And btw I did not get an answer for this code even on the university supercomputer, so its not a computer problem. The available matlab numerical integral functions are up to three integrals, so my question was if there is a numerical function that can handle more than three integrals.

Answer 1

[John D'Errico]

(I've got a bit more time now, and since Mike has not chimed in, I'll offer an answer.)

The point is, there are things you MIGHT do, but how can we possibly help you? What do you know or not know about numerical integration? What characteristics does your function have? Are the limits really infinite, or are they finite? The problem is, one could write a book about integration. Oh, thats right, there are already lots of books on the subject. But without knowing something about what you have to solve, the best answer is to buy a faster computer.

There are no n-fold integrations in MATLAB, but I recall that Mike Hosea has posted a higher degree integration tool on the FEX, here:

http://www.mathworks.com/matlabcentral/fileexchange/47919-integraln-m

It handles 4,5,6 fold integrations, I think even allowing infinite limits.

That is the good news. The bad news is to expect it to be seriously slow. Iterated integration tools like this tend to take a lot of function evaluations. The reason is simple. To understand, lets make a simple test...

function [f,nx] = testfun(x)
persistent nevals
if isempty(nevals)
  nevals = 0;
end
nevals = nevals + numel(x);

f = 1/sqrt(2*pi)*exp(-x.^2/2);

if nargout > 1
  nx = nevals;
end

That little function counts the number of function evaluations integral makes for a simple Gaussian.

integral(@testfun,-inf,inf)
ans =
            1
[f,nx] = testfun(0)
f =
      0.39894
nx =
   211

So to compute the integral of ONE gaussian in one dimension, it took 211 function evaluations.

But you have a 5-fold integral. Each of these integrations is of similar complexity, so you might expect to need

211^5
ans =
   4.1823e+11

so about 4e11 function evaluations. That is 400 billion evals. If your function has 600 lines of code, how long will it take to execute? If it is a bit slow, and you get only 1 eval per second, then it would take

211^5/30e6
ans =
        13491

YEARS to finish. Even if one call to your 600 line function takes 0.001 seconds to execute, that is still over 13 YEARS to terminate.

Even if some of those integrations are simpler, and take only half as many calls, go get a cup of coffee and a good book, as you will be waiting for a bit. "War and Peace" might be a good choice of books for this one, or better yet, the complete works of Shakespeare, and expect to re-read whatever book you have quite a few times over. The good news is when it gets done, you will be a true Shakespeare scholar.

My point about not telling us anything is still valid. Are there symmetries you can take advantage of? Or, if your function has characteristics like that of a Gaussian as you show, do you know anything about Gauss-Hermite quadrature? This could easily be a boon for you, or a bust. I cannot know, since you have told us nothing.

[hadi]

John I really want to thank you for all your efforts with me, and I can see your point very clear.

I will explain more why I am doing this, actually I have variable limits that can take any value from -inf to inf, except one variable that I want the final answer to be in terms if it.

Ex:

Assume I have 3 random variables: X1, X2, X3 which have the Gaussian joint density f(x1,x2,x3), now what I am concerned about is to have a mathematical close (or approximated) expression for P(X1<x1, X2<4, X3<inf) while this is an example of what I might have, so as we can see that the limits for X2 and X3 are known for a specific case but not X1.

So that my first question was how we can do these integrations (i.e I1, I2, I3,...) faster so that I can have at the end a mathematical expression in terms of X1, thats why I used this clumsy method, but now I can see from your answer that this is not the way to solve the problem.

But now I can change the question and say how does matlab find the value using "mvncdf"? because if I know the procesure I may be able to find a mathematical expression in terms of X1 because "mvncdf" gives only numbers and does not give or accept "syms" expressions.

And maybe another possible solution could be the Riemann's sum.