Do you have an idea about how to prove the following? (references will be useful)
Let $x$ a random variable with probability distribution $ρ(x)$, and $y=f(x)$ another random variable, then the probability distribution of $y$ is
$$P(y)=\int \rho (x)\, \delta(y-f(x))\,dx$$
I just haven't seen it in my life, so any help would be really appreciated.
The $\delta$ function, popular in physics, is defined axiomatically by $$ \int_{-\infty}^\infty f(x)\delta(x)\,dx=f(0). $$ In other words, the integral is evaluated by setting the argument to $\delta$ to 0, and making the corresponding substitution in $f(x)$.
But here it seems that instead $\delta(x)=1_{\{0\}}(x)$, the function that is 1 when $x=0$ and 0 otherwise.
It is helpful to use capital letters for random variables. Here $Y=f(X)$, and $$ \int \rho (x)\, \delta(y-f(x))\,dx = \int_{\{x:f(x)=y\}} \rho(x)\,dx = P(f(X)=y) = P(Y=y) $$ This is correct if $X$ is a continuous random variable with density $\rho$, and $Y$ is a discrete random variable with $Y=f(X)$.
It is sort of consistent with the first $\delta$ function definition I mentioned: using that, your integral would evaluate to "$\rho(f^{-1}(y))$", which thinking of $\rho$ as the probability of an infinitesimal region intuitively should be $\sum_{x\in f^{-1}(y)}\rho(x)$, i.e., $\int_{x:f(x)=y}\rho(x)\,dx$.
It also makes sense if $\delta$ is considered formally to be a measure with $\delta(A)=1_{A}(0)$, 1 if $0\in A$, 0 otherwise. Then $\int f(x)\delta(x)dx$ is supposed to mean $\int f(x)d\delta(x)=f(0)$, and according to the definition of the Lebesgue integral, for a partition $\Pi=\{0=a_0<a_1<a_2<\dots\}$, $$ \int \rho(x)d\delta(y-f(x)) = \lim_{\|\Pi\|\rightarrow 0}\sum a_i \delta(\{y-f(x): a_i\le\rho(x)<a_{i+1}\})=\int \rho(x)1_{\{0\}}(y-f(x))dx $$