Average angle between two randomly chosen vectors in a unit square Consider two randomly chosen vectors $(a,b)$ and $(c,d)$ within the unit square, where $a, b, c,$ and $d$ are chosen uniformly from $[0,1]$. What is the expected angle between the vectors?
Here's what I have so far. The angle between any two positive vectors $<a,b>$ and $<c, d>$ is $\arccos{\frac{ac + bd}{\sqrt{(a^2 + b^2)(c^2 + d^2)}}}$. We just need to find the average value of this function over $a, b, c,d$ in range $[0,1]$. This is equivalent to the quadrupal integral $$\iiiint_V \arccos{\frac{ac + bd}{\sqrt{(a^2 + b^2)(c^2 + d^2)}}} \,da\,db\,dc\,dd$$.
I can't find any way to compute this. I entered this into Mathematica and it wasn't able to output even a decimal approximation. I tried all the usual substitions to reduce this into a double integral but none of them worked. Converting to polar coordinates didn't work out nicely either. For one, the bounds become harder to work with and the jacobian is pretty nasty, so converting doesn't seem like it would help me. Probably the difficult part is working with the $\arccos$ function. I have no idea how to deal with it.
Is it possible that there is an analytic solution to this integral? Even if someone found a numerical approximation, it would help me out.
 A: I’m not sure why you gave up on polar coordinates and why you say that the Jacobian is nasty – the Jacobian is well-known and simple, and the integrals all work out nicely.
We have $\mathrm dx\mathrm dy=r\mathrm dr\mathrm d\phi$. Denote the polar coordinates of the two points by $r_1,\phi_1$ and $r_2,\phi_2$. We need two separate cases according as the points are in the same or different octants.
For the same octant, the integral over the region $\phi_2\lt\phi_1\lt\frac\pi4$ is
\begin{eqnarray}
&&
\int_0^\frac\pi4\mathrm d\phi_1\int_0^{\phi_1}\mathrm d\phi_2\int_0^\frac1{\cos\phi_1}\mathrm r_1dr_1\int_0^\frac1{\cos\phi_2}\mathrm r_2dr_2\left(\phi_1-\phi_2\right)
\\
&=&
\frac14\int_0^\frac\pi4\mathrm d\phi_1\frac1{\cos^2\phi_1}\int_0^{\phi_1}\mathrm d\phi_2\frac{\phi_1-\phi_2}{\cos^2\phi_2}
\\
&=&
\frac14\int_0^\frac\pi4\mathrm d\phi_1\frac1{\cos^2\phi_1}\left[\left(\phi_1-\phi_2\right)\tan\phi_2-\log\cos\phi_2\right]_0^{\phi_1}
\\
&=&
-\frac14\int_0^\frac\pi4\mathrm d\phi_1\frac{\log\cos\phi_1}{\cos^2\phi_1}
\\
&=&
-\frac14\left[\tan\phi_1(1+\log\cos\phi_1)-\phi_1\right]_0^\frac\pi4
\\
&=&
\frac\pi{16}+\frac18\log2-\frac14\;.
\end{eqnarray}
For different octants, the integral over the region $\phi_2\lt\frac\pi4\lt\phi_1$ is
\begin{eqnarray}
&&
\int_\frac\pi4^\frac\pi2\mathrm d\phi_1\int_0^\frac\pi4\mathrm d\phi_2\int_0^\frac1{\cos\left(\frac\pi2-\phi_1\right)}\mathrm r_1dr_1\int_0^\frac1{\cos\phi_2}\mathrm r_2dr_2\left(\phi_1-\phi_2\right)
\\
&=&
\frac14\int_\frac\pi4^\frac\pi2\mathrm d\phi_1\frac1{\cos^2\left(\frac\pi2-\phi_1\right)}\int_0^\frac\pi4\mathrm d\phi_2\frac{\phi_1-\phi_2}{\cos^2\phi_2}
\\
&=&
\frac14\int_\frac\pi4^\frac\pi2\mathrm d\phi_1\frac1{\cos^2\left(\frac\pi2-\phi_1\right)}\left[\left(\phi_1-\phi_2\right)\tan\phi_2-\log\cos\phi_2\right]_0^\frac\pi4
\\
&=&
\frac14\int_\frac\pi4^\frac\pi2\mathrm d\phi_1\frac{\phi_1-\frac\pi4+\frac12\log2}{\cos^2\left(\frac\pi2-\phi_1\right)}
\\
&=&
\frac14\int_0^\frac\pi4\mathrm d\phi\frac{\frac\pi4+\frac12\log2-\phi}{\cos^2\phi}
\\
&=&
\frac14\left[\left(\frac\pi4+\frac12\log2\right)\tan\phi-\left(\phi\tan\phi+\log\cos\phi\right)\right]_0^\frac\pi4
\\
&=&
\frac14\log2\;.
\end{eqnarray}
There are $4$ symmetric copies of the first region and $2$ of the second, for a total of
$$
\boxed{\frac\pi4+\log2-1\approx0.4785}\;.
$$
This is not too different from the value $\frac\pi6\approx0.5236$ if the points are picked from the first quadrant of the unit disk.
Here’s Java code that checks the result by simulation.
A: Attempt to compute the 4-d integral is a wrong approach. 
For this sort of problem, the first thing to do is minimize the numbers of variable you need to work with as much as possible. In this case, you can work with the angles of the two vectors directly.

