Numerical multivariate definite integration

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')
• The best method probably depends on the form of f, which is a question that may be too difficult for computer algebra packages. I suggest reading up on Monte-Carlo integration. A quick search brought up venus.unive.it/r.casarin/Summer2012/SlidesCasarinMC.pdf
– user25199
Jun 18 '13 at 9:02
• Up to scaling you want to find the expectation of $f(x)=\prod\exp(-|x_i|)$ when $x$ has a normal distribution. Look up how to sample random points from that normal distribution, then the average value of $f(x)$ for such points will give you an estimate of the expectation. Jul 27 '13 at 7:15
• If I recall correctly, similar integrals appear in so-called expectation propagation algorithms for Bayesian inference, and their efficient evaluation is one of the key ingredients. Maybe you can get some ideas from those papers. Jul 27 '13 at 8:41

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.
• I have tried to diagonalize A with a change of variable so $\sum_ {i=1}^n \sum_ {j=1}^n ai,aj xixj$ is just $\sum_{i=1}^n d_i y_i$ The problem is that with this transformation %|x_i|% become $|\sum_{i=1}^n vi yi|$ Jun 20 '13 at 12:10
• $y_i$ are the new integration variables and d and v vectors of real numbers Jun 20 '13 at 12:27