I got used to treat transformation of random variables in the following way. Example: let $x, y \sim P_{\{x,y\}} (x, y)$. Then if $z = f(x, y)$, the probability distribution for $z$ would be*
$$ P_z (z) = \int P (x, y) \, \delta (z - f (x, y)) $$
So I have a following problem: I have a distribution $P (r)$ for magnitudes of vector $\vec{r}$ and also $P_\theta (\theta)$ for the angle between $\vec{r}$ and $\hat{z}$, i.e. $\cos \theta = z / r$ and $r$ and $\theta$ are uncorrelated. If this is true, then $z = r \cos \theta$ and the PDF for $z$-component only is
$$ P_z (z) = \int \limits_0^\infty \mathrm{d} r \, P_r (r) \int \limits_0^\pi \mathrm{d} \theta \, P_\theta (\theta) \, \delta (z - r \cos \theta) $$
This can be simplified further by letting $u = \cos \theta$ and integrating over $\theta$ using the $\delta$-function, but that's not important now. In this integral, everything is clear to me.
However, when I try to invert this logic and get the distribution for $\theta$ from knowing distributions for $z$ and $r$, I run into problems. Naively, the formula is as follows
$$ P_\theta (\theta) = \int \limits_0^\infty \mathrm{d} r \, P_r (r) \int \limits_{-\infty}^\infty \mathrm{d} z \, P_z(z) \, \delta (\theta - \cos^{-1} (z/r)) $$
Problem: the argument of the $\delta$-function is complex for $z > r$. That makes sense: in everyday geometry, $|z| < r$, because no component can be bigger than the magnitude of the vector, however, statistically, there is nothing that would prevent to have $|z| > r$ in this integral: in fact, if $P_r (r)$ doesn't have a sharp cutoff, $P_z (z)$ extends to infinity as well (can be seen from the expression for $P_z (z)$ above).
So where is the catch in all this? Am I supposed to just restrict the integral over $z$ between $-B$ and $B$ so that the expression inside the $\delta$-function is real only? Or, there might be a mistake in my assumptions: while $r$ and $\theta$ can be statistically independent (as it happens in isotropic systems), $r$ and $z$ are always correlated, so $P (r, z)$ cannot be separated and must be treated together, thus $|z| < r$ would be true by the nature of the function of $P (r, z)$ (it would be zero if $|z| > r$)? If this is true, how is the joint distribution for $r$ and $z$ obtained, given the distributions for $r$ and $\theta$ (uncorrelated)?
*Derivation: the CDF for $z$ is given by $$ F_z (z) = \int P (x, y) \, \theta(f (x, y) - z) $$ because of the simple counting: $\theta$-function chooses the part of the $x,y$ space that satisfies $z < f (x, y)$. The formula for the PDF is obtained by differentiating w.r.t. $z$
