N维空间两随机向量的问题

[guofei9987]

https://github.com/MathAndAlgo/DiscussionBoard/blob/master/%E9%9A%8F%E6%9C%BA%E5%90%91%E9%87%8F%E7%9A%84%E4%B8%80%E4%BA%9B%E9%97%AE%E9%A2%98.docx

[shouldsee]

摘自docx

N=3时，近似均匀分布，有没有更深层的原因的？
提一个民科一点的想法：量子力学有MWI诠释，大概意思是每次量子行为都分裂成两个宇宙。因为宏观来说，维度太高，因此近似垂直，也就不相互影响。维度较低，例如只有一个光子时，两个宇宙会互相干涉，形成条纹。
提问：N=3的均匀分布，是否与现实世界的3维有关？

[shouldsee]

N=2 的情况比较好解释
V = <v1,v2> = cos(w1)cos(w2) + sin(w1)sin(w2) = cos(w1-w2).
w1-w2 \approx U(-pi, pi)

但是要推广的话还是用n-sphere上的均匀分布的余弦会比较好算。

[shouldsee]

若用旋转法考虑三维空间上球面的面积元。

$$\begin{align}
x^2 + y^2 + z^2 &= 1\\
z &= 0\\
dS &= 2 \pi ( y ) ( \sqrt{ 1 + (dy/dx)^2 } dx ) \\
dy/dx &= - \sqrt{x/(1-x^2)}\\
dS &= 2\pi (\sqrt{1 - x^2})( \sqrt{1/(1-x^2)}  ) dx \\
dS &= 2\pi dx \,\,\text{(independent of x)}\\
\end{align}$$

由于P(S)dS = P(x)dx,且P(S)是均匀分布，所以P(x) = 1 * dS/dx。因此只需考察高维球面的dS/dx便可解出向量夹角cosine的分布。高维球面面积和体积的讨论参见wiki-n-sphere，且与Gamma函数有关。所以n=3的情况可以从这个角度再作理解。

[Joe-C-Ding]

空间里两个向量 X Y 的余弦可以写成内积形式

如果向量各分量间独立的话，这个东西就只跟 Z = x1 * y1 和 W = (x1)^2 的分布有关了。
如果假设 x1 (y1) 服从 [-M, M] 上的均匀分布的话，刚算了一下它们的分布是

它们的均值是 EZ=0，EW=M^2/3。因此维数较高的时候，用中心极制定理近似和分布，比较容易得出两向量大概率垂直的结果。

但三维的时候正好是均匀分布不太好验证，因为这两个分布都没可加性，感觉只能硬算。
(我感觉最后的结果应该和 M 的取值没关系，所以取 M=1 大概就行，分布的形式最简单)。

[Joe-C-Ding]

若用旋转法考虑三维空间上球面的面积元。

$$\begin{align}
x^2 + y^2 + z^2 &= 1\\
z &= 0\\
dS &= 2 \pi ( y ) ( \sqrt{ 1 + (dy/dx)^2 } dx ) \\
dy/dx &= - \sqrt{x/(1-x^2)}\\
dS &= 2\pi (\sqrt{1 - x^2})( \sqrt{1/(1-x^2)}  ) dx \\
dS &= 2\pi dx \,\,\text{(independent of x)}\\
\end{align}$$

由于P(S)dS = P(x)dx,且P(S)是均匀分布，所以P(x) = 1 * dS/dx。因此只需考察高维球面的dS/dx便可解出向量夹角cosine的分布。高维球面面积和体积的讨论参见wiki-n-sphere，且与Gamma函数有关。所以n=3的情况可以从这个角度再作理解。

我感觉 @guofei9987 的本意是空间里任意两个向量吧，不一定是球面上的。不过球面上是不是有类似的结论倒是也可以讨论看看。

[shouldsee]

@joe 高维球面是考虑随机向量之间的余弦的最自然的情景，因为向量长度恒定了，而且可以得出一些简化的答案，所以我直接考虑了这个问题。因此，投影到超球面上可以更系统地处理余弦的分布的计。

不过@guofei的确没有正式描述随机向量。 你所给出的分布是高维立方体的体积上的分布，而立方体对应的是L1-norm。L1-norm下的内积我并没有计算过。

[guofei9987]

