Consider an arbitrary state: $$|\psi\rangle = a|0\rangle+b|1\rangle,$$
where $a=cos\left(\frac{\theta}{2}\right), b=sin\left(\frac{\theta}{2}\right)e^{i\phi}$ (neglecting global phase), $\phi$ is the polar angle, $\theta$ is the azimutal angle on the Bloch sphere.
The Bloch sphere vectors can be found as: $$r_z=cos\left(\theta\right), r_x=sin\left(\theta\right)cos\left(\phi\right), r_y=sin\left(\theta\right)sin\left(\phi\right)$$
If I want to project this state $|\psi\rangle$ into the state $|1\rangle\langle1|$, I can get $p_z$ component of the Bloch vector:
$$p_z = |\langle1|\psi\rangle|^2 = |\langle1|\left(a|0\rangle+b|1\rangle\right)|^2 = |b|^2 = sin^2\left(\frac{\theta}{2}\right) = \frac{1-cos(\theta)}{2} = \frac{1-r_z}{2}$$
To find another two expectation values $p_x$ and $p_y$ I need to change the basis of initial state (which is $|0\rangle$) by rotating the Bloch sphere vector by $\frac{\pi}{2}$ around $X$ or $Y$ axis of the Bloch sphere.
Here is the derivation for the expectation value $p_x$ (after rotating the Bloch sphere vector by $\frac{\pi}{2}$ around $Y$ axis):
X-axis basis states:$$ |+\rangle = \frac{|0\rangle+|1\rangle}{\sqrt{2}}, |-\rangle = \frac{|0\rangle-|1\rangle}{\sqrt{2}},$$then
$$\psi_x = a|+\rangle+b|-\rangle$$,
projecting into the state $|1\rangle\langle1|:$
$$p_x = |\langle1|\psi_x\rangle|^2 = \left(|a\frac{1}{\sqrt{2}}+b\frac{1}{\sqrt{2}}|\right)^2 = \frac{1}{2}\left(|cos\left(\frac{\theta}{2}\right)-sin\left(\frac{\theta}{2}\right)cos\left(\phi\right)-isin\left(\frac{\theta}{2}\right)sin\left(\phi\right)|\right)^2,$$
after expanding the module and finding the square:
$$p_x = \frac{1+r_x}{2}$$
For $p_y$ I get almost the same result:
Y-axis basis states:
$$|R\rangle = \frac{|0\rangle+i|1\rangle}{\sqrt{2}}, |L\rangle = \frac{|0\rangle-i|1\rangle}{\sqrt{2}},$$
repeating the same routine as for $p_x$:
$$\psi_y = a|R\rangle+b|L\rangle$$,
$$p_y = |\langle1|\psi_y\rangle|^2 = \left(|ai\frac{1}{\sqrt{2}}+b(-i)\frac{1}{\sqrt{2}}|\right)^2 = \frac{1}{2}\left(|cos\left(\frac{\theta}{2}\right)i-sin\left(\frac{\theta}{2}\right)cos\left(\phi\right)i+sin\left(\frac{\theta}{2}\right)sin\left(\phi\right)|\right)^2,$$
$$p_y = \frac{1-r_x}{2}$$
The result is strange, because I expected to get $\frac{1-r_y}{2}$
I think I made the mistake somewhere in basis definition, but I can't understand, where could it be.
Here's the explanation:
You have $|R\rangle=\frac{|0\rangle+i|1\rangle}{\sqrt 2}$, $|L\rangle=\frac{|0\rangle-i|1\rangle}{\sqrt 2}$ and you just created a state $\psi_y=a|R\rangle+b|L\rangle$, then you project it into $|1\rangle\langle 1|$ . But what does the $z$ basis look like from the $y$ basis? Let's see!
We can tell that $|0\rangle=\frac{|R\rangle+|L\rangle}{\sqrt 2}$ and $|1\rangle=\frac{|R\rangle-|L\rangle}{i\sqrt 2}$. This doesn't look like the form of the $y$ basis... instead, it looks exactly like the $x$ basis with respect to the $z$ basis, except for the global phase $i$ which doesn't do anything in the probability!
Your angles are defined with respect to your $y$ basis, and from the point of view of the $y$ basis, the $z$ basis looks like the $x$ basis relatively. So it's expected that you get $r_x$ again.
If you do the same thing to the $x$ basis, you'll notice that your $z$ basis looks like $x$ basis from $x$ basis as well. So where does the $r_y$ go, though? Well, if we add a global phase $i$ on the $|L\rangle$, i.e, $|L\rangle=\frac{-i|0\rangle-|1\rangle}{\sqrt 2}$, then the $z$ basis would look like $y$ basis from $y$, and we'll get the $r_y$ term. The same thing can be done for $x$ basis. But I'm not sure if there's any point in doing either of that...