I've seen some post addressing this but still I didn't find the answer I looking for, if there's some already existing question where this is discussed in detail please let me know.
Let $(\Omega,\mathscr F,P)$ be a probability space and let $(\mathbb R,\mathscr B,\lambda)$ the a measure space where the $\sigma$-algebra is the Borel $\sigma$-algebra, and $\lambda$ stands for the Lebesgue measure.
Let be $X:\Omega\rightarrow \mathbb R$ be a $\mathscr F/\mathscr B$-measurable function, aka. a random variable.
My aim is to prove rigorously that
$$\int_{\Omega}XdP=\int_{\mathbb R}Xfd\lambda$$ where $f$ is the density function (assume for now that it exists).
I have that the distribution of my random variable $X$ can be written as the image measure of $P$ under $X$. $$P(X\in B)=P(X^{-1}(B)), \forall B\in\mathscr B$$
Assume furthermore that $P\circ X^{-1}<<\lambda$ and hence we can find the Radon Nikodym derivative and call it the density function $f$ of $X$.
Now in this set of slides it's stated that (I am changing the notation to make it compatible with mine):
Let be $P\circ X^{-1}:=\nu=f\cdot\lambda $ be a measure on $(\mathbb R,\mathscr B)$, then forall $g\in\mathcal M(\mathscr B)$ it holds that if $$\int |gf|d\lambda<\infty$$
then $$\int gd\nu=\int g P\circ X^{-1} (dx)=\int gf\lambda(dx)$$ (I assume from now on that the condition holds)
I can take $g(x)=1_Bx$,$B\in\mathscr B$ which is a Borel measurable function
$$\int_B x P\circ X^{-1} (dx)=\int_B xf\lambda(dx)$$ But now I can see $x$ as a particular realization of $X$, i.e. $x=X(\omega)$
$$\int_{X^{-1}(B)} X(\omega) P\circ X^{-1} (dX(\omega))=\int_B xf\lambda(dx)$$ Now if I let $B=\mathbb R$
$$\int_{\Omega} X(\omega) P\circ X^{-1} (dX(\omega))=\int_{\mathbb R} xf\lambda(dx)$$
And there is where I get stucked, I don't really know how to go from the integral w.r.t. the image measure to the integral w.r.t the probability measure. I feel that I might be close to the solution (or completely off). Any hint or suggestion will be greatly appreciated.
After some reading and thanks to the Kabo Murphy's answer I came up with an answer to my doubt, I am going to post it here in case someone could find this helpful.
We have a random variable $X:\Omega\rightarrow \mathbb R$ then one could define a measure in $(\mathbb R,\mathscr B)$ as $$P(X^{-1}(B))=(PX^{-1})(B), B\in\mathscr B$$ (from now on $B$ is some Borel measurable set).
Take the left hand side and notice that it's by definition equal to $$P(X^{-1}(B))=\int 1_{X^{-1}(B)}(\omega)P(d\omega)=\int 1_{B}\big(X(\omega)\big)P(d\omega)$$
Now the right side it's equal by definition to
$$(PX^{-1})(B)=\int 1_{B}(x)(PX^{-1})(dx)$$
Then we have that $$\int 1_{B}\big(X(\omega)\big)P(d\omega)=\int 1_{B}(x)(PX^{-1})(dx)$$
This will hold for any simple function $f_n$ by linearity, and if we assume that $0\leq f_n\uparrow f$ then $0\leq f_n(X) \uparrow f(X)$ we can apply Beppo-Levi theorem to show that this will hold also for the limit
$$\int f(X(\omega))P(d\omega)=\int f(x)(PX^{-1})(dx)$$
If we take an f that is not nonnegative the previous must hold for $|f|$ and then if $f\circ X$ is integrable w.r.t $P$ then $f$ is integrable w.r.t. $(PX^{-1})$.
Now assume that $f$ is the identity map, $f(x)=x$ and assume that $E[X]<\infty$ then we have that
$$\int X(\omega)P(d\omega)=\int x (PX^{-1})(dx).$$
Now assume that $(PX^{-1})$ has a density $\delta$ w.r.t. the Lebesgue measure $\lambda$, since $f:x\mapsto x$ is integrable w.r.t. $(PX^{-1})$ (by the previous assumption of $E[X]<\infty$) then $\delta x$ is integrable w.r.t. $\lambda$, and the integrals must be equal.
Finally putting everything together we have
$$\int_{\Omega} X(\omega) P(d\omega)=\int_{\mathbb R} x (PX^{-1})(dx)=\int_{\mathbb R} x\delta\lambda(dx)$$