补充一下，一开始，我把“随机向量”定义为“多维正方形”中任取一点，但后来发现把“随机变量”定义为“多维球内任意一点”或者定义为“多维球面上任意一点”，从图像上看，结论似乎不变。

也就是说，从以下三种集合中做随机模拟实验，看起来结论一致：

1. “多维正方形”
2. “多维球体”
3. “多维球面”

下面是代码

多维正方形的：

# %%方形随机
import numpy as np
import matplotlib.pyplot as plt

n_dim = 3

cos_all = []
for i in range(10000):
    vec_1 = np.random.rand(n_dim) - 0.5
    vec_2 = np.random.rand(n_dim) - 0.5
    cos_tmp = (sum(vec_1 * vec_2)) / np.sqrt(sum(vec_1 ** 2) * sum(vec_2 ** 2))
    cos_all.append(cos_tmp)

plt.hist(cos_all)
plt.show()

多维球面不太容易模拟，似乎应该再开一个 issue 讨论一下“如何在任意曲面上随机选取一点，以达成曲面上的随机分布”
好在球面上的点和球体内的点可以有一种一致性的映射，也就是说，在讨论夹角余弦时，可以把样本总体定位球体内的点集，其结论应当等价于“球面上的点集”结果。
代码：

import numpy as np
import matplotlib.pyplot as plt

n = 3
cos_all = []
for i in range(10000):
    vec_1 = 2 * np.random.rand(n) - 1
    vec_2 = 2 * np.random.rand(n) - 1
    print(sum(vec_1 ** 2))
    if sum(vec_1 ** 2) <= 1 and sum(vec_2 ** 2) <= 1:
        cos_tmp = (sum(vec_1 * vec_2)) / np.sqrt(sum(vec_1 ** 2) * sum(vec_2 ** 2))
        cos_all.append(cos_tmp)

print(len(cos_all))
plt.hist(cos_all)
plt.show()

奇怪之处是，球体上的结论，与正方体上的结论似乎一致：n_dim=3时，余弦值近似均匀分布。

@shouldsee @Joe-C-Ding

[guofei9987]

用公式做一下恒等变换：
假设两个向量分别是 $X=(x_1,...,x_n), Y=(y_1,...,y_n)$
可以不妨假设 $Y=(1,0,...,0)$,
那么$<x,y>=x_1/|X|$
做一个平方，得到 $x_1^2/(x_1^2+x_2^2+...+x_n^2)$,其中，$x_i$是相互独立的均匀分布

这样，问题转化为，上面这个式子的分布律如何计算的问题

作为验证：

import numpy as np
import matplotlib.pyplot as plt

n = 4
cos_all = []
for i in range(10000):
    vec_1 = 2 * np.random.rand(n) - 1
    if sum(vec_1 ** 2) <= 1:
        cos_tmp = vec_1[0] / np.sqrt(sum(vec_1 ** 2))
        cos_all.append(cos_tmp)

print(len(cos_all))
plt.hist(cos_all)
plt.show()

[guofei9987]

继续推导，
$$\pm\sqrt{x_1^2/(x_1^2+x_2^2+...+x_n^2)}
=\pm \sqrt{1/(1+(x_2^2+...+x_n^2)/x_1^2)}
$$

发现，
$$(x_2^2+...+x_n^2)/x_1^2$$ 这一部分可以用卷积+拉普拉斯变幻来求解概率密度函数，需要一定计算量，为了规避这个，我们再做一个变换，来让问题简化：

1. 我们不妨假设$x_i$不再服从均匀分布，而是正态分布，结论大概不会变化。
2. 目标分布是对称的，所以可以忽略负半部分

研究 $\sqrt{1/(1+(n-1)/t^2(n-1))}$
当n=3时，
$\sqrt{1/(1+2/t^2(2))}$，记为K

$P(K\leq k)=P(t^2(2)\leq 2/(1/k^2-1))=2P(0\leq t(2)\leq \sqrt{2/(1/k^2-1)})$
$=2P(t(2)\leq \sqrt{2/(1/k^2-1)})-1$

而，
$P(t(2)\leq \sqrt{2/(1/k^2-1)})=\int_0^{\sqrt{2/(1/k^2-1)}} f_T(k) dk$
求导数，（记$b(k)=\sqrt{2/(1/k^2-1)}$）
得到 $f_K(k)=b'(k)f_T(b(k))$
把T分布的概率密度函数代入上式，化简得到一个常数。

