Apex angle of a triangle as a random variable I am not an expert in Probability Theory and I apologize if I make some mistakes in "translating" the following problem into the language of random variables. Any help also to improve the formulation of the problem is very much welcome.
Problem. Given two angles $\alpha, \beta \in (0,2\pi)$, consider the isosceles triangle with the apex angle $\theta = \alpha - \beta$ (or $\pi -(\alpha -\beta)$ if $\alpha - \beta > \pi$) and leg lengths 1. By well known elementary geometry, (twice) the area (without sign) of the triangle is given by $\vert \sin \theta \vert$. 
Now consider a probability measure $\mu$ and suppose $\alpha, \beta$ are independent random variables w.r.t. to this probability measure. 

Q. Which is the probability measure $\mu$ which maximizes the mean area of the triangle, i.e. 
  $$
\int_0^{2\pi}\int_0^{2\pi} | \sin(\alpha - \beta)| d\mu(\alpha) d \mu(\beta)?
$$ 

Thanks to @JeanMarie's comments below, we can get rid of an angle (by sticking one of the sides of the triangle to the x-axis) so the question above can be equivalently formulated as follows: 

Q. Which is the probability measure $\mu$ which maximizes the mean area of the triangle, i.e. 
  $$
\int_0^{2\pi} | \sin(\alpha)| d\mu(\alpha)?
$$ 

For instance, if $\mu$ is the normalized Lebesgue measure on $(0,2\pi)$ the area of the triangle is 
$$
\frac{1}{4\pi}\int_0^{2\pi} \int_0^{2\pi} \vert \sin (\alpha-\beta) \vert d\alpha d \beta = \frac{1}{2}. 
$$
Is this the maximum value?
 A: Here is a proof that the uniform distribution over $[0, 2\pi]$ (i.e. the normalized Lebesgue measure) maximizes the integral, based on variational argument. The idea in this proof is based on my other answer for the similar question.
On the space $\mathcal{M}$ of signed Borel measures on $[0, 2\pi]$, define $\langle \cdot, \cdot \rangle$ by
$$ \forall \mu, \nu \in \mathcal{M}, \quad  \langle \mu, \nu \rangle = \int_{[0,2\pi]^2} \lvert \sin(x-y)\rvert \, \mu(\mathrm{d}x)\nu(\mathrm{d}y). $$
Obviously $\langle \cdot, \cdot \rangle$ is a symmetric bilinear form on $\mathcal{M}$. Using this, OP's question is rephrased as the problem of finding maximizers of the functional $I(\mu) = \langle \mu, \mu \rangle$ over the set of all probability measures on $[0, 2\pi]$. In this regard, we have the following characterization:

Proposition. Write $\mathcal{M}_1$ for the set of all $\mu \in \mathcal{M}$ satisfying $\mu([0,2\pi]) = 1$. Then, for $\mu \in \mathcal{M}_1$, the followings are equivalent:

*

*$\mu$ maximizes the functional $I$.

*$s \mapsto \langle \mu, \delta_s \rangle$ is constant for $s \in [0, 2\pi]$.


