根据概率密度函数生成随机数的代码

我这里并不是要讲“伪随机”、“真随机”这样的问题，而是关于如何生成服从某个概率分布的随机数（或者说 sample）的问题。比如，你想要从一个服从正态分布的随机变量得到 100 个样本，那么肯定抽到接近其均值的样本的概率要大许多，从而导致抽到的样本很多是集中在那附近的。当然，要解决这个问题，我们通常都假设我们已经有了一个生成 0 到 1 之间均匀分布的随机数的工具，就好像 random.org 给我们的结果那样，事实上许多时候我们也并不太关心它们是真随机数还是伪随机数，看起来差不多就行了。

现在再回到我们的问题，看起来似乎是很简单的，按照概率分布的话，只要在概率密度大的地方多抽一些样本不就行了吗？可是具体要怎么做呢？要真动起手来，似乎有不是那么直观了。实际上，这个问题曾经也是困扰了我很久，最近又被人问起，那我们不妨在这里一起来总结一下。为了避免一下子就陷入抽象的公式推导，那就还是从一个简单的具体例子出发好了，假设我们要抽样的概率分布其概率密度函数为，并且被限制在区间上，如右上图所示。

好了，假设现在我们要抽 100 个服从这个分布的随机数，直观上来讲，抽出来的接近 3 的数字肯定要比接近 0 的数字要多。那究竟要怎样抽才能得到这样的结果呢？由于我们实际上是不能控制最原始的随机数生成过程的，我们只能得到一组均匀分布的随机数，而这组随机数的生成过程对于我们完全是透明的，所以，我们能做的只有把这组均匀分布的随机数做一些变换让他符合我们的需求。找到下手的点了，可是究竟要怎样变换呢？有一个变换相信大家都是很熟悉的，假设我们有一组之间的均匀分布的随机数，那么令的话，就是一组在之间均匀分布的随机数了，不难想象，等于某个数的概率就是等于的概率（“等于某个数的概率”这种说法对于连续型随机变量来说其实是不合适的，不过大概可以理解所表达的意思啦）。似乎有一种可以“逆转回去”的感觉了。

于是让我们来考虑更一般的变换。首先，我们知道的概率密度函数是，假设现在我们令，不妨先假定是严格单调递增的函数，这样我们可以求其逆函数（也是严格单调递增的）。现在来看变换后的随机变量会服从一个什么样的分布呢？

这里需要小心，因为这里都是连续型的随机变量，并不像离散型随机变量那样可以说成“等于某个值的概率”，因此我们需要转换为概率分布函数来处理，也就是求一个积分啦：

那么的概率分布函数为。很显然小于或等于某个特定的值这件事情是等价于这件事情的。换句话说，等于。于是，的概率分布函数就可以得到了：

再求导我们就能得到的概率密度函数：

这样一来，我们就得到了对于一个随机变量进行一个映射之后得到的随即变量的分布，那么，回到我们刚才的问题，我们想让这个结果分布就是我们所求的，然后再反推得即可：

经过简单的化简就可以得到，亦即。也就是说，把得到的随机数带入到到函数中所得到的结果，就是符合我们预期要求的随机数啦！让我们来验证一下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16

#!/usr/bin/python   import numpy as np import matplotlib.pyplot as plot   N = 10000 X0 = np.random.rand(N) X1 = 3*X0 Y = np.power(9*X1, 1.0/3)   t = np.arange(0.0, 3.0, 0.01) y = t*t/9  plot.plot(t, y, 'r-', linewidth=1) plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)plot.show()

这就没错啦，目的达成啦！让我们来总结一下。问题是这样的，我们有一个服从均匀分布的随机变量，它的概率密度函数为一个常数，如果是上的分布，那么常数就直接等于 1 了。现在我们要得到一个随机变量使其概率密度函数为，做法就是构造出一个函数满足（在这里加上了绝对值符号，这是因为如果不是递增而是递减的话，推导的过程中有一处就需要反过来）

反推过来就是，对目标的概率密度函数求一个积分（其实就是得到它的概率分布函数 CDF ，如果一开始就拿到的是 CDF 当然更好），然后求其反函数就可以得到需要的变换了。实际上，这种方法有一个听起来稍微专业一点的名字：Inverse Transform Sampling Method 。不过，虽然看起来很简单，但是实际操作起来却比较困难，因为对于许多函数来说，求逆是比较困难的，求积分就更困难了，如果写不出解析解，不得已只能用数值方法来逼近的话，计算效率就很让人担心了。可事实上也是如此，就连我们最常见的一维标准正态分布，也很难用这样的方法来抽样，因为它的概率密度函数

的不定积分没有一个解析形式。这可真是一点也不好玩，费了这么大劲，结果好像什么都干不了。看来这个看似简单的问题似乎还是比较复杂的，不过也不要灰心，至少对于高斯分布来说，我们还有一个叫做 Box Muller 的方法可以专门来做这个事情。因为高斯分布比较奇怪，虽然一维的时候概率分布函数无法写出解析式，但是二维的情况却可以通过一些技巧得出一个解析式来。

首先我们来考虑一个二维的且两个维度相互独立的高斯分布，它的概率密度函数为

这个分布是关于原点对称的，如果考虑使用极坐标（其中）的话，我们有，这样的变换。这样，概率密度函数是写成：