如此证明了 目标分布确实是一个均匀分布

[guofei9987]

证明是完事了，但还有问题：

1. 只是证明了球体和球面上的集合成立，但没有证明或者证否正方形上的情况
2. 没有回答我最初的问题，只有3维向量互相投影为均匀分布，这是否与现实的三维空间有关？

[shouldsee]

@guofei9987
随机曲面上的采样我之前在RandomUnitaryVectors.ipynb有稍做讨论。但的确最简单的采样办法也是利用球体到球面的映射。

正方体的我还要求解一下。PS:GH对latex太不友好了。
正方形缺少对称性，太难算了，我还是放弃吧。
不过通过我的模拟，得出的结论是三维正方体体积上的均匀分布上的随机向量之间的夹角cosine并不服从均匀分布。

import numpy as np
N = int(50E4)
### 正方体
x = np.random.random(size=(N,3)) -0.5
y = np.random.random(size=(N,3)) -0.5

### 球体
# x = np.random.normal(size=(N,3))
# y = np.random.normal(size=(N,3))

def l2__norm(x,axis=1):
    return np.sqrt(np.sum(x**2,axis=axis,keepdims=True))

cosine = np.sum(x*y,axis=1,keepdims=1) / l2__norm(x) / l2__norm(y)

import matplotlib.pyplot as plt
plt.hist(cosine,bins=40);
plt.title("N={N}".format(**locals()))

[guofei9987]

@shouldsee 感谢，之前用的默认bins做随机模拟，没能看出来cos值的变化，
另外补充一个倒数第二步的模拟

from scipy import stats
n_dim=3
rv=stats.t(n_dim-1).rvs(1000000)
a=np.sqrt(1/(1+(n_dim-1)/(t_rvs**2)))
plt.hist(a,bins=40)
plt.show()

ps，怎么插图呀

[shouldsee]

你的公式建议用CodeCog渲染一下不然我实在没法读。

插图用的是 Ctrl+C -> 选定GH输入框 -> Ctrl+V -> 等待上传完成。
我稍微改了一下。

import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
n_dim=3
t_rvs=stats.t(n_dim-1).rvs(1000000)
a=np.sqrt(1/(1+(n_dim-1)/(t_rvs**2)))
plt.hist(a,bins=40)
plt.title("$a=\sqrt{1/(1+2/t^2)}$")
plt.show()

[guofei9987]

整理了一下
另外，随机曲面上的采样的代码不太好读，我这里装不了那个包，能不能描述一下？@shouldsee

[shouldsee]

随机曲面那个本子是很久以前的了，但是采样的想法也就只是从多维高斯分布采样然后用L2度量正规化投射到球面。本子更多地在讨论球面上的非均匀分布的采样以及写一些可视化。你有兴趣的话我可以抽空新开一本写，但是我最后直接用了vmf分布做超球面上的建模。

我没有看懂最后一步那里求导的过程，希望能够给一电对定积分求导的教程。

总的来说这种硬算读起来还是有一种诡异的暴力美感....而t分布和chi分布中的gamma函数本身可能也来自于超球面的几何。毕竟超球可以看作高维空间的体积源所以还是能体现一些本质问题的。有时间不妨把n=N时余弦的现实表达写出来做个比较。

[guofei9987]

定积分求导那个是定理，用定积分的定义也不难证明。
下面是那个定理，没写证明。
http://www.guofei.site/2017/10/26/multivariablecalculus.html#%E5%AE%9A%E7%A7%AF%E5%88%86%E6%B1%82%E5%AF%BC

随机曲面哪个，有空就写写呗，感觉挺有趣的

[shouldsee]

我查了一下你给的那个是Lebniz定理 https://en.wikipedia.org/wiki/Leibniz_integral_rule
高维球上的随机采样有空的话我会写一下rejection sampling或者reparametrisation，但要等到我哪个周末闲下来了。如果不是球面的话，我也不知道怎么采样的。

[guofei9987]

@shouldsee 针对一般曲面，我已经想出一个数学上可行，但编程上不太可取的方案。这周末可以整理完。
再开个 issue 吧