Integrate standard normal pdf over union of intersection of half-spaces

136 Views Asked by At

Let $\varphi(x) = {1 \over \sqrt{2\pi}^n} e^{-\|x\|^2/2}$ be the probability density function for the standard normal distribution in $n$ dimensions. Let $x_1\ge c_1,\ldots,x_n \ge c_n$ be $n$ inequalities in a single variable. Let $S_{\ge k}$ denote the subset of $\mathbb{R}^n$ that satisfies at least $k$ out of $n$ of these inequalities. I would like to compute

$$\int_{x \in S_{\ge k}} \varphi(x) \; dx.$$

Is there an efficient way to compute this integral? Or, failing that, a good way to estimate it?

The best I can come up with is to use the inclusion-exclusion formula. I end up with $\sum_{i=0}^k {n \choose i}$ terms, where each term has the integral $\int_{x \in T} \varphi(x) \; dx$ where $T$ is a subset of $\mathbb{R}^n$ that satisfies some specific subset of (at most $k$) specified inequalities and does not satisfy the rest. This kind of integral can be computed in closed form (using the error function erf), so each term can be computed exactly. Unfortunately, the number of terms obtained in this way can be exponential in $k$. Is there a better way?

1

There are 1 best solutions below

0
On

This problem reduces to the following discrete problem:

Given independent Bernoulli variables $X_1,\dots,X_n$ where $X_i \sim \text{Bernoulli}(p_i)$, compute the probability that $X_1 + \dots + X_n \ge k$.

In particular, we let $X_i$ be 1 if $x_i \ge c_i$, or 0 otherwise; then it is easy to compute $p_i=\Pr[X_i=1]$ in closed form (using the error function erf).

Now the discrete problem can be solved using generating functions. Let $Y=X_1 + \dots + X_n$, let $f_i(x) = (1-p_i) + p_i x$ be the generating function of $X_i$, and $f(x) = \sum_i \Pr[Y=i] x^i$ be the generating function of $Y$. Then

$$f(x) = f_1(x) \cdots f_n(x).$$

Moreover, we can compute $f(x)$ efficiently. Then, $\Pr[Y \ge k]$ can be read off from the coefficients of $f(x)$.

The entire procedure can be done efficiently: in $O(n \log^2 n)$ time, using FFT. It can even be computed in $O(n \log^2 k)$ time by keeping the low $k$ coefficients, using them to compute $\Pr[Y < k]$, and then using $\Pr[Y \ge k] = 1 - \Pr[Y < k]$.

My thanks to Michael Hardy for suggesting the key insight that this can be reduced to a discrete problem.