问题
单位正方形内两个随机点的距离期望是多少？
求解因由
何昆师兄博士毕业之际，小龙师兄发朋友圈称赞其深刻的数学思维时，提到这个问题，引起了我的兴趣。我的数学天分一般，小时候没出现机缘参加数学竞赛的训练，但对欣赏数学形式的优美这件事情从来都有着强烈的共鸣，特别是通过复杂的推导过程得到简洁的解析解之后。
上学读书时候，我对随机和概率这样的观念和工具总是不太认同，这大概也源于我认知天性中对殊相比共相更贴近真相的一种倾向。不过，在阅读了《复杂》和一些量子力学的科普读物后，我渐渐感受到随机和概率是当下和未来最具确定性的工具，虽然局限性非常明显，即需要大量的观察和经验，且多适用于集体行为的预测。不过话又说回来，这本也是当下科学方法论的局限。语言文字和符号的局限性数学语言无法躲避，如今科学方法论的局限性，数学也难以躲避，可是要说发现形式上的美感这件事情，随机计算领域应该蕴藏丰富，聊以慰藉受困于苦痛中的人们。
解后感受
思考过程
单位长度上两个随机点的距离期望是多少？
先看单位长度上一个点的随机。仔细考量，发现，能从两个角度来思考：分布律和概率密度函数。
这两个角度能够通过几何方法产生关联和替换，但是彼此之间有着鲜明的鸿沟，即不是一样的东西。人们常讲“点构成线，线构成面，面构成体”，这是分布律的眼光，不是概率密度函数的眼光。通过分布律的眼光，才产生了极限、无穷多的势这样的观念。以概率密度函数的眼光来看，线和点完全就是不一样的两种东西，任意一段线可以由单位的线段计数组成，但绝对不是由点组成。再举例来讲，数字1既用来表示单位长度的线段，也用来表示单位大小的正方形，也用来表示单位大小的立方体。集合和映射的出现，丰富了代数手段，却容易让人忘记数字所代表的含义。数字再抽象，脱离了所表示的东西就显得毫无意义，推导本身也变得荒唐。
角度一：分布律
将单位长度转化为间距相等的\(n\)个点，这样随机两个点就容易找样本空间和组合数了。随机抛点的问题转变为了从\(n\)个数中有放回的取两次的问题。样本空间是
$$n^2$$
转换后，距离也容易度量，相邻两点间的距离是
$$\frac1 {n-1}$$
有了这两个基础，两个随机点的距离取值范围就清楚了，即
$$0, \frac1 {n-1}, \frac2 {n-1}, …, \frac{n-2} {n-1}, 1$$
令\(L\)为随机事件——取任意两点——所得的距离，那么
$$P(L=\frac{k}{n-1})=\frac{2(n-k)}{n^2}$$
其中
$$1 \le k \le n-1$$
对于距离为\(0\)的情况，其概率为
$$P(L=0)=\frac{n}{n^2}=\frac{1}{n}$$
根据分布律的规定，必须有
$$\sum_{k=0}^{n-1}{P(L=\frac{k}{n-1})}=1$$
验证
$$\begin{align*}
\sum_{k=0}^{n-1}{P(L=\frac{k}{n-1})}
&=\frac{1}{n^2}[n+2\sum_{k=1}^{n-1}(n-k)] \\
&=\frac{1}{n^2}[n+2\cdot\frac{n(n-1)}{2}] \\
&=1
\end{align*}$$
数学期望是对分布律上随机变量均值的描述，形式化为
$$\begin{align*}
\sum_{k=0}^{n-1}{L \cdot P(L=\frac{k}{n-1})}&=\sum_{k=1}^{n-1}\frac{2k(n-k)}{(n-1)n^2} \\
&=\frac{2}{(n-1)n^2}\sum_{k=1}^{n-1}{k(n-k)} \\
&=\frac{n+1}{3n}
\end{align*}$$
以”点密成线“的观点来看，当\(n\)趋于无穷时候，所得的期望值，即是单位长度上随机两点距离的期望值，为
$$\lim_{n \to \infty}{\frac{n+1}{3n}}=\frac{1}{3}$$
角度二：概率密度函数
连续型的概率只能讲同等规模的描述对象事件发生的概率，所有低于此等规模的描述对象的概率都是零。比如单位长度上，只能问点落在某个长度区域内的概率是多少，这样的概率值会是一个非零的数，因为这样的计算的描述对象是线段，样本空间是线段长度，不是点。倘若问点落在1/2处的概率，那一定是零，因为点相对于线段是不同的东西，甚至都很难说是不同量级的东西，定义这样事件发生的概率为零也没什么坏处。连续型的概率函数的事件通常都写成区间的形式，比如单位长度上一个随机点落在某个区域
$$[x_1, x_2]$$
这个事件表示为
$$x_1\le X\le x_2$$
基于前面的讨论，这个事件也等于
$$x_1< X < x_2$$
略去区域左边的边界（一般是零或者负的无穷大），概率函数为
$$P(X\le x)=x$$
其中
$$X \in [0,1], x \in [0, 1], X\le x \subseteq [0, 1]$$
概率密度函数表征概率函数的变化率，形式上为概率函数的导数，因而单位长度上随机均匀分布的概率密度函数是
$$f(x)=1$$
验证
$$\int_{0}^{1}f(x)\mathrm{d}x=1$$
令\(L\)为随机事件——任意两点随机抛落在单位长度区间内——所得的距离，那么随机事件可以表示为
$$L=|X_1-X_2|\le l$$
其中
$$0\le l \le 1$$
且\(X_1\)与\(X_2\)服从前述的均匀分布，即满足
$$0\le X_1 \le 1, 0\le X_2 \le 1$$
事件所在的样本空间为单位正方形，事件发生区域为二维图像上的一块面积，如图所示灰色区域。
令\(y=X_2, x=X_1\)用积分的方法求灰色区域的面积为
$$\begin{align*}
&P(L\le l) \\
&=\int_0^l{(x+l)\mathrm{d}x}+\int_{l}^{1-l}{(2l)\mathrm{d}x}+\int_{1-l}^{1}(1-x+l)\mathrm{d}x\\
&=2l-l^2
\end{align*}$$
验证当\(l=1\)时，有
$$P(L\le 1)=1$$
对分布函数求导数，得到距离的概率密度函数为
$$f(l)=2-2l$$
求数学期望，得
$$\mathrm{E}(f(l))=\int_0^1{l\cdot (2-2l)\mathrm{d}l}=\frac{1}{3}$$
不论用分布律还是概率密度函数推导都能得出一样的期望值，即1/3. 同理，原问题也能应用同样的两个思路求解，不过求解难度却明显增加，两种方法增加难度之处不大一样。
解题过程
角度一：分布律
假设单位正方形内均匀分布了
$$n^2$$
个点，沿笛卡尔横纵坐标轴方向的相邻两点的距离均为
$$\frac{1}{n-1}$$
每个点可以通过唯一的坐标定位，即有矩阵
$$\frac{1}{n-1}\cdot \left[ \begin{array}{c}
(0, 0) & (0, 1) & … & (0, n-1) \\
(1, 0) & (1, 1) & … & (1, n-1) \\
… & … & … & … \\
(n-1,0) & (n-1,1) & … & (n-1, n-1)
\end{array} \right ]
$$
令随机变量\(L\)表示为随机从\(n^2\)个点中有放回的抽取两次所得两点之间的欧几里得距离，于是开始考察距离取值的可能性。假设抽取的两个点，坐标分别为\((x_1, y_1)\)与\((x_2,y_2)\)，根据一维问题的经验，知道横坐标和纵坐标差值的取值可能性的数量各有\(n\)种，即
$$\vec{d}=\frac{1}{n-1}\cdot \left[ \begin{array}{c}
0 \\
1 \\
… \\
n-1
\end{array} \right ]
$$
两点欧式距离取值如一下矩阵所示，即
$$\begin{align*}
\mathbf{D}&=\frac{1}{n-1} \\
&\cdot \left[ \begin{array}{c}
0 & 1 & … & n-1 \\
1 & \sqrt{2} & … & \sqrt{1^2+(n-1)^2}\\
… & … & … & … \\
n-1 & \sqrt{(n-1)^2+1^2} & … & (n-1)\sqrt{2}
\end{array} \right ]
\end{align*}$$
矩阵的每个格子中距离\(d\)的计算方法，即
$$\begin{align*}
L&=d[(x_1, y_1), (x_2, y_2)] \\
&=d(|x_1-x_2|, |y_1-y_2|) \\
&=d(\frac{k_1}{n-1}, \frac{k_2}{n-1}) \\
&=\frac{\sqrt{k_1^2+k_2^2}}{n-1}
\end{align*}$$
其中
$$0\le k_1, k_2 \le n-1$$
观察距离取值的矩阵\(\mathbf{D}\)，发现其具有对称性，只看右上三角矩阵即可，也就是只看
$$0\le k_1 \le k_2 \le n-1$$
的那些项。至此，事件发生的样本空间清楚了，取值有
$$\frac{n(n+1)}{2}$$
种可能性。虽然样本空间只有平方的量级，但抽样方法却是在四次方的样本空间下发生的，因而还得通过有放回的抽取四个数值（\(x_1, y_1, x_2, y_2\)）的事件或者有放回的抽取两个数值（\(k_1, k_2\)）的事件来模拟。不管是哪种，都需要分四类情况计算事件发生的概率。
- 情况一：\(L=0\)
也即
$$x_1=x_2 \wedge y_1=y_2$$
或是
$$k_1=k_2=0$$
这种情况表明两次抽到了同一个点，样本空间是\(n^4\)，事件发生次数是\(n^2\)，因而事件发生概率为
$$P(L=0)=\frac{n^2}{n^4}=\frac{1}{n^2}$$
- 情况二：\(L=\frac{k}{n-1}\)
也即
$$x_1=x_2 \wedge y_1 \ne y_2 \vee x_1 \ne x_2 \wedge y_1 = y_2 $$
或是
$$k=k_1 \ne 0 \wedge k_2 = 0 \vee k=k2 \ne 0 \vee k_1=0$$
这种情况表明两次抽到的两个点来自同一横坐标或者同一纵坐标上的不同点。差值为
$$L=\frac{k}{n-1}$$
其中
$$1\le k \le n-1$$
的事件发生次数可以这样计算：假设两点横坐标相同，先选横坐标，有
$$n$$
种选法，满足差值的纵坐标有
$$2(n-k)$$
种选法。同理若两点纵坐标相同，根据对称性，得到相同的结论。综合得到事件发生的概率为
$$P(L=\frac{k}{n-1})=\frac{4n(n-k)}{n^4}$$
其中
$$1\le k \le n-1$$
- 情况三：\(L=\frac{\sqrt{2}k}{n-1}\)
也即
$$|x_1-x_2|=|y_1-y_2| \ne 0$$
或者
$$k=k_1=k_2 \ne 0$$
这种情况表明两个点的横坐标差值和纵坐标差值相等，给定差值，那么横纵坐标各只需要选一次。对于差值为
$$L=\frac{\sqrt{2}k}{n-1}$$
其中
$$1\le k \le n-1$$
的事件发生次数可以这样计算：选横坐标，有
$$2(n-k)$$
种选法；选纵坐标，也有这样多种选法。依据乘法原理，综合得到事件发生的概率为
$$P(L=\frac{\sqrt{2}k}{n-1})=\frac{4(n-k)^2}{n^4}$$
其中
$$1\le k \le n-1$$
- 情况四：\(L=\frac{\sqrt{k_1^2+k_2^2}}{n-1}\)
也即
$$|x_1-x_2|\ne |y_1-y_2| \ne 0$$
或者
$$k_1 \ne k_2 \ne 0$$
这种情况对应到矩阵\(\mathbf{D}\)上，就是除去首行、首列、对角线外的其余部分，因为横纵坐标差在距离的计算上具有对称性，因而最终算得的概率系数为8. 计算方法类似情况二和三。
$$P(L=\frac{\sqrt{k_1^2+k_2^2}}{n-1})=\frac{8(n-k_1)(n-k_2)}{n^4}$$
其中
$$1\le k_1 < k_2 \le n-1$$
综合四种情况，得到分布律为
$$\left\{\begin{array}{llr}
P(L=0)=\frac{n^2}{n^4}=\frac{1}{n^2} \\
P(L=\frac{k}{n-1})=\frac{4n(n-k)}{n^4} & 1\le k \le n-1 \\
P(L=\frac{\sqrt{2}k}{n-1})=\frac{4(n-k)^2}{n^4} & 1\le k \le n-1\\
P(L=\frac{\sqrt{k_1^2+k_2^2}}{n-1})=\frac{8(n-k_1)(n-k_2)}{n^4} & 1\le k_1 < k_2 \le n-1
\end{array}\right.
$$
验证所求分布律是否正确
$$\begin{align*}
&P(L=0)+\sum_{k=1}^{n-1}{P(L=\frac{k}{n-1})} \\
&+\sum_{k=1}^{n-1}{P(L=\frac{\sqrt{2}k}{n-1})} \\
&+\sum_{k_1=1}^{n-2}{\sum_{k_2=k_1+1}^{n-1}{P(L=\frac{\sqrt{k_1^2+k_2^2}}{n-1})}} \\
&=\frac{1}{n^2}+\frac{4}{n^3}\sum_{k=1}^{n-1}{(n-k)} \\
&+\frac{4}{n^4}\sum_{k=1}^{n-1}{(n-k)^2} \\
&+\frac{8}{n^4}\sum_{k_1=1}^{n-2}{\sum_{k_2=k_1+1}^{n-1}(n-k_1)(n-k_2)} \\
&=\frac{1}{n^2}+\frac{4}{n^3}\cdot \frac{n(n-1)}{2} \\
&+\frac{4}{n^4}\cdot \frac{n(n-1)(2n-1)}{6}\\
&+\frac{8}{n^4}\cdot \sum_{k_1=1}^{n-2}{(n-k_1)\sum_{k_2=k_1+1}^{n-1}(n-k_2)} \\
&=\frac{1}{n^2}+\frac{4}{n^3}\cdot \frac{n(n-1)}{2} \\
&+\frac{4}{n^4}\cdot \frac{n(n-1)(2n-1)}{6} \\
&+\frac{8}{n^4}\cdot \sum_{k_1=1}^{n-2}{(n-k_1)^2(n-k_1-1)} \\
&=\frac{2n-1}{n^2}+\frac{4}{n^4}\cdot \frac{n(n-1)(2n-1)}{6} \\
&+\frac{4}{n^4}\cdot [(\frac{n(n-1)}{2})^2-\frac{n(n-1)(2n-1)}{6}] \\
&=\frac{2n-1}{n^2}+\frac{(n-1)^2}{n^2} \\
&=1
\end{align*}
$$
该分布律的数学期望为
$$\begin{align*}
\mathrm{E_n}(L) &= 0\times \frac{1}{n^2} \\
&+ \sum_{k=1}^{n-1}\frac{k}{n-1}\frac{4n(n-k)}{n^4} \\
&+ \sum_{k=1}^{n-1}\frac{\sqrt{2}k}{n-1}\frac{4(n-k)^2}{n^4} \\
&+ \sum_{k_1=1}^{n-2}{\sum_{k_2=k_1+1}^{n-1}\frac{\sqrt{k_1^2+k_2^2}}{n-1}\frac{8(n-k_1)(n-k_2)}{n^4}}
\end{align*}$$
所求的单位正方形中随机两点的数学期望为
$$\begin{align*}
\mathrm{E}(L)
&=\lim_{n\to \infty}\mathrm{E_n}(L) \\
&=\lim_{n\to \infty}\sum_{k_1=1}^{n-2}{\sum_{k_2=k_1+1}^{n-1}\frac{\sqrt{k_1^2+k_2^2}}{n-1}\frac{8(n-k_1)(n-k_2)}{n^4}}
\end{align*}$$
以上等式之所以成立在于当\(n\)趋于无穷大时，中间两项的分子都只是\(n\)的三次幂和四次幂，极限值都为零。对于第四项的和，我曾经一度寻求不等式放缩，但发现无论怎样放缩缝隙都非常大，并且不可能找到一个\(k_1\)和\(k_2\)构成的线性组合去替换\(\sqrt{k_1^2+k_2^2}\)。其实当把极限的数学形式写出来后，应该容易联想到转换为积分计算能将问题简化，毕竟在\(n\)趋向于无穷大时，\(n\)的最高幂次的系数决定了最终的结果，于是
$$\begin{align*}
&\lim_{n\to \infty}\mathrm{E_n}(L) \\
&=\lim_{n\to \infty}\frac{8}{(n-1)n^4}\sum_{k_1=1}^{n-2}{\sum_{k_2=k_1+1}^{n-1}\sqrt{k_1^2+k_2^2}(n-k_1)(n-k_2)} \\
&=\lim_{n\to \infty}\frac{8}{(n-1)n^4}\int_0^n\int_x^n\sqrt{x^2+y^2}(n-x)(n-y)dydx
\end{align*}$$
根据对称性，有
$$\begin{align*}
&\int_0^n\int_x^n\sqrt{x^2+y^2}(n-x)(n-y)dydx \\
&=\int_0^n\int_0^x\sqrt{x^2+y^2}(n-x)(n-y)dydx
\end{align*}$$
观察积分区域，发现其为90度的等腰三角形，转换到极坐标系更容易计算，于是
$$x=\rho\cos\theta$$
$$y=\rho\sin\theta$$
积分区域从
$$0\le x \le n$$
$$0\le y \le x$$
变为
$$0\le \theta \le \frac{\pi}{4}$$
$$0\le \rho \le n\sec\theta$$
有
$$\begin{align*}
&\int_0^n\int_0^x\sqrt{x^2+y^2}(n-x)(n-y)\mathrm{d}y\mathrm{d}x \\
&=\int_0^{\frac{\pi}{4}}\int_0^{n\sec\theta}\rho(n-\rho\cos\theta)(n-\rho\sin\theta)\rho\mathrm{d}\rho\mathrm{d}\theta
\end{align*}$$
积分的内容展开后一共四项，根据积分的线性可叠加性，可以拆开逐项计算
- 第一项：
$$
\begin{align*}
&n^2\int_0^{\frac{\pi}{4}}\int_0^{n\sec\theta}\rho^2\mathrm{d}\rho\mathrm{d}\theta \\
&=\frac{n^5}{3}\int_0^{\frac{\pi}{4}}\sec^3\theta\mathrm{d}\theta \\
&=\frac{n^5}{3}\cdot \frac{1}{2}[\frac{\sin\theta}{\cos^2\theta}+\ln|\tan\theta+\sec\theta|]_0^{\frac{\pi}{4}} \\
&=\frac{n^5}{6}(\sqrt{2}+\ln(1+\sqrt{2}))
\end{align*}$$
- 第二项：
$$
\begin{align*}
n\int_0^{\frac{\pi}{4}}\cos\theta\int_0^{n\sec\theta}\rho^3\mathrm{d}\rho\mathrm{d}\theta
&=\frac{n^5}{4}\int_0^{\frac{\pi}{4}}\sec^3\theta\mathrm{d}\theta \\
&=\frac{n^5}{8}(\sqrt{2}+\ln(1+\sqrt{2}))
\end{align*}$$
- 第三项：
$$
\begin{align*}
n\int_0^{\frac{\pi}{4}}\sin\theta\int_0^{n\sec\theta}\rho^3\mathrm{d}\rho\mathrm{d}\theta
&=\frac{n^5}{4}\int_0^{\frac{\pi}{4}}\sin\theta\sec^4\theta\mathrm{d}\theta \\
&=\frac{n^5}{4}\int_0^{\frac{\pi}{4}}\cos^{-4}\theta\mathrm{d}\cos\theta \\
&=\frac{n^5}{12}(2\sqrt{2}-1)
\end{align*}$$
- 第四项：
$$
\begin{align*}
\int_0^{\frac{\pi}{4}}\cos\theta\sin\theta\int_0^{n\sec\theta}\rho^4\mathrm{d}\rho\mathrm{d}\theta
&=\frac{n^5}{5}\int_0^{\frac{\pi}{4}}\sin\theta\sec^4\theta\mathrm{d}\theta \\
&=\frac{n^5}{15}(2\sqrt{2}-1)
\end{align*}$$
将计算所得带回原式，合并化简得到
$$\begin{align*}
&\lim_{n\to \infty}\mathrm{E_n}(L) \\
&=\lim_{n\to \infty}\frac{8}{(n-1)n^4}\cdot \frac{[5\ln(1+\sqrt{2})+\sqrt{2}+2]n^5}{120} \\
&= \frac{5\ln(1+\sqrt{2})+\sqrt{2}+2}{15} \\
&\approx {0.5214}
\end{align*}$$
角度二：概率密度函数
类比一维的解法，令\(L\)为随机事件——任意两点随机抛落在单位正方形区间内——所得的距离，那么随机事件可以表示为
$$L=\sqrt{(X_1-X_2)^2+(Y_1-Y_2)^2}\le l$$
其中
$$0\le l \le 1$$
且\((X_1, Y_1)\)与\((X_2, Y_2)\)服从均匀分布
$$P(X\le x, Y\le y)=xy$$
其中
$$X, Y \in [0, 1]$$
类比前面的思路，下一步要求概率分布
$$P(L\le l)=\iiiint\limits_{\sqrt{(x_1-x_2)^2+(y_1-y_2)^2}\le l}f(x_1,x_2,y_1,y_2)\mathrm{d}x_1\mathrm{d}x_2\mathrm{d}y_1\mathrm{d}y_2$$
