Gaussian integral of a function with non-zero mean (generalizing Wick theorem) From the wikipedia article, for a Gaussian integral of an analytic function we have that

This is equivalent to the Wick theorem when $f(x)$ is a polynomial.
Now I'm trying to obtain a similar formula when there is a linear term in the Gaussian (ie the Gaussian has a nonzero mean).
My guess is that
$$
\int f(x) \exp \left( - \frac{1}{2} x^T A x + B^T x \right) d^n x = \\
= \sqrt{\frac{(2 \pi)^n}{\det A}} \exp \left[ \frac{1}{2} \left( B^T + \frac{\partial}{\partial x_i} \right) A^{-1} \left( B + \frac{\partial}{\partial x_j} \right) \right] f(x) \bigg|_{x=0}
\tag{2}$$
but I can't prove it. Is this equation correct? How can I prove it?
 A: The case with the linear term is obtained from the original one by a simple shift, i.e. the substitution
$$ x = X + A^{-1} B $$
Substitute it to the exponent in your more general integral:
$$ -\frac 12 x^T A x + B^T x = -\frac 12 (X^T+B^T A^{-1}) A(X+A^{-1}B)+B^T (X+A^{-1}B)=\dots $$
I used $A=A^T$. Now, all the terms that are schematically $BX$ i.e. linear in $X$ cancel: the coefficient is $-1/2-1/2+1=0$. The remaining terms give
$$ - \frac 12 X^T A X + \frac 12 B^T A^{-1} B $$
The coefficient $+1/2$ in the second term came from $-1/2+1$. The second term only gives a simple factor (the exponential of that), and it's a part of your result – except that the last $B^T$ in your result should be simply $B$.
The hard, quadratic/Gaussian part of the expression may be rewritten in the Wick way from your first identity. It could be enough if you were satisfied with the evaluation of the $x$-derivatives not at $X=0$ but at the right original value $x=0$ which means, thanks to my substitution
$$ X = -A^{-1}B. $$
However, if you want to use the values of the derivatives at $X=0$, you have to Taylor-expand the shift operator from $X=0$ to $X=-A^{-1}B$. The shift is the operator 
$$ \exp(B^T A^{-1} \frac d{dx}) $$
which is exactly what you get from the mixed terms in your guessed exponent, up to an overall sign perhaps that you will surely be able to catch correctly.
A: Lubos Motl has already provided a correct answer. This answer uses a different approach in the spirit of perturbation theory with $j$-sources:
$$\begin{align}\int_{\mathbb{R}^n} \! d^nx ~f(x)~&e^{-\frac{1}{2}x^TAx +j^Tx}\cr
~=~~& f\left(\frac{\partial}{\partial j}\right)
\int_{\mathbb{R}^n} \! d^nx ~e^{-\frac{1}{2}x^TAx +j^Tx}\cr
~\stackrel{\begin{matrix}\text{Gauss.}\\ \text{int.}\end{matrix}}{=}&C~ f\left(\frac{\partial}{\partial j}\right)e^{\frac{1}{2}j^TA^{-1}j}\cr
~\stackrel{\text{Taylor}}{=}&\left. C~ e^{\left(\frac{\partial}{\partial j}\right)^T
\frac{\partial}{\partial x}}f(x)\right|_{x=0}e^{\frac{1}{2}j^TA^{-1}j}\cr
~=~~&\left. C~ e^{\left(\frac{\partial}{\partial x}\right)^T
\frac{\partial}{\partial j}}e^{\frac{1}{2}j^TA^{-1}j}f(x)\right|_{x=0}\cr
~\stackrel{\text{Taylor}}{=}&\left. C~e^{\frac{1}{2}\left(j+\frac{\partial}{\partial x}\right)^T
A^{-1} \left(j+\frac{\partial}{\partial x}\right)} f(x) \right|_{x=0},\end{align}\tag{A}$$
where the constant
$$C~:=~\sqrt{\frac{(2\pi)^n}{\det A}},\tag{B}$$
where $A$ is a symmetric $n\times n$ matrix with positive-definite real part ${\rm Re}(A)>0$,
where $f$ is an analytic function, and where we have used the Taylor formula
$$f(a+x)~=~ e^{a^T\frac{\partial}{\partial x}}f(x)\tag{C}$$
twice.
Eq. (A) is OP's eq. (2) with sources $j\equiv B$. To obtain OP's first eq. remove the sources.
