Taylor approximation of expected value of multivariate function I have a multivariate function $f:\mathbb{R}^n \to \mathbb{R}$ that takes as input $n$ random variables, each i.i.d. according to a uniform distribution on the interval $(a,b)$. Denote the set of these random variables as $\mathbf{X}.$
I would like to derive an approximate closed-form solution of $E[f(\mathbf{X})]$. 
According to Variance of a function of one random variable, a Taylor series expansion can be done to approximate $E[f(X)]$ for the single variable case:
$$E[f(X)] = f(E[X]) + \frac{f''(E[X])}{2} E[(X- E[X])^2] + R^3$$
where $R^3$ is some remainder term (see link for more details).
Does this extend to the multivariate case as well?
The multivariate taylor for $f(\mathbf{X})$ around $\mathbf{a}$ is
$$f(\mathbf{a})+ (\mathbf{X} - \mathbf{a})^T \nabla f(\mathbf{a}) + (\mathbf{X} - \mathbf{a})^T(H_f(\mathbf{X})^T(\mathbf{X} - \mathbf{a})) + \dots $$
How would the extension for $E[f(\mathbf{X})]$ look? Maybe like this:
$$E[\mathbf{X}] = f(E[\mathbf{X}]) + E[\mathbf{X}- E[\mathbf{X}]]^T\frac{H_f(E[\mathbf{X}])^T}{2} E[\mathbf{X}- E[\mathbf{X}]] + R^3$$
where $E[\mathbf{X}_i] = (a+b)/2, \forall i$.
But doesn't $E[\mathbf{X}- E[\mathbf{X}]]=0$ so the second term is also 0?
How would the extension look and is an extension even possible?
 A: Taylor series approximation of multivariate function $f$ around $x_0$ is
$$ 
f(x) \approx f(x_0) + \nabla f(x_0)'(x-x_0) + \frac{1}{2} (x-x_0)' H_f(x_0) (x-x_0). 
$$
If you substitute $x=X$ and $x_0 = \mathbb{E}X$ you get
$$ 
f(X) \approx f(\mathbb{E}X) +  \nabla f(\mathbb{E}X)'(X-\mathbb{E}X) + \frac{1}{2} (X-\mathbb{E}X)' H_f(\mathbb{E}X) (X-\mathbb{E}X).
$$
Taking expectation on both sides gives
$$
\mathbb{E}f(X) \approx f(\mathbb{E}X) +  \nabla f(\mathbb{E}X)' \mathbb{E}(X-\mathbb{E}X) + \frac{1}{2} \mathbb{E}[(X-\mathbb{E}X)' H_f(\mathbb{E}X) (X-\mathbb{E}X)].
$$
As you noticed, $\mathbb{E}(X-\mathbb{E}X) = 0$, so the expression simplifies to
$$
\mathbb{E}f(X) \approx f(\mathbb{E}X) + \frac{1}{2} \mathbb{E}[(X-\mathbb{E}X)' H_f(\mathbb{E}X) (X-\mathbb{E}X)].
$$
This is as far as you can get without assumptions on X. However in your specific case the second term can be further simplified. Rewriting quadratic form using sums gives
$$
 \mathbb{E}[(X-\mathbb{E}X)' H_f(\mathbb{E}X) (X-\mathbb{E}X)] = \sum_{i=1}^n \sum_{j=1}^n \mathbb{E}[(X_i - \mathbb{E}X_i) H_f(\mathbb{E}X)_{ij} (X_j - \mathbb{E}X_j)] = (*).
$$
If $i \neq j$ then $X_i$ and $X_j$ are independent, therefore
$$
\mathbb{E}[(X_i−\mathbb{E}X_i)Hf(\mathbb{E}X)_{ij}(X_j−\mathbb{E}X_j)] = \mathbb{E}[X_i−\mathbb{E}X_i]Hf(\mathbb{E}X)_{ij}\mathbb{E}[X_j−\mathbb{E}X_j] = 0.
$$
Using that fact, the double sum simplifies to a single sum
$$
(*) = \sum_{i=1}^n \mathbb{E}[(X_i - \mathbb{E}X_i) H_f(\mathbb{E}X)_{ii} (X_i - \mathbb{E}X_i)] = (*).
$$
Expression $H_f(\mathbb{E}X)_{ii}$ is constant (not random), so it can be extracted from expectation
$$
(*) = \sum_{i=1}^n  H_f(\mathbb{E}X)_{ii}\mathbb{E}[(X_i - \mathbb{E}X_i)^2] = \sum_{i=1}^n H_f(\mathbb{E}X)_{ii} Var(X_i).
$$
Summing up, the Taylor series approximation simplifies to
$$
\mathbb{E}f(X) \approx f(\mathbb{E}X) + \frac{1}{2} \sum_{i=1}^n H_f(\mathbb{E}X)_{ii} Var(X_i).
$$
In your case $\mathbb{E}X_i = \frac{a+b}{2}$ and $Var(X_i) = \frac{1}{12}(b-a)^2$. Also, you don't need to compute the whole Hessian matrix, because only its diagonal elements are used in the formula.
