  We prove a formula for the evaluation of expectations containing a scalar
function of a Gaussian random vector multiplied by a product of the random
vector components, each one raised at a non-negative integer power. Some powers
could be of zeroth-order, and, for expectations containing only one vector
component to the first power, the formula reduces to Stein's lemma for the
multivariate normal distribution. On the other hand, by setting the said
function inside expectation equal to one, we easily derive Isserlis theorem and
its generalizations, regarding higher order moments of a Gaussian random
vector. We provide two proofs of the formula, with the first being a rigorous
proof via mathematical induction. The second is a formal, constructive
derivation based on treating the expectation not as an integral, but as the
consecutive actions of pseudodifferential operators defined via the
moment-generating function of the Gaussian random vector.