注意到在给定的情况下其概率密度是不依赖于的，也就是说对于来说是一个均匀分布，这和我们所了解的标准正态分布也是符合的：在一个圆上的点的概率是相等的。确定了的分布，让我们再来看，用类似于前面的方法：

根据前面得出的结论，我现在得到了的概率分布函数，是不是只要求一下逆就可以得到一个了？亦即。

现在只要把这一些线索串起来，假设我们有两个相互独立的平均分布在上的随机变量和，那么就可以得到了，而就得到了（实际上，由于和实际上是相同的分布，所以通常直接写为）。再把极坐标换回笛卡尔坐标：

这样我们就能得到一个二维的正态分布的抽样了。可以直观地验证一下，二维不太好画，就画成 heatmap 了，看着比较热的区域就是概率比较大的，程序如下：

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19

#!/usr/bin/python   import numpy as np import matplotlib.pyplot as plot   N = 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()

画出来的图像这个样子：

不太好看，但是大概的形状是可以看出来的。其实有了二维的高斯分布，再注意到两个维度在我们这里是相互独立的，那么直接取其中任意一个维度，就是一个一维高斯分布了。如下：

如果即服从标准正态分布的话，则有，也就是说，有了标准正态分布，其他所有的正态分布的抽样也都可以完成了。这下总算有点心满意足了。不过别急，还有最后一个问题：多元高斯分布。一般最常用不就是二元吗？二元不是我们一开始就推出来了吗？推出来了确实没错，不过我们考虑的是最简单的情形，当然同样可以通过这样的方式来处理每一个维度，不过高维的情形还有一个需要考虑的就是各个维度之间的相关性——我们之前处理的都是两个维度相互独立的情况。对于一般的多维正态分布，如果各个维度之间是相互独立的，就对应于协方差矩阵是一个对角阵，但是如果在非对角线的地方存在非零元素的话，就说明对应的两个维度之间存在相关性。

这个问题还是比较好解决的，高斯分布有这样的性质：类似于一维的情况，对于多维正态分布，那么新的随机变量将会满足

所以，对于一个给定的高斯分布来说，只要先生成一个对应维度的标准正态分布，然后令即可，其中是对进行 Cholesky Decomposition 的结果，即。

结束之前让我们来看看 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 要好一点啦！

然后，到这里为止，我们算是把高斯分布弄清楚了，不过这只是给一个介绍性的东西，里面的数学推导也并不严格，而 Box Muller 也并不是最高效的高斯采样的算法，不过，就算我们不打算再深入讨论高斯采样，采样这个问题本身也还有许多不尽人意的地方，我们推导出来的结论可以说只能用于一小部分简单的分布，连高斯分布都要通过 trick 来解决，另一些本身连概率密度函数都写不出来或者有各种奇怪数学特性的分布就更难处理了。所以本文的标题里也说了，这是上篇，如果什么时候有机会抽出时间来写下篇的话，我将会介绍一些更加通用和强大的方法，诸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比较感兴趣，可以先自行 Google 一下解馋！

根据概率密度函数生成随机数的代码相关推荐

计算机的随机数函数,使用rand()函数生成随机数
rand()函数生成随机数.计算机无法生成真正的随机数,它只产生所谓的伪随机数. 生成伪随机数,需要stdlib.h头文件. 参考以下代码: #include #include int main() ...
C++ 中生成随机数的代码
[知识点] 在C++中,可直接调用rand()函数生成随机数.不过,在调用rand()函数之前,需要先使用srand(time(0))函数设置随机数种子. 如果没有使用srand(time(0))函数 ...
vb.net产生随机数Random代码实例
Private Sub Button1_Click(sender As Object, e As EventArgs) Handles Button1.Click ' 随机数Dim a As Rand ...
arcgis字段计算器使用rnd函数生成随机数
arcgis字段计算器使用rnd函数生成随机数 VB Script里 rnd应该是个函数,取0到1中的随机数,乘10再减去5就能保证结果在基准值左右.截图中为了要整数的结果,又取了个整,如果小数结果也 ...
Java不重复的随机数获取_java获取 1--N 的不重复随机数程序代码
文章分享一篇关于java获取 1--N 的不重复随机数程序代码,有感兴趣的同学可以参考一下. 代码如下复制代码 import java.util.ArrayList; import java.ut ...
C语言rand函数生成随机数详解和示例
文章目录 1.生成随机数 2.生成一定范围随机数 3.获取视频教程 4.版权声明在C/C++程序开发中,会经常用到随机数这个功能,例如编写游戏类(纸牌)的程序时就需要用到随机数. 1.生成随机数在 ...
PHP内置函数生成随机数的方法汇总
PHP内部生成随机数的方法相比其他方法简单,不需要额外配置,是生成随机数的首选方案. 1 rand函数 rand() 函数可以不加任何参数,就可以生成随机整数.如果要设置随机数范围,可以在函数中设置 ...
Math函数生成随机数用法
1.Math.random() 生成0~0.9999的随机数 [0,1) 代码: @Test public void random(){double random = Math.random();// ...
php中生成随机数种子的函数有哪些,PHP内置函数生成随机数的方法汇总
PHP内部生成随机数的方法相比其他方法简单,不需要额外配置,是生成随机数的首选方案. 1 rand函数 rand() 函数可以不加任何参数,就可以生成随机整数.如果要设置随机数范围,可以在函数中设置 ...

根据概率密度函数生成随机数的代码

根据概率密度函数生成随机数的代码相关推荐

最新文章

热门文章