Integrating density over measurable sets yields probability? I have a bit of trouble proving a statement expressing the expectation of a function of random variables as the conditional expectation, namely:
$$\mathbb{E}g(X,Y) = \int_{\mathbb{R}^n} \int_{\mathbb{R}^m} g(x,y) f(y|x) \, dy \cdot f_X(x) \, dx = \mathbb{E}(\mathbb{E}[g(X,Y)|X])$$
with $X,Y$ appropriately dimensioned real-valued random variables, $f(y|x) := \frac{f_{X,Y}(x,y)}{f_X(x)}$ being the conditional density of $f_{X,Y}$ with respect to $(X,Y)$ and $f_{X}$ to $X$.
There is a statement leading up to this claim in my notes:
$$\forall A \in \mathcal{B}^n, B \in \mathcal{B}^m: \mathbb{P}(X \in A, Y \in B) = \int_{A}\int_{B} f(y|x) \, dy f_{X}(x) \, dx$$
I have already seen and proven a similar statement for a single random variable, but how does it work out in the multivariate case? Can it also be done using the fact that $\mathcal{B}^i$ is generated by half open box sets? I played around with the idea and using the integral representation of the densities cdf's, but the fact that there are multiple variables to be considered leaves me stumped.
As a bonus: How do any of the equalities in the first paragraph make any sense?? They are stated as is in my notes without any further exposition. I know all of the listed quantities and a couple of identities for the conditional expectation.
Thank you for reading!
 A: The comment of @NapD.Lover provides a useful hint. Let $X:\Omega \to \mathbb{R}^n,Y:\Omega \to \mathbb{R}^m$ be measurable random variables. Let $A \in X^{-1}(\mathcal{B}(\mathbb{R}^n))=:\sigma(X)$. Then if $A \in \sigma(X),\,\exists B \in \mathcal{B}(\mathbb{R}^n):A=X^{-1}(B)$. Let $E[|g(X,Y)|]<\infty$. So, to prove that $h(X)=\int_{\mathbb{R}^m}g(X,y)f_{Y|X}(X,y)dy$ is s.t. $h(X)=E[g(X,Y)|X]$ a.e. we use:
$$\begin{aligned}
\int_Ag(X,Y)dP&=\int_{X^{-1}(B)}g(X,Y)dP=\\
&=\int_{\Omega}\mathbf{1}_{X^{-1}(B)}(\omega)g(X(\omega),Y(\omega))P(d\omega)\\
&=\int_{\mathbb{R}^n\times \mathbb{R}^m} \mathbf{1}_{B\times \mathbb{R}^m}(x,y)g(x,y)P_{X,Y}(dx,dy)=\\
&=\int_{\mathbb{R}^n}\int_{\mathbb{R}^m} \mathbf{1}_{B\times \mathbb{R}^m}(x,y)g(x,y)f_{X,Y}(x,y)dxdy=\\
&=\int_{\mathbb{R}^n}\int_{\mathbb{R}^m} \mathbf{1}_{B\times \mathbb{R}^m}(x,y)g(x,y)f_{Y|X}(x,y)f_X(x)dxdy=\\
&=\int_{B}\bigg(\int_{\mathbb{R}^m} g(x,y)f_{Y|X}(x,y)dy\bigg)f_X(x)dx=\\
&=\int_{X^{-1}(B)}h(X)dP=\\
&=\int_Ah(X)dP
\end{aligned}$$
where we used Fubini-Tonelli for the measurability and integrability of $h(X)$ and for the double integral representation. The claim follows from $E[g(X,Y)]=E[E[g(X,Y)|X]]=E[h(X)]$. To prove the first equality, we use the tower property and the fact that $\sigma(X)\supseteq \{\emptyset,\Omega\}$ and $E[g(X,Y)|\{\emptyset,\Omega\}]=E[g(X,Y)]$.
