5
$\begingroup$

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')
$\endgroup$
3
  • 2
    $\begingroup$ 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 $\endgroup$
    – user25199
    Commented Jun 18, 2013 at 9:02
  • 2
    $\begingroup$ 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. $\endgroup$ Commented Jul 27, 2013 at 7:15
  • 1
    $\begingroup$ 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. $\endgroup$ Commented Jul 27, 2013 at 8:41

1 Answer 1

1
$\begingroup$

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

$\endgroup$
6
  • $\begingroup$ 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|$ $\endgroup$
    – Javier
    Commented Jun 20, 2013 at 12:10
  • $\begingroup$ $y_i$ are the new integration variables and d and v vectors of real numbers $\endgroup$
    – Javier
    Commented Jun 20, 2013 at 12:27
  • $\begingroup$ Besides in my case A is not always full rank and det A may be zero $\endgroup$
    – Javier
    Commented Jun 20, 2013 at 12:29
  • $\begingroup$ @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. $\endgroup$ Commented Jun 21, 2013 at 21:52
  • $\begingroup$ 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. $\endgroup$
    – Javier
    Commented Jun 26, 2013 at 12:46

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .