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]