Given two positive real numbers $a$ and $b$, one can define their arithmetic-geometric mean as $M(a,b):=\ell$, where $$ (\ell,\ell)=\lim_{k\to\infty}(a_k,b_k)\quad\text{and}\quad\begin{cases}a_{k+1}=\frac{a_k+b_k}{2} \\ b_{k+1}=\sqrt{a_kb_k} \\ a_0=a,\ b_0=b.\end{cases} $$ On the other hand, one can define the elliptic integral $$ I(a,b):=\int_0^{\frac{\pi}{2}}\frac{d\theta}{\sqrt{a^2\cos^2\theta+b^2\sin^2\theta}}. $$ According to this interesting blog post by Paramanand, the famous formula $$\boxed{M(a,b)I(a,b)=\frac{\pi}{2}}$$ was discovered by Gauss, thanks to some (actually very hidden) numerical evidence. Using a clever change of variables, it can be shown that $I(a,b)=I(\frac{a+b}{2},\sqrt{ab})$, from which the identity follows. Actually, Gauss' first proof was a rather tedious computation with power series.
Is there a simple, natural proof of the identity, explaining how one could guess it (without numerical evidence) in the first place?