python概率分布函数 python概率分布函数拟合( 九 )

\displaystyle \begin{aligned} P(rR) = \int_0^{2\pi}\int_0^R\frac{1}{2\pi}e^{\frac{r^2}{2}}rdrd\theta \ = \int_0^Re^{-\frac{r^2}{2}}rdr \ = 1-e^{-\frac{R^2}{2}} \end{aligned}
根据前面得出的结论，我现在得到了 r 的概率分布函数，是不是只要求一下逆就可以得到一个 \phi(\cdot) 了？亦即 \phi(t) = \sqrt{-2\log (1-t)}。
现在只要把这一些线索串起来，假设我们有两个相互独立的平均分布在 [0,1] 上的随机变量 T_1 和 T_2 ，那么 2\pi T_1 就可以得到 \theta 了，而 \phi(T_2) = \sqrt{-2\log(1-T_2)} 就得到 r 了（实际上，由于 T_2 和 1-T_2 实际上是相同的分布，所以通常直接写为 \sqrt{-2\log T_2}）。再把极坐标换回笛卡尔坐标：
\displaystyle \begin{aligned} x = r\cos\theta= \sqrt{-2\log T_2}\cdot \cos(2\pi T_1) \ y = r\sin\theta = \sqrt{-2\log T_2}\cdot \sin(2\pi T_1) \end{aligned}
这样我们就能得到一个二维的正态分布的抽样了。可以直观地验证一下，二维不太好画，就画成 heatmap 了，看着比较热的区域就是概率比较大的，程序如下：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19
#!/usr/bin/pythonimport numpy as np import matplotlib.pyplot as plotN = 50000 T1 = np.random.rand(N) T2 = np.random.rand(N)r = np.sqrt(-2*np.log(T2)) theta = 2*np.pi*T1 X = r*np.cos(theta) Y = r*np.sin(theta)heatmap, xedges, yedges = np.histogram2d(X, Y, bins=80) extent = [xedges[0], xedges[-1], yedges[0], yedges[-1]]plot.imshow(heatmap, extent=extent) plot.show()
画出来的图像这个样子：
不太好看，但是大概的形状是可以看出来的。其实有了二维的高斯分布，再注意到两个维度在我们这里是相互独立的，那么直接取其中任意一个维度，就是一个一维高斯分布了。如下：
如果 X\sim N(0,1) 即服从标准正态分布的话，则有 \sigma X+\mu \sim N(\mu, \sigma^2) ，也就是说，有了标准正态分布，其他所有的正态分布的抽样也都可以完成了。这下总算有点心满意足了。不过别急，还有最后一个问题：多元高斯分布。一般最常用不就是二元吗？二元不是我们一开始就推出来了吗？推出来了确实没错，不过我们考虑的是最简单的情形，当然同样可以通过 \sigma X+\mu 这样的方式来处理每一个维度，不过高维的情形还有一个需要考虑的就是各个维度之间的相关性——我们之前处理的都是两个维度相互独立的情况。对于一般的多维正态分布 X\sim N(\mathbf{\mu}, \Sigma) ，如果各个维度之间是相互独立的，就对应于协方差矩阵 \Sigma 是一个对角阵，但是如果 \Sigma 在非对角线的地方存在非零元素的话，就说明对应的两个维度之间存在相关性。
这个问题还是比较好解决的，高斯分布有这样的性质：类似于一维的情况，对于多维正态分布 X\sim N(\mathbf{\mu}, \Sigma)，那么新的随机变量 X_1=\mathbf{\mu}_1 + LX 将会满足
\displaystyle X_1 \sim N(\mathbf{\mu}_1+L\mu, L\Sigma L^T)
所以，对于一个给定的高斯分布 N(\mathbf{\mu}, \Sigma) 来说，只要先生成一个对应维度的标准正态分布 X\sim N(0, I) ，然后令 X_1 = \mu+LX 即可，其中 L 是对 \Sigma 进行 Cholesky Decomposition 的结果，即 \Sigma = LL^T。
结束之前让我们来看看 matlab 画个 3D 图来改善一下心情：
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
N = 50000; T1 = rand(1, N); T2 = rand(1, N);r = sqrt(-2*log(T2)); theta = 2*pi*T1; X =[r.*cos(theta); r.*sin(theta)]; mu = [1; 2]; Sigma = [5 2; 2 1]; L = chol(Sigma);X1 = repmat(mu,1, N) + L*X;nbin = 30;hist3(X1', [nbin nbin]); set(gcf, 'renderer', 'opengl'); set(get(gca,'child'), 'FaceColor', 'interp', 'CDataMode', 'auto');[z c] = hist3(X1', [nbin nbin]); [x y] =meshgrid(c{1}, c{2});figure; surfc(x,y,-z);
下面两幅图，哪幅好看一些（注意坐标比例不一样，所以看不出形状和旋转了）？似乎都不太好看，不过感觉还是比前面的 heatmap 要好一点啦！