Parametric approximation of the expectation of a function Given a function $f(x):R^{n} \to R$ where $x \in R^n$, and $p(x)=\mathcal{N}(x|m,V); V=diag(v)$ is a multivariate Gaussian distribution with diagonal covariance matrix $V$. We want to compute $E_{p}f(x)$ based on $m,V$ and function $f$.
One possible way is to use Taylor expansion around $m$, i.e $f(x)=f(m)+\nabla f(x)_{x=m}(x-m) $, then compute:$E_{p}f(x) \approx \int \big( f(m)+\nabla f(x)_{x=m}(x-m)  \big) p(x)dx$.
However, Taylor expansion is accurate when $x$ is close to $m$. When taking expectation, each dimension of $x$ runs from $-\infty $ to $\infty$, leading $x$ is probably far from $m$.
How to estimate $E_p f(x)$ more accurate and we still obtain a parametric form (the estimation is a function of $m$ and $V$) ?
 A: You can use Laplace's method, an introduction here: https://en.wikipedia.org/wiki/Laplace%27s_method
With your notation, $f$ is a function from $\mathbb{R}^n \mapsto \mathbb{R}$, and I will assume it is a positive function.
$$
  E_p f(X) = \idotsint_{\mathbb{R}^n} e^{-g(x)} p(x) \; dx
$$
where $g(x)=-\log f(x)$ (so here we do need that $f$ is positive ...).  Then inserting an expression for the multivariate normal density $p(x)$, and taylor-expanding $g$, around the minimizer of $g$, $x_0$:
$$
  g(x)=g(x_0)+g'(x_0)'(x-x_0) + \frac12 (x-x_0)'g''(x_0)(x-x_0) + \dots
$$
(and then forgetting about the dots for the approximation).  This becomes
$$
   E_p f(X) \approx \idotsint_{\mathbb{R}^n} \exp(-g(x_0)-g'(x_0)'(x-x_0)-\frac12 (x-x_0)'g''(x_0) (x-x_0) ) \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \exp(-\frac12 (x-v)'V^{-1} (x-v)) \; dx
$$
here $g'$ is the gradient, $g''$ is the hessian matrix, and elsewhere ' denotes differentiation. Putting together the two exponentials, this becomes
$$
\frac1{(2\pi)^{n/2} |V|^{(1/2)}} \idotsint_{\mathbb{R}^n} \exp(-g(x_0)-g'(x_0)'(x-x_0)-\frac12 (x-x_0)'g''(x_0) (x-x_0) -\frac12 (x-v)'V^{-1} (x-v)) \; dx
$$
Now, to evaluate the integral we need to "complete the square":
In the following, we concentrate only on the argument of the exponential in the integral. ultiply out everything in the argument, and write the resulting form as 
$$d+a'x -\frac12 x'(V^{-1}+g''(x_0))x
$$
and then require to write it in the form
$$
-\frac12 (x-\mu)'(V^{-1}+g''(x_0)) (x-\mu)
$$
we find that
$d=-\frac12 \mu' (V^{-1}+g''(x_0)) \mu$,
$a' = -\mu'  (V^{-1}+g''(x_0)) $ by solving that we find that 
$\mu= -a (V^{-1}+g''(x_0))^{-1}$. Then, our integral has become
$$
  \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \idotsint_{\mathbb{R}^n} 
     \exp( -\frac12 (x-\mu)' A^{-1} (x-\mu) ) \; dx
$$
where we have introduced $A= (V^{-1}+g''(x_0))^{-1}$. Now, that has the form as the integral of a multivariate normal density, so can be found to be
$$
  \frac1{(2\pi)^{n/2} |V|^{(1/2)}} \cdot (2\pi)^{n/2} |A|^{1/2}
$$
which I leave for you to simplify further. 