Assuming this proposition, we see that the normalized Lebesgue measure $\mu = \frac{1}{2\pi} \operatorname{Leb}|_{[0,2\pi]}$ satisfies the second condition, and so, it maximizes $I$ not only among probability measures but on all of $\mathcal{M}_1$. (On the other hand, I have not proved the uniqueness of the maximizer. I guess that the maximizer is unique, based on the idea that the map $s \mapsto \langle \mu, \delta_s \rangle$ behaves like Fourier transform, and so, it should have inverse transform.)
As for the proof, we resolve an easier part first. This is a typical application of the variational argument.
Proof of $(1)\Rightarrow(2)$. Assume that $\mu \in \mathcal{M}_1$ maximizes the functional $I$. Then for any $s, t \in [0, 2\pi]$, the map $ \epsilon \mapsto I(\mu + \epsilon \delta_s - \epsilon \delta_t) $ attains the maximum at $\epsilon = 0$. Differentiating with respect to $\epsilon$ and plugging $\epsilon = 0$ proves the desired implication. $\square$
So we turn to proving the opposite direction.
Proof of $(2)\Rightarrow(1)$. We begin by citing the following Fourier series.
$$ \lvert \sin x \rvert = \frac{2}{\pi} - \frac{4}{\pi} \sum_{n=1}^{\infty} \frac{\cos (2nx)}{4n^2 - 1}. $$
For instance, see this answer, for the proof. To make use of this, write $\hat{\mu}(\xi) = \int_{[0,2\pi]} e^{i \xi x} \, \mu(\mathrm{d}x)$ for the Fourier transform of $\mu$. Now, since this series converges uniformly, we can plug this formula to $\langle \mu, \nu \rangle$ and interchange the order of the integration and summation. This yields
$$ \langle \mu, \nu \rangle
= \frac{2}{\pi} \hat{\mu}(0)\overline{\hat{\nu}(0)} - \frac{4}{\pi} \sum_{n=1}^{\infty} \frac{1}{4n^2 - 1}\operatorname{Re}\left( \hat{\mu}(2n)\overline{\hat{\nu}(2n)} \right). $$
From this, we immediately observe that
$$ \text{if} \quad \mu([0,2\pi]) = 0, \qquad \text{then} \quad \langle \mu, \mu \rangle
= -\frac{4}{\pi} \sum_{n=1}^{\infty} \frac{1}{4n^2 - 1}\left|\hat{\mu}(2n)\right|^2 \leq 0. $$
Now we are ready for the final blow. Let $\mu \in \mathcal{M}_1$ be such that $\langle \mu, \delta_s \rangle = \int_{[0,2\pi]} \lvert \sin(x-s)\rvert \, \mu(\mathrm{d}x)$ does not depend on $s \in [0,2\pi]$. If $c$ denotes this constant value, then $ \langle \mu, \nu \rangle = c $ holds for any $\nu \in \mathcal{M}_1$. So,
$$ \forall \nu \in \mathcal{M}_1, \qquad
I(\nu)
= I(\mu) + 2\underbrace{\langle \mu, \nu-\mu \rangle}_{c-c = 0} + \underbrace{I(\nu - \mu)}_{\leq 0} \leq I(\mu). $$
Therefore $\mu$ maximizes $I$ over $\mathcal{M}_1$ as desired. $\square$
A: @Sangchul Lee: The proof is very nice and natural. In general, it is well known that the Lebesgue measure minimizes the energy  $\int\limits_{S^{d-1}}  \int\limits_{S^{d-1}} F (\langle x,y \rangle ) \, d\mu (x) d\mu (y)$ over probability measures  on the sphere $S^{d-1}$ iff the underlying function  $F:[-1,1] \rightarrow \mathbb R$ is positive definite on $S^{d-1}$, up to an additive constant, which according to a famous characterization due to Schoenberg is equivalent to saying that the expansion of $F$ in Gegenbauer polynomials has positive coefficients. 
In this case, we are on $S^1$, so Gegenbauer polynomials become Chebyshev polynomials (of the first kind), and through the relation $T_n (\cos \theta) = \cos (n\theta)$ this can be identified with Fourier series. So it is quite natural to use them. Here we are maximizing, NOT minimizing,  hence the important property is that all coefficients in the Fourier expansion are negative (non-positive).
In fact the proof could have been made a bit shorter: after $\langle\mu, \nu \rangle$ is computed, it is obvious that our energy satisfies  $\langle \mu,\mu \rangle  \ge \frac{2}{\pi} |\hat{\mu} (0)|^2 = \frac{2}{\pi}$  (since $\hat{\mu} (0) = \mu (S^1) = 1$), and the equality is achieved for the Lebesgue measure. 
One should not expect uniqueness here: first, since the energy doesn't distinguish between $\alpha$ and $-\alpha$, so putting all the mass uniformly in $[0,\pi]$ would give the same energy. More subtly, since the expansion of $|\sin x|$ only has even frequencies of the cosine, one can perturb the Lebesgue measure by a cosine with an odd frequency, and the energy will not "see" it, i.e. give the same value as for the Lebesgue measure. 
Lastly, I think the statement is still true on the higher-dimensional spheres $S^{d-1}$. It gets just a bit more technical, one would need to show that the function $(1-t^2)^{1/2}$ has negative coefficients in its Gegenbauer polynomial expansion, which seems to work. (This function arises from taking $\langle x, y \rangle = \cos \theta$, so that $|\sin \theta| = \sqrt{1 - \langle x,y \rangle^2}$.) 
A: let us consider the second (reduced) case
$$\int_0^{2\pi}| \sin(\alpha)| \underbrace{f(\alpha) d\alpha}_{d\mu(\alpha)}\tag{1}$$
(I assume we are looking for absolutely continuous measures) This quantity (where $f$ is the looked-for density and d$\alpha$ is Lebesgue measure) is maximized by considering the case of equality in Cauchy-Schwarz identity, ($\int_a^b f(x)g(x)dx$ is maximal when $f$ and $g$ are proportional) which occurs if :
$$f(\alpha)=K| \sin(\alpha)| \ \text{with} \ K=\frac14 \tag{2}$$
(the value of normalizing constant $K$ comes from the fact that $\int_0^{2\pi}| \sin(\alpha)|d\alpha=4$).
Of course, (2) has to be understood "almost everywhere".
Now, what is this maximal value ?
$$\int_0^{2\pi}\frac14 |\sin(\alpha)|^2 d\alpha=2 \int_0^{\pi}\frac14 \sin(\alpha)^2 d\alpha=2 \frac14 \frac12 [\alpha - \frac12\sin(2 \alpha)]_0^{\pi}=\pi/4$$
which is very slightly bigger than the value we obtain with Lebesgue mesure i.e., $\frac{2}{\pi}$.
Instead of using a.c. measures, one could include measures with diracs by using Riemann-Stieljes integral : https://math.stackexchange.com/q/263474
A: Here is the explanation of the connection  between positive definite functions and energy minimization on the sphere,  that I promised  earlier, as well as an alternative solution to the  posted question, which extends to all dimensions. I'll try to keep preliminaries to a minimum. 
Let $F:[-1,1]\rightarrow \mathbb R$ be a continuous function. We say that the function $F$ is positive definite on the sphere $S^{d-1}$ if for any finite subset $Z = \{ z_1, \dots, z_N\} \subset S^{d-1}$ the matrix $\big[ F ( \langle z_i, z_j \rangle ) \big]_{i,j=1}^N$ is positive semidefinite. 
Let $\{P_n^d\}_{n=0}^\infty$ denote the sequence of Gegenbauer (ultraspherical)  polynomials, i.e.  polynomials on $[-1,1]$, orthogonal with respect to the weight $(1-t^2)^{\frac{d-3}2}$ (there are various normalizations and notations in the literature -- this is not the most common one, but it is simpler for our purposes).  Notice that the weight is the one that arises when integrating a zonal function on $S^{d-1}$, i.e.
\begin{equation}\label{e.int}
\int\limits_{S^{d-1}}   f ( \langle p, x \rangle ) d\sigma (x) =  \frac{\omega_{d-2}}{\omega_{d-1}} \int_{-1}^1 f (t ) \big(1-t^2\big)^{\frac{d-3}2}\, dt,
\end{equation}
where $\sigma$ is the normalized uniform surface measure on the sphere $S^{d-1}$,  $p \in S^{d-1}$ is an arbitrary pole, and $\omega_{d-1}$ is the Lebesgue surface measure of $S^{d-1}$.
When  $d=2$ these are Chebyshev polynomials (of the first kind), and for $d=3$ these are the Legendre polynomials. Since this is an orthogonal sequence, for any function $F$ as above, one can consider its Gegenbauer expansion  $\sum\limits_{n=0}^\infty  \widetilde{F}_n P_n^d (t)$. We take $P_0^d (t) \equiv 1$, and one can see from the first equation, that 
\begin{equation}
\widetilde{F}_0 =  \int\limits_{S^{d-1}}   F ( \langle p, x \rangle ) d\sigma (x) =   \int\limits_{S^{d-1}} \int\limits_{S^{d-1}}   F ( \langle  x, y  \rangle ) d\sigma (x) d\sigma (y), 
\end{equation}
where the second identity follows by rotational invariance (the expression in the middle is independent of $p$). The aforementioned theorem of Schoenberg states that $F$ is positive definite on the sphere if and only if $\widetilde{F}_n \ge 0$ for all $n \ge 0$.
In some special cases there is an easier way to check for positive definiteness. Assume $F$ is also analytic in $(-1,1)$ and  in its Taylor series has non-negative coefficients, then $F$ is positive definite on $S^{d-1}$ for any $d\ge 2$. This can be proved in a couple of ways: either by integration by parts and using Rodrigues formula, or alternatively as follows. Schur's theorem says that if matrices $A$ and $B$ are positive semidefinite, then their Hadamard (elementwise) product $A\circ B = [ a_{ij} b_{ij} ]$ is also positive semidefinite. Since the Gram matrix $\big[ \langle z_i, z_j \rangle  \big]_{i,j=1}^N$ is always positive semidefinite, by Schur's theorem so are its Hadamard powers $\big[ \langle z_i, z_j \rangle ^k \big]_{i,j=1}^N$. And since all coefficients in the Taylor series of $F$ are non-negative, then the matrix $\big[ F ( \langle z_i, z_j \rangle ) \big]_{i,j=1}^N$ is also positive-semidefinite. 
Now define the energy
\begin{equation}
I_F (\mu) =  \int\limits_{S^{d-1}} \int\limits_{S^{d-1}}   F ( \langle  x, y  \rangle ) d\mu (x) d\mu (y),
\end{equation}
and assume we want to minimize it over the class $\mathcal M$ of probability measures on the sphere $S^{d-1}$. 
Claim: if $F$ is positive definite on $S^{d-1}$, up to an additive constant, then the uniform distribution $\sigma$ is a minimizer of $I_F$ (this is sort of folklore). 
Here is one proof. Assume that  $F$ is positive definite on $S^{d-1}$, up to an additive constant, i.e. $\widetilde{F}_n \ge 0$ for all $n \ge 1$. Recall that $\widetilde{F}_0 = I_F (\sigma)$, as discussed above.  We shall use the ``addition formula" for spherical harmonics. Let $M_{n,d}$ be the dimension of the space of spherical harmonics on $S^{d-1}$ of degree $n$, and let $\{ H_{n,k} \}_{k=1}^{M_{n,d}}$ be an orthonormal basis of this space (with respect to $\sigma$). Then for any $x$, $y \in S^{d-1}$
\begin{equation}
\sum_{k=1}^{M_{n,d}}  H_{n,k} (x) H_{n,k} (y) = P_n^d ( \langle x, y \rangle),
\end{equation}
where $P_n^d$ is the $n^{th}$ Gegenbauer polynomial.  (In $d=2$, i.e. on the circle, this is equivalent to the formula $\cos (nx ) \cos (ny) + \sin (ny) \sin (ny) = \cos n (x-y)$, since spherical harmonics can be associated with the Fourier basis, and Chebyshev polynomials satisfy $T_n (\cos \theta) = \cos n \theta$.) We than have, using the addition formula,
\begin{align*}
I_F (\mu ) & = \widetilde{F}_0 + \sum_{n\ge1} \widetilde{F}_n \int\limits_{S^{d-1}} \int\limits_{S^{d-1}}   P_n^d ( \langle  x, y  \rangle ) d\mu (x) d\mu (y) \\
&= I_F (\sigma)  +  \sum_{n\ge1} \widetilde{F}_n \sum_{k=1}^{M_{n,d}} \bigg(  \int\limits_{S^{d-1}} H_{n,k} (x) d\mu (x) \bigg)^2 \ge I_F (\sigma),
\end{align*}
i.e. $\sigma$ is a minimizer. (Caution: I have used Fubini above, i.e. interchanged integration and summation, without justifying it -- but in fact there is an interesting self-improving property here: for positive definite $F$, its Gegenbauer series converges absolutely, this could be proved either using Mercer's theorem from spectral theory or directly using  properties of Gegenbauer polynomials.)
This prove is the most natural, but for completeness I include another, completely elementary, proof of the Claim. Recall that positive definiteness means that for all $Z = \{ z_1, \dots, z_N\} \subset S^{d-1}$  and all $c_1, \dots, c_N \in \mathbb R$ we have 
\begin{equation*}
I_F \big( \sum_{i} c_i \delta_{z_i}  \big) =   \sum_{i,j} F ( \langle z_i, z_j \rangle ) c_i c_j  \ge 0,
\end{equation*}
(it is easy to see that the expressions above are equal). Since linear combinations of Dirac delta masses are weak-$*$ dense in the space of Borel measures, this condition is equivalent to saying that $I_F (\nu) \ge 0$ for all signed Borel measures. Then observe the following, for any $\mu \in \mathcal M$
\begin{align*}
0  &\le  I_F ( \mu - \sigma) = I_F (\mu) + I_F (\sigma) - 2  \int\limits_{S^{d-1}} \int\limits_{S^{d-1}}   F ( \langle  x, y  \rangle ) d\sigma (x) d\mu (y)\\
& = I_F (\mu) + I_F (\sigma) - 2  \int\limits_{S^{d-1}} I_F (\sigma)  d\mu (y) =  I_F (\mu) - I_F (\sigma) ,
 \end{align*}
hence $\sigma$ is a minimizer. In the last line I used the fact that the inner integral in $\sigma$ is independent of $y$ and equals $I_F (\sigma)$. 
Finally, returning to our problem, our function is $| \sin \angle(x,y) | = \sqrt{1-\langle x,y \rangle^2}$, i.e. the underlying function $F(t) = (1-t^2)^{1/2}$. We can write the Taylor series 
\begin{equation}
(1-t^2)^{1/2}  = \sum_{n=0}^\infty  (-1)^n {{1/2} \choose n} t^{2n},
\end{equation}
and it is easy to see that all the coefficients (except $n=0$) are non-positive. Since we are maximizing (instead of minimizing), it follows from the discussion above that $\sigma$ is a maximizer in all dimensions $d\ge 2$. 
