Statement
If
is a zero-mean multivariate normal random vector, then
where the sum is over all the pairings of
, i.e. all distinct ways of partitioning
into pairs
, and the product is over the pairs contained in
. [5] [6]
More generally, if
is a zero-mean complex-valued multivariate normal random vector, then the formula still holds.
The expression on the right-hand side is also known as the hafnian of the covariance matrix of
.
Odd case
If
is odd, there does not exist any pairing of
. Under this hypothesis, Isserlis's theorem implies that
This also follows from the fact that
has the same distribution as
, which implies that
.
Even case
In his original paper, [7] Leon Isserlis proves this theorem by mathematical induction, generalizing the formula for the
order moments, [8] which takes the appearance

If
is even, there exist
(see double factorial) pair partitions of
: this yields
terms in the sum. For example, for
order moments (i.e.
random variables) there are three terms. For
-order moments there are
terms, and for
-order moments there are
terms.
Example
We can evaluate the characteristic function of gaussians by the Isserlis theorem:
Proof
Since both sides of the formula are multilinear in
, if we can prove the real case, we get the complex case for free.
Let
be the covariance matrix, so that we have the zero-mean multivariate normal random vector
. Since both sides of the formula are continuous with respect to
, it suffices to prove the case when
is invertible.
Using quadratic factorization
, we get

Differentiate under the integral sign with
to obtain
.
That is, we need only find the coefficient of term
in the Taylor expansion of
.
If
is odd, this is zero. So let
, then we need only find the coefficient of term
in the polynomial
.
Expand the polynomial and count, we obtain the formula. 
Generalizations
Gaussian integration by parts
An equivalent formulation of the Wick's probability formula is the Gaussian integration by parts. If
is a zero-mean multivariate normal random vector, then
This is a generalization of Stein's lemma.
The Wick's probability formula can be recovered by induction, considering the function
defined by
. Among other things, this formulation is important in Liouville conformal field theory to obtain conformal Ward identities, BPZ equations [9] and to prove the Fyodorov-Bouchaud formula. [10]
Non-Gaussian random variables
For non-Gaussian random variables, the moment-cumulants formula [11] replaces the Wick's probability formula. If
is a vector of random variables, then
where the sum is over all the partitions of
, the product is over the blocks of
and
is the joint cumulant of
.
Consider
uniformly distributed on the unit sphere
, so that
almost surely. In this setting, the following holds.
If
is odd, 
If
is even,
where
is the set of all pairings of
,
is the Kronecker delta.
Since there are
delta-terms, we get on the diagonal:
Here,
denotes the double factorial.
These results are discussed in the context of random vectors and irreducible representations in the work by Kushkuley (2021). [12]