PeterDonis said:
If x, y, and z do not commute, is there even an operator that corresponds to this?
The position components ##(x,y,z)## of course always commute. It's part of the basic realizations of the Galileo group (or rather it's quantum extension, but that doesn't matter so much here) in terms of the Heisenberg algebra. Writing ##(x_1,x_2,x_3)## it's given by
$$[\hat{x}_j,\hat{x}_k]=0, \quad [\hat{p}_j,\hat{x}_k]=0, \quad [\hat{x}_j,\hat{p}_k]=\mathrm{i} \delta_{jk} \hat{1}.$$
The hydrogen energy eigenfunctions are of course no eigenstates of ##\hat{\vec{x}}##, because ##\hat{\vec{x}}## doesn't commute with the Hamiltonian, of which you want eigenstates. Due to rotational symmetry to get unique energy eigenstates you need to determine the common eigenstates of a complete set of compatible observables, which (from rotational symmetry again) can be chosen as ##(H,\vec{L}^2,L_3)##, leading to the here discussed energy eigenfunctions.
The probability density for position is for any wave function, i.e., also for these energy-eigenfunctions
$$P(\vec{x})=|\psi(\vec{x})|^2.$$
In this case these eigenfunctions take the form
$$u_{n,l,m}(r,\vartheta,\varphi)=R_{nl}(r') \text{Y}_{lm}(\vartheta,\varphi).$$
Of course, the probability distribution for the distance, ##r## of the electron from the proton is given by
$$\tilde{P}(r)=\int_{\mathbb{R}^3} \mathrm{d}^3 x P(\vec{x}) \delta(|\vec{x}|-r),$$
as for any probability distribution following from "coarse graining". You get
$$\tilde{P}(r)=\int_{0}^{\infty} \mathrm{d} r' \int_{S_2} \mathrm{d}^2 \Omega r^{\prime 2} \delta(r'-r) |R_n(r)|^2 |Y_{lm}(\vartheta,\varphi)|^2 = r^2 |R_{nl}(r)|^2.$$