Randomly plotting the angle

I've managed to get a closed formula that matches computed values, but it's not pretty. I can share it here, or shall I post it as a solution to your next problem? ;-)

Well, the other solution poster pointed out that:

$\theta = \pi - 2\tan^{-1} (\tfrac yx) ,$ where $x,y$ are independent and identical continuous uniform distribution, $U \sim [0,1]$ . That's all I can gather.

Well, I daren't post it, as we seem to have a discrepancy - any ideas?? I'll have a look at the mode, but meantime I found it useful to look at this on Desmos - that's when I realised it was just a bit of a circle all along. If you paste this...

\arctan\left(\frac{y}{x}\right)+\arctan\left(\frac{1-y}{x}\right)=c\left{0<y<1\right}\left{0<x<1\right}

...into Desmos you should get a slider that lets you play with the graph.

Using the double variable substitution $u =xy$ , $v = y/x$ , we have the Jacobian $\frac{\partial(u,v)}{\partial(x,y)} \; = \; \frac{2y}{x} \; = \; 2v$ and hence the integral $\begin{aligned} J & = \; \int_0^1 \int_0^1 \left(\tan^{-1}\tfrac{y}{x}\right)^2\,dx\,dy \; = \; \tfrac12\int_0^\infty\,dv \int_0^{\mathrm{min}(v,v^{-1})}\,du \,v^{-1}\big(\tan^{-1}v\big)^2 \\ & = \; \tfrac12\int_0^1 \big(\tan^{-1}v\big)^2\,dv + \tfrac12\int_1^\infty v^{-2}\big(\tan^{-1}v\big)^2\,dv \\ & = \; \tfrac12\int_0^1 \left(\big(\tan^{-1}v\big)^2 + \big(\tfrac12\pi - \tan^{-1}v\big)^2\right)\,dv \\ & =\; \int_0^1 \big(\tan^{-1}v\big)^2\,dv - \tfrac12\pi \int_0^1 \tan^{-1}v\,dv + \tfrac18\pi^2 \end{aligned}$ Integrating by parts, $\int_0^1 \tan^{-1}v\,dv \; = \; \Big[v \tan^{-1}v\Big]_0^1 - \int_0^1 \frac{v}{1+v^2}\,dv \; = \; \tfrac14\pi - \Big[\tfrac12\ln(1+v^2)\Big]_0^1 \; = \; \tfrac14\pi - \tfrac12\ln2$ so that $J \; = \; \int_0^1 \big(\tan^{-1}v\big)^2\,dv + \tfrac14\pi \ln2$ Integrating by parts twice more: $\begin{aligned} J & = \; \Big[v\big(\tan^{-1}v\big)^2\Big]_0^1 - \int_0^1 \frac{2v \tan^{-1}v}{1+v^2}\,dv + \tfrac14\pi\ln2 \\ & = \; \tfrac{1}{16}\pi^2 + \tfrac14\pi\ln2 - \Big[\tan^{-1}v\,\ln(1+v^2)\Big]_0^1 + \int_0^1 \frac{\ln(1+v^2)}{1+v^2}\,dv \\ & = \; \tfrac{1}{16}\pi^2 + \int_0^{\frac14\pi}\ln(\sec^2\theta)\,d\theta \; = \; \tfrac{1}{16}\pi^2 - 2\int_0^{\frac14\pi}\ln(\cos\theta)\,d\theta \\ & = \; \tfrac{1}{16}\pi^2 + \tfrac12\pi \ln 2 - G \end{aligned}$ where $G$ is the Catalan constant.

Similarly, but more simply, we have $\begin{aligned} K & = \; \int_0^1 \int_0^1 \tan^{-1}\tfrac{y}{x}\,dx\,dy \; = \; \tfrac12\int_0^\infty\,dv \int_0^{\mathrm{min}(v,v^{-1})}\,du \,v^{-1}\,\tan^{-1}v \\ & = \; \tfrac12\int_0^1 \tan^{-1}v\,dv + \tfrac12\int_1^\infty v^{-2}\tan^{-1}v\,dv \\ & = \; \tfrac12\int_0^1 \left(\tan^{-1}v + \tfrac12\pi - \tan^{-1}v\right)\,dv \\ & =\; \tfrac14\pi \end{aligned}$

BTW Wolfram Alpha may not be able to handle evaluating $J$ , but Mathematica makes short work of it.

Wait... to be honest, I still don't understand your solution.

@Pi Han Goh . I think I made a bad decision to write the angle expression as $\pi-arctan(\frac{y}{x}) - arctan\left(\frac{y}{1-x}\right)$ . Simply writing it as $arccot\left(\frac{x}{y}\right)+arccot\left(\frac{1-x}{y}\right)$ makes more sense.

Thanks! I wonder if it might be useful for the variance question too (though I suspect not - I think the way you're going with integrating over the square is more likely to work). If you've any specific questions let me know and I'll try to make it clearer.