Take the 2-minute tour ×
MathOverflow is a question and answer site for professional mathematicians. It's 100% free, no registration required.

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')
share|improve this question
2  
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 –  Carl Jun 18 '13 at 9:02
2  
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. –  Brendan McKay Jul 27 '13 at 7:15
1  
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. –  Christian Clason Jul 27 '13 at 8:41

1 Answer 1

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

share|improve this answer
    
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|$ –  Javier Jun 20 '13 at 12:10
    
$y_i$ are the new integration variables and d and v vectors of real numbers –  Javier Jun 20 '13 at 12:27
    
Besides in my case A is not always full rank and det A may be zero –  Javier Jun 20 '13 at 12:29
    
@Javier: Sorry, I just kind of put that out there as a hope it could help. I wonder if you could so some of the integrations analytically and the rest numerically; i. e. split off the part of the basis corresponding to the non-full-rank part. But you also had your problem with the absolute value term mentioned above. –  Tom Dickens Jun 21 '13 at 21:52
    
Sorry if i sound rude, it was not my intention @Tom dickens, i really appreciate all comments and your answer really goes in the direction i took to find an analytical solution. –  Javier Jun 26 '13 at 12:46

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

Not the answer you're looking for? Browse other questions tagged or ask your own question.