Using integrals to prove that the mean of the sampling distribution is the population mean Let the random variables $X_1, X_2, \dots X_n$ denote a random sample from a population.
The sample mean of these random variables is:
$\overline{X}=\frac{1}{n}\sum\limits_{i=1}^{n}X_i$
I would like to show that the mean of the sampling distribution of the sample mean is $\mu$, the population mean.
Here's what I have done:

$$\begin{align} 
E(\overline{X}) &= \int\limits_{\overline{X}} \bar{x} f(\bar{x})\,\, d\bar{x} \\
&=\int\limits_{\overline{X}} \left(\frac{1}{n} \sum\limits_{i=1}^n X_i \right) f(\bar{x}) \, d\bar{x}
\end{align}$$

From here, I am not sure what to do anymore but anyway I end up with:

$$\begin{array} {cc} 
&=& \frac{1}{n} \left( \int\limits_{\overline{X}}X_1f(\bar{x}) \, d\bar{x} + \int\limits_{\overline{X}}X_2f(\bar{x}) \, d\bar{x}+ \dots + \int\limits_{\overline{X}}X_nf(\bar{x}) \, d\bar{x}\right)
\end{array}$$

Now, I don't know how to complete this as I am unsure how to interpret the last equation. Somehow, the $\int\limits_{\overline{X}}X_if(\bar{x}) \, d\bar{x}$ is suppose to equal to $\mu$ but I don't see how that can be true.
I know the answer will be $\mu$ because of here but I would like to arrive at the answer using integrals instead.
 A: You should not confuse the argument to the probability density function with the random variable.  Often one uses a lower-case letter for the former and a capital letter for the latter.  For example, if the density function for the random variable $X$ is $f$, then one can speak of $\int_{-\infty}^\infty xf(x)\,dx$, and the lower-case "$x$" is not the random variable $X$.  It is analogous to such things as "$P(X\le x)$": the "$X$" and the "$x$" mean two different things.  One can also write $f(3)$, and it doesn't mean the density function of a random variable called "$3$"; it means the value at the number $3$, of the random variable $X$.
Thus one can write $f_X(4)$ and $f_Y(4)$ and they're the values at $4$ of the density functions of two different random variables $X$ and $Y$.
If you're writing
$$
\mathbb E(\bar X) = \int_{-\infty}^\infty \bar x f(\bar x)\,d\bar x,
$$
then that means the same thing as
$$
\mathbb E(\bar X) = \int_{-\infty}^\infty x f_{\bar X}(x)\,dx.
$$
You can't put the random variable $\bar X$ in place of the bound variable $\bar x$.
To understand this, it may help to realize that
$$
\sum_{j=1}^3 (i^2 j^3) \text{ and } \sum_{k=1}^3 (i^2 k^3)
$$
both mean
$$
i^2 1^3 + i^2 2^3 + i^2 3^3,
$$
and thus both depend on the value of the "free variable" $i$ but not on any value of the "bound variable" that is either $j$ or $k$.  There's no $j$ or $k$ for them to depend on.  One can freely change the name of the bound variable from $j$ to $k$ without changing the value of the expression.
The same thing applies to
$$
\mathbb E(X) = \int_{\mathbb R} x f_X(x)\,dx = \int_{\mathbb R} w f_X(w)\,dw.
$$
One can freely change the name of the bound variable from $x$ to $w$.  But the capital $X$ still refers to the same random variable.
If you intend $f$ to be the density function of $\bar X$, then in your second displayed line in your question, $f$ is still the density function of the random variable $\bar X$.
If I wanted to do this by using the density function of the sample mean $\bar X$, I'd actually need to show how that function depends on the density function of $X_1$.  That would be an $n$-gold convolution.
I wouldn't do it that way if I could help it.  Instead I'd use the linearity of expectation, thus:
$$
\mathbb E\left(\frac{X_1+\cdots+X_n}{n}\right) = \frac1n\mathbb E(X_1+\cdots+X_n) = \frac1n((\mathbb E X_1)+\cdots+(\mathbb E X_n))
$$
There is of course the problem of how to prove the linearity of expectation.  Maybe that's posted here somewhere as a separate question already.
A: No  problem giving a complete answer. Let's see
As already stated, the sample mean is a function of many random variables, and so the symbol $E$ refers to the expected value with respect to their joint distribution. Denoting $\mathbf X$ the multivariate vector of the $n$ r.v.'s, their joint density can be written as $f_{\mathbf X}(\mathbf x)= f_{X_1,...,X_n}(x_1,...,x_n)$ and their joint support 
$D = S_{X_1} \times ...\times S_{X_n}$
The sample mean is a function of this multivariate vector, $\bar X = \frac 1n \sum_{i=1}^{n}X_i = g(\mathbf X)$. Using the Law of Unconcscious Statistician  We have
$$E[\frac 1n \sum_{i=1}^{n}X_i] = \int_D g(\mathbf x)f_{\mathbf X}(\mathbf x)d\mathbf x$$.
Under convergence regularity conditions we can decompose the multidimensional integral into an n-iterative integral: 
$$E[\frac 1n \sum_{i=1}^{n}X_i] = \int_{S_{X_n}}...\int_{S_{X_1}}\left[\frac 1n \sum_{i=1}^{n}x_i\right]f_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_n $$
and using the linearity of integrals we can decompose into
$$ = \frac 1n\int_{S_{X_n}}...\int_{S_{X_1}}x_1f_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_n \; + ...\\ ...+\frac 1n\int_{S_{X_n}}...\int_{S_{X_1}}x_nf_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_n $$
For each n-iterative integral we can re-arrange the order of integration so that, in each, the outer integration is with respect to the variable that is outside the joint density. Namely,
$$\frac 1n\int_{S_{X_n}}...\int_{S_{X_1}}x_1f_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_n = \\\frac 1n\int_{S_{X_1}}x_1\int_{S_{X_n}}...\int_{S_{X_2}}f_{X_1,...,X_n}(x_1,...,x_n)dx_2...dx_ndx_1$$ 
and in general
$$\frac 1n\int_{S_{X_n}}...\int_{S_{X_j}}...\int_{S_{X_1}}x_jf_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_j...dx_n =$$
$$=\frac 1n\int_{S_{X_j}}x_j\int_{S_{X_n}}...\int_{S_{X_{j-1}}}\int_{S_{X_{j+1}}}...\int_{S_{X_1}}f_{X_1,...,X_n}(x_1,...,x_n)dx_1...dx_{j-1}dx_{j+1}......dx_ndx_j$$
As we calculate one-by-one the integral in each n-iterative integral (starting from the inside), we "integrate out" a variable and we obtain in each step the "joint-marginal" distribution of the other variables. Each n-iterative integral therefore will end up as $\frac 1n\int_{S_{X_j}}x_jf_{X_j}(x_j)dx_j$.
Bringing it all together we arrive at 
$$E[\frac 1n\sum_{i=1}^{n} X_i ] = \frac 1n\int_{S_{X_1}}x_1f_{X_1}(x_1)dx_1 +...+\frac 1n\int_{S_{X_n}}x_nf_{X_n}(x_n)dx_n $$
But now each simple integral is the expected value of each random variable separately, so 
$$= E[\frac 1n\sum_{i=1}^{n} X_i ] = \frac 1nE(X_1) + ...+\frac 1nE(X_n) =  \frac 1nE(X) + ...+\frac 1nE(X)$$
$$= \frac 1n nE(X) = E(X)$$
