Intersection of random line segments in the plane Let a point on the plane be randomly chosen via $(\sqrt{\frac{t}{1-t}}\cos(2\pi\theta),\sqrt{\frac{t}{1-t}}\sin(2\pi\theta))$, where $t$ and $\theta$ are uniformly randomly chosen on $[0,1]$ (equivalently, choose a point uniformly randomly on the surface of the sphere and then project stereographically). Then, what is the probability that two random line segments (determined by their endpoints) will intersect?
This is a repost of a subproblem in a previous post that never got answered. Monte Carlo simulation suggests that the answer is precisely $1/5$, but I have no fruitful ideas left how to prove it.
 A: This is not a finished solution, just a collection of ideas, but with a bit of luck it will get you there.


*

*Switch to Cartesian coordinates. Expressing intersections there will be easier. To achieve this, you need a probability density function $p(x,y)$. It should be proportional to the ratio of sphere surface area per plane surface area after stereographic projection, for infinitesimally small areas. It should only depend on the (squared) radius $x^2+y^2$. And of course it should sum up to one, as in $$\iint_{-\infty}^{+\infty}p(x,y)\,\mathrm dx\,\mathrm dy=1$$ Unless I made a mistake, the probability density function you want should be $$p(x,y)=\frac1{\pi\left(x^2+y^2+1\right)^2}$$ This is based not on your formula for $t$ but on my considerations for stereographic projection of the unit sphere onto the equatorial plane. Please double-check this.

*With probability $1$ any three random points do not lie on a line. In that case you can express the fourth point as a linear combination of these, namely $$P_4=\lambda_1P_1+\lambda_2P_2+\lambda_3P_3\qquad\text{with }\lambda_1+\lambda_2+\lambda_3=1$$ Then segment $(P_1,P_2)$ will intersect segment $(P_3,P_4)$ iff $\lambda_1>0,\lambda_2>0,\lambda_3<0$.

*Combine these. Three points are random in the plane, the fourth random but satisfying these constraints. \begin{align*}
f_1&=\iint_{-\infty}^{+\infty}f_2\,p(x_1,y_1)\,\mathrm dx_1\,\mathrm dy_1 \\
f_2&=\iint_{-\infty}^{+\infty}f_3\,p(x_2,y_2)\,\mathrm dx_2\,\mathrm dy_2 \\
f_3&=\iint_{-\infty}^{+\infty}f_4\,p(x_3,y_3)\,\mathrm dx_3\,\mathrm dy_3 \\
f_4&=\int_0^1\int_{1-\lambda_1}^{+\infty} q(x_4,y_4)\,\mathrm d\lambda_2\,\mathrm d\lambda_1
+\int_1^{+\infty}\int_0^{+\infty} q(x_4,y_4)\,\mathrm d\lambda_2\,\mathrm d\lambda_1 \\
x_4 &= \lambda_1x_1+\lambda_2x_2+(1-\lambda_1-\lambda_2)x_3 \\
y_4 &= \lambda_1y_1+\lambda_2y_2+(1-\lambda_1-\lambda_2)y_3
\end{align*}

*The formulation above is using a different probability density function $q$ in the last step, due to the different parametrization. You will need to express $q$ in terms of $p$, using regular rules for integration by substitution. Essentially $\mathrm dx\,\mathrm dy$ describes a rectangular area in the plane. The corresponding area $\mathrm d\lambda_1\,\mathrm d\lambda_2$ projects onto the plane as a parallelogram-shaped area which you can compute as a the absolute value of a determinant $$q(x,y)=\left\lvert\det\begin{pmatrix}x_1-x_3&x_2-x_3\\y_1-y_3&y_2-y_3\end{pmatrix}\right\rvert\,p(x,y)$$ You can of course move that determinant in front of the integral.

*Now you only have to hope that some combination of brain power and computer algebra system can compute these integrals without building up too much complexity.
By now I gave that a try, and the results are rather discouraging. Seems the terms will become already quite complicated at the innermost integral, and Sage asks me for complicated case distinctions. So I'm no longer optimistic this will be a suitable approach without any fundamental insights into the computation of the integrals.
