如何生成随机数（上）

快三个月没有写日志了，大概是我开始认真写 blog 来第一次，也是因为发生了一些预料之外的事情，中断了许久，到后来又一直非常非常忙，不过我终于又爬上来冒个泡了，表明我还活着。

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

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

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

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

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

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

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

#!/usr/bin/pythonimport numpy as np
import matplotlib.pyplot as plotN = 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/9plot.plot(t, y, 'r-', linewidth=1)
plot.hist(Y, bins=50, normed=1, facecolor='green', alpha=0.75)
plot.show()

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

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

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

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

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

#!/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()

画出来的图像这个样子：

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

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

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

<img src="http://blog.pluskid.org/latexrender/pictures/bf449b4ed053041e7429be818334130c.png" _xhe_src="http://blog.pluskid.org/latexrender/pictures/bf449b4ed053041e7429be818334130c.png" title="" \displaystyle"="" alt="" align="absmiddle" style="border: 0px; margin-left: auto; margin-right: auto;">

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

结束之前让我们来看看 matlab 画个 3D 图来改善一下心情：

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 一下解馋！

末了，还得废话两句，虽然是冒了一个泡，但是也是冒得很艰难，发了这篇日志也并不代表 blog 会恢复往日的发文频率，实际上这篇 blog 也差不多凑了半个多月的时间才折腾完，文章开头的“三个月”我想是不是也该更新为“四个月”了。不过我想，如果有时间和有趣的 idea 的话，我还是会写的，毕竟这也是我的乐趣之一呢！

如何生成随机数（上）相关推荐

java随机数语句_Java语言程序设计（七）Math类生成随机数及if语句
Java有几种类型的选择语句,单向if语句,双向if语句,嵌套if语句,switch语句和条件表达式. 1.单向if语句 if(radius>=0){ area = radius*radius* ...
vrf名称_如何使用VRF(可验证随机函数)在以太坊上生成随机数
Chainlink 如何解决以太坊"随机数问题" 随机数和区块链一直很难达到"一致"(译者注:区块链要求确定性,而随机数正相反).到目前为止,区块链上还没有可验 ...
汇总|C++常见知识点总结，涉及文本输出、排序、生成随机数、异常处理、关联容器、printf重定向、sprintf用法、cout重定向
文章目录一将程序运行结果输出到txt文本文件中二排序算法三生成随机数四异常处理六关于GitHub上zip与tar.gz的区别七容器中查找最大值所在的位置八 C++中关联容器的 ...
Java学习之生成随机数
1.Java中的方法random()可用于生成随机数,称为伪随机数生成器,它返回一个大于等于0.0.小于1.0的数(double类型),即0.0<=X<1.0 .之所以产生的数称为伪随机数 ...
单片机生成随机数的方法总结
去年冬天在帮学校附近一家密室逃脱店做一些电子机关,其中一个打地鼠项目需要用到单片机产生随机数,用于实现随机让几个地鼠"钻"出来.一开始想法很单纯,不就是随机函数么,之前C语言课上就 ...
java 随机数生成实现_Java中生成随机数的实现方法总结
搜索热词在实际开发工作中经常需要用到随机数.如有些系统中创建用户后会给用户一个随机的初始化密码.这个密码由于是随机的,为此往往只有用户自己知道.他们获取了这个随机密码之后,需要马上去系统中更改.这就 ...
python 编程一日一练-「每日一练」巧用python生成随机数
原标题:「每日一练」巧用python生成随机数随机数在我们的生产和生活中有很多的应用场景,比如说登录验证的随机数字等等,那么你知道在Python中怎么生成随机数吗? 往下看,就是这么简单! 题目 p ...
python产生5个随机数_Python和numpy生成随机数
http://blog.csdn.net/pipisorry/article/details/39086463 随机数种子要每次产生随机数相同就要设置种子,相同种子数的Random对象,相同次数生成 ...
Numpy 生成随机数和乱序
Numpy 生成随机数和乱序 Numpy官网:http://www.numpy.org/ 一.生成随机数 1. numpy.random.rand(d0, d1, -, dn) 生成在 [0, 1) ...
从最基础的讲起如何做到均匀的生成随机数
题目描述:470. 用 Rand7() 实现 Rand10() 因为是第一次接触到这样的题目,毫无思绪,对官方题解也是"不知道为什么要这么做".看过一些题解之后才逐渐明白,现在让我 ...

如何生成随机数（上）

如何生成随机数（上）相关推荐

最新文章

热门文章