Let 


*

*$\theta_1 = \tan^{-1}\frac{b}{a}$ and $\theta_2 = \tan^{-1}\frac{d}{c}$.

*$\theta_m = \min(\theta_1,\theta_1)$ and $\theta_M = \max(\theta_1,\theta_2)$.

*$f(\theta)$ be the common CDF for $\theta_1, \theta_2$, i.e.
$$f(\theta) 
= \verb/CDF/_{\theta_k}(\theta) \stackrel{def}{=} \verb/P/[ \theta_k \le \theta ]\quad\text{ for }\quad k = 1, 2$$
The average angle you want is the expected value of $|\theta_1 - \theta_2| = \theta_M - \theta_m$,
To compute this, we need the CDF for $\theta_m$ and $\theta_M$.
$$\begin{align}
\verb/CDF/_{\theta_m}(\theta) &\stackrel{def}{=} \verb/P/[ \theta_m \le \theta ]
= \verb/P/[ \theta_1 \le \theta \lor \theta_2 \le \theta ]\\
&= 1 - \verb/P/[\theta_1 > \theta \land \theta_2 > \theta ]
 = 1 - \verb/P/[\theta_1 > \theta]\verb/P/[\theta_2 > \theta]\\
&= 1 - (1 - f(\theta))^2\\
\verb/CDF/_{\theta_M}(\theta) &\stackrel{def}{=} \verb/P/[ \theta_M \le \theta ]
= \verb/P/[ \theta_1 \le \theta \land \theta_2 \le \theta ]
= \verb/P/[\theta_1 \le \theta]\verb/P/[\theta_2 \le \theta ]\\
&= f(\theta)^2
\end{align}$$
This allows us to express the average angle as an integral
$$\begin{align}\mathcal{A} \stackrel{def}{=} \verb/E/\big[ |\theta_1 - \theta_2| \big]
& = \verb/E/[\theta_M - \theta_m]
 = \int_0^\frac{\pi}{2} \theta (\verb/CDF/'_{\theta_M} - \verb/CDF/'_{\theta_m}) d\theta\\
&= \int_0^{\frac{\pi}{2}} \theta [ f^2(\theta) - 1 + (f(\theta)-1)^2 ]' d\theta\\
&= 2 \int_0^{\frac{\pi}{2}} \theta [ f(\theta)(f(\theta)-1) ]' d\theta
\end{align}
$$
Integrate by part and notice $f(\theta) + f(\frac{\pi}{2} - \theta) = 1$, we find
$$\mathcal{A} = 2 \int_0^{\pi/2} f(\theta)(1 - f(\theta)) d\theta
= 4\int_0^{\pi/4} f(\theta)(1 - f(\theta)) d\theta$$
It is easy to see $\;f(\theta) = \frac12\tan\theta\;$ for $\;\theta \in [0,\frac{\pi}{4}]$. Change variable to $t = \tan\theta$ and integrate, the end result is:
$$\mathcal{A} = \int_0^1 \frac{t(2-t)}{1+t^2} dt = \frac{\pi}{4} + \log(2) - 1$$
