8

Given a multi-dimensional gaussian function, defined by $$f(\boldsymbol{x})=\exp\left\{-\frac{1}{2} \boldsymbol{x}^T\boldsymbol{x} \right\}=\exp\left\{-\frac{1}{2} \sum_{i=1}^nx_i^2 \right\}$$ And an integration region as the form of a $n$-dimensional parallelepiped, defined by $$\mathcal{D} = \left\{\boldsymbol{l \le Lx\le u} \right\}$$ with

  • the lower triangular matrix $\boldsymbol{L}\in \Bbb R^{n\times n}$ where all lower elements and diagonal equal to $1$, all upper elements equal to $0$ $$ \boldsymbol{L} = \left( \begin{matrix} 1&0&0&\ldots&0\\ 1&1&0&\ldots&0\\ 1&1&1&\ldots&0\\ \vdots&\vdots&\vdots&\ddots&0\\ 1&1&1&\ldots&1\\ \end{matrix} \right) $$
  • the vectors $\boldsymbol{l,u} \in \Bbb R^n$: $\boldsymbol{l} = (l_1,...,l_n)'$ and $\boldsymbol{u} = (u_1,...,u_n)'$

The integration region

Are there any methods/algorithms that we can use to approximate the integral of $f(\boldsymbol{x})$ over $\mathcal{D}$ $$\int_{\mathcal{D}}f(\boldsymbol{x})d\boldsymbol{x}=\int_{\{\boldsymbol{l \le Lx\le u} \}}\exp\left\{-\frac{1}{2} \boldsymbol{x}^T\boldsymbol{x} \right\}d\boldsymbol{x}$$ satisfying

  • Fast computation (because later I must compute many integrals with different values of $\boldsymbol{l,u}$)
  • The accuracy doesn't need to be high (absolute error less than $10^{-3}$ is sufficient)

My attempt: we may use Monte Carlo simulation to approximate this integral but given the very specific form of the integration region and also the integrand, I hope there may exist a faster numerical method/algorithm/closed-form approximation.

Besides, we notices that the inversion matrix $\boldsymbol{L}^{-1}$ is an upper bi-diagonal matrix $$ \boldsymbol{L}^{-1} = \left( \begin{matrix} 1&-1&0&\ldots&0\\ 0&1&-1&\ldots&0\\ 0&\ddots&\ddots&\ddots&0\\ \vdots&\ddots&\ddots&\ddots&-1\\ 0&\ldots&\ddots&\ddots&1\\ \end{matrix} \right) $$ So, by making a change of variable $\boldsymbol{y = Lx}$, we can transform the integral into $$\int_{\mathcal{D}}f(\boldsymbol{x})d\boldsymbol{x}= \int_{\{\boldsymbol{l' \le y \le u'}\}} \exp \left\{-\frac{1}{2} \left(y_n^2+\sum_{i=1}^{n-1} (y_i-y_{i+1})^2 \right) \right\} d\boldsymbol{y}$$ with $\mathcal{D}' = \{\boldsymbol{l' \le y \le u'}\}$ is a rectangular region.

I tried to use multivariae Mill ratio but it doesn't work. In fact, multivariae Mill ratio is for the calculation of $\mathbb{P}(\boldsymbol{X \le u})$, not $\mathbb{P}(\boldsymbol{l \le X \le u})$ with $\boldsymbol{X}$ following a multivariate normal distribution $\mathcal{N}(\boldsymbol{\mu,\Sigma})$.

Thank you in advance!

NN2
  • 20,162
  • 1
    have you tried to approximate it with Mill's ratio? see, https://nvlpubs.nist.gov/nistpubs/jres/66B/jresv66Bn3p93_A1b.pdf and https://nvlpubs.nist.gov/nistpubs/jres/68B/jresv68Bn1p3_A1b.pdf – hyportnex Jun 26 '21 at 23:48
  • @hyportnex Thank you for the Mill's ratio! I haven't known this method yet. I'll try it. – NN2 Jun 27 '21 at 00:03
  • 1
    I think you will find this useful: STECK "LOWER BOUNDS FOR THE MULTIVARIATE NORMAL MILLS' RATIO" The Annals of Probability, 1979, Vol. 7, No. 3, pp547-551. As far as I can remember once I used this publication a long, long time ago ... – hyportnex Jun 28 '21 at 00:54
  • @hyportnex Thank you very much. This paper is really useful. I'll use it. As the inverse of my $\Sigma$ is a tridiagonal matrix, the bounds must be easy to compute. – NN2 Jun 28 '21 at 18:06
  • Unluckily, the paper of Steck can't help me. Because it presents the bounds of $\mathbb{P}(\boldsymbol{X} \le \boldsymbol{a})$ and not $\mathbb{P}(\boldsymbol{l}\le \boldsymbol{X} \le \boldsymbol{u})$ which is more difficult. PS: if I apply the result of $\mathbb{P}(\boldsymbol{X} \le \boldsymbol{a})$, I need to decompose $\mathbb{P}(\boldsymbol{l}\le \boldsymbol{X} \le \boldsymbol{u})$ in not less than $2^n$ probabilities of type $\mathbb{P}(\boldsymbol{X} \le \boldsymbol{a})$ – NN2 Jul 08 '21 at 22:57
  • Have you you looked at the jresv66Bn3p93_A1b paper? Savage has both upper and lower exponentially tight bounds for tail events. – hyportnex Jul 09 '21 at 00:57
  • It may be of some help if you could specify whether $|\boldsymbol u|\gg1$, $|\boldsymbol v|\gg1$ or $|\boldsymbol u-\boldsymbol v|\ll1$. – Ѕᴀᴀᴅ Jul 09 '21 at 02:31
  • @Saad Unluckily, $\boldsymbol{l}$ and $\boldsymbol{u}$ can take any value in $\mathbb{R}^n$ – NN2 Jul 09 '21 at 07:42
  • Can you exploit the fact that different $l,u$ have overlaps, so that MC can produce the value of several integration regions in parallel? I mean that a single sample contributes to the integral in all parallelepipeds that contain it. – user619894 Jul 15 '21 at 06:35

0 Answers0