[Math] Numerical multivariate definite integration

integrationna.numerical-analysis

I need to compute a set of multivariate definite integrals with infinite integration domain

$$\displaystyle \int_{-\infty}^{\infty} \cdots \int_{-\infty}^{\infty} f(x_1,x_2, \ldots , x_n)\;\;dx_1 \cdots dx_n $$

With $n$ typically taking values between 80-150.

I have not found a closed form solution. Besides the function $f$ cannot be decomposed as product of unidimensional functions in $x_i$.

Does there exist a matlab toolbox or any software package with a matlab interface to do that?

Note: I am totally newbie in numerical integration. This problem comes from a totally different field.
Thanks.

Edit: after some minor manipulation the integral look like this

Let $\boldsymbol{x},\boldsymbol{b}\in\mathbb{R}^n$ and $\boldsymbol{A}\in\mathbb{R}^{n\times n}$

$$ I = \int_{-\infty}^{\infty} exp ( -\frac{1}{2} x^T A x +
b^T x ) \; \prod_{i=1}^n
exp(-|x_i|) \; dx = \int_{-\infty}^{\infty}
exp ( -\frac{1}{2} x^T A x +
b^T x –
\sum_{i=1}^n |x_i| )
\; dx
$$

$$I=\int_{-\infty}^{\infty} \cdots\int_{-\infty}^{\infty} exp(-\frac{1}{2}\sum_{i=1}^n \sum_{j=1}^n a_{i,j}\;x_ix_j + \sum_{i=1}^n b_ix_i-\sum_{i=1}^n |x_i|)\; dx_1 \cdots dx_n$$

Whereas I try to find an analytical solution I have been playing with numerical integration. I wrote this script to integrate a multivariate normal distribution(because is an easy function and we know the result is 1). The problem is that the computation time grows with the dimension and the script is at moment too slow for my requirements. If anyone can give me advise to speed up the code I will be grateful.

clc

close all

clear all

pmax=30;

t=zeros(pmax,1);
result=zeros(pmax,1);
sigma=2;

for iter=1:pmax

    str = [' Iteration  ', num2str(iter), ' of ',num2str(pmax)];
    disp(str)

    p=iter;
    S=(sigma^2)*eye(p);
    invS=inv(S);
    detS=det(S);

    %x=[x1;...; xp]
    x=sym('x',[p,1]);

    % x'=  complex conjugate transpose
    % x.'= transpose without conjugation
    xT=x.';
    f=exp(-0.5*xT*invS*x);
    fk=f;
    tstart=tic;
    for k=1:p
        fk=int(fk,x(k),-inf,inf);
    end
    t(iter)=toc(tstart);
    aux=1/sqrt((2*pi)^p);
    aux=aux/sqrt(detS);
    aux=aux*fk;
    result(iter)=aux;
end

figure

plot(result)

title('Result (Ideally always 1)')

figure

plot(t)

title('time elapsed')

Best Answer

Look up "gaussian integral" on Wikipedia. I think you can handle the absolute value portion by splitting the range, and you will get the integral in terms of $ \det A$. Take a look.

Tom

Related Question