Integration over a function of random variables Let $ X=(X_1,\ldots ,X_n)\colon (\Omega',\mathcal{A}',P)\rightarrow (\bigotimes_{i=1}^n\Omega_i,\bigotimes_{i=1}^n\mathcal{A}_i)$ be a random vector and $f\colon (\bigotimes_{i=1}^n\Omega_i,\bigotimes_{i=1}^n\mathcal{A}_i)\rightarrow \overline{\mathbb{R}}$ be an integrable function.
It should hold that (I can prove it for $n=1$, but don't know how to generalize):
$$\int_{\Omega'} f\circ X\,\mathrm{d}P=\int_{\bigotimes_{i=1}^n\Omega_i}f\,\mathrm{d}(P\circ X^{-1})$$
Ok, now what I always took for granted was that I can do the following:
\begin{equation}
\int f(X,Y)\mathrm{d}P = \int f(x,y)\mathrm{d}P_X\mathrm{d}P_Y\quad (1)
\end{equation}
meaning when taking the expectation $\mathbb{E}[f(X,Y)]=\mathbb{E}_X[\mathbb{E}_Y[f(X,Y)]]$.
But how can I rigorously proof this?
I read that the distribution of a random vector of independent variables was given by the product measures of the individual distributions and therefore I could set $P\circ X^{-1}=\bigotimes_{i=1}^nP\circ X_i^{-1}$ at least for independent variables. Then by using Fubini ($f$ was assumed to be integrable) I could have a proof.
Is my reasoning correct? Why can I integrate like I did in (1) e.g. split the integral?
 A: You can prove something more general.
Let $(\Omega, \mathcal{F})$ and $(\Omega', \mathcal{F}')$ be measurable spaces. Let $\mu$ be a measure on $\mathcal{F}$ and let $T: \Omega \to \Omega'$ be an $\mathcal{F}, \mathcal{F}'$-measurable map. Then we have a measure
$$\mu \circ T^{-1}: \mathcal{F}' \to [0, \infty]: A \mapsto \mu(T^{-1}(A))$$
You can then prove that if $h: \Omega' \to \overline{\mathbb{R}}$ is $\mathcal{F}'$-measurable, we have that
$$\int_{\Omega'} h d (\mu\circ T^{-1}) = \int_\Omega h \circ Td \mu$$
in the sense that if one side exists, so does the other and both sides are then equal (for positive functions the equality always holds).
You can prove this by proving it first for indicator functions, then for simple functions by linearity of the integral. Subsequently, use monotone convergence theorem to lift the result to arbitrary positive functions and use the decomposition $h= h^+-h^-$ to get the result for general functions.
Let me know if you need more details for the last part.
