by pluskid, on 2010-04-28, in 
快三个月没有写日志了，大概是我开始认真写 blog 来第一次，也是因为发生了一些预料之外的事情，中断了许久，到后来又一直非常非常忙，不过我终于又爬上来冒个泡了，表明我还活着。 🙂
第二点要澄清的是，我这里并不是要讲“伪随机”、“真随机”这样的问题，而是关于如何生成服从某个概率分布的随机数（或者说 sample）的问题。比如，你想要从一个服从正态分布的随机变量得到 100 个样本，那么肯定抽到接近其均值的样本的概率要大许多，从而导致抽到的样本很多是集中在那附近的。当然，要解决这个问题，我们通常都假设我们已经有了一个生成 0 到 1 之间均匀分布的随机数的工具，就好像 random.org 给我们的结果那样，事实上许多时候我们也并不太关心它们是真随机数还是伪随机数，看起来差不多就行了。 :p
现在再回到我们的问题，看起来似乎是很简单的，按照概率分布的话，只要在概率密度大的地方多抽一些样本不就行了吗？可是具体要怎么做呢？要真动起手来，似乎有不是那么直观了。实际上，这个问题曾经也是困扰了我很久，最近又被人问起，那我们不妨在这里一起来总结一下。为了避免一下子就陷入抽象的公式推导，那就还是从一个简单的具体例子出发好了，假设我们要抽样的概率分布其概率密度函数为 $p(x) = \frac{1}{9}x^2$ ，并且被限制在区间 $[0, 3]$ 上，如右上图所示。
好了，假设现在我们要抽 100 个服从这个分布的随机数，直观上来讲，抽出来的接近 3 的数字肯定要比接近 0 的数字要多。那究竟要怎样抽才能得到这样的结果呢？由于我们实际上是不能控制最原始的随机数生成过程的，我们只能得到一组均匀分布的随机数，而这组随机数的生成过程对于我们完全是透明的，所以，我们能做的只有把这组均匀分布的随机数做一些变换让他符合我们的需求。找到下手的点了，可是究竟要怎样变换呢？有一个变换相信大家都是很熟悉的，假设我们有一组 $[0,1]$ 之间的均匀分布的随机数 $X_0$ ，那么令 $X_1=3X_0$ 的话，$X_1$ 就是一组在 $[0,3]$ 之间均匀分布的随机数了，不难想象，$X_1$ 等于某个数 $x^*$ 的概率就是 $X_0$ 等于 $x^*/3$ 的概率（“等于某个数的概率”这种说法对于连续型随机变量来说其实是不合适的，不过大概可以理解所表达的意思啦）。似乎有一种可以“逆转回去”的感觉了。
于是让我们来考虑更一般的变换。首先，我们知道 $X_1$ 的概率密度函数是 $f(x) = 1/3, x\in[0,3]$ ，假设现在我们令 $Y = \phi (X_1)$ ，不妨先假定 $\phi(\cdot)$ 是严格单调递增的函数，这样我们可以求其逆函数 $\phi^{-1}(\cdot)$ （也是严格单调递增的）。现在来看变换后的随机变量 $Y$ 会服从一个什么样的分布呢？
这里需要小心，因为这里都是连续型的随机变量，并不像离散型随机变量那样可以说成“等于某个值的概率”，因此我们需要转换为概率分布函数来处理，也就是求一个积分啦：
\[ \displaystyle F(x) = P(X \leq x) = \int_{-\infty}^x f(t)dt \]那么 $X_1$ 的概率分布函数为 $F(x) = \frac{1}{3}x$ 。很显然 $Y$ 小于或等于某个特定的值 $y^*$ 这件事情是等价于 $X_1=\phi^{-1}(Y)\leq\phi^{-1}(y^*)$ 这件事情的。换句话说，$P(Y\leq y^*)$ 等于 $P(X_1 \leq \phi^{-1}(y^*))$ 。于是，$Y$ 的概率分布函数就可以得到了：
\[ \displaystyle G(y) = P(Y \leq y) = P(X_1 \leq \phi^{-1}(y)) = F(\phi^{-1}(y)) \]再求导我们就能得到 $Y$ 的概率密度函数：
\[ \displaystyle g(y) = \frac{dG(y)}{dy} = f(\phi^{-1}(y))\frac{d}{dy}\phi^{-1}(y) \]这样一来，我们就得到了对于一个随机变量进行一个映射 $\phi(\cdot)$ 之后得到的随即变量的分布，那么，回到我们刚才的问题，我们想让这个结果分布就是我们所求的，然后再反推得 $\phi(\cdot)$ 即可：
\[ \displaystyle \frac{1}{9}y^2 = g(y) = f(\phi^{-1}(y))\frac{d}{dy}\phi^{-1}(y) = \frac{1}{3}\frac{d}{dy}\phi^{-1}(y) \]经过简单的化简就可以得到 $\phi^{-1}(y) = \frac{1}{9} y^3$ ，亦即 $\phi(x) = (9x)^{1/3}$ 。也就是说，把得到的随机数 $X_1$ 带入到到函数 $\phi(\cdot)$ 中所得到的结果，就是符合我们预期要求的随机数啦！ 😀 让我们来验证一下：
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() |
这就没错啦，目的达成啦！让我们来总结一下。问题是这样的，我们有一个服从均匀分布的随机变量 $X$ ，它的概率密度函数为一个常数 $f(x)=C$ ，如果是 $[0,1]$ 上的分布，那么常数 $C$ 就直接等于 1 了。现在我们要得到一个随机变量 $Y$ 使其概率密度函数为 $g(y)$ ，做法就是构造出一个函数 $\phi(\cdot)$ 满足（在这里加上了绝对值符号，这是因为 $\phi(\cdot)$ 如果不是递增而是递减的话，推导的过程中有一处就需要反过来）
\[ \displaystyle g(y) = f(\phi^{-1}(y))\left|\frac{d}{dy}\phi^{-1}(y)\right| = C\left|\frac{d}{dy}\phi^{-1}(y)\right| \]反推过来就是，对目标 $y$ 的概率密度函数求一个积分（其实就是得到它的概率分布函数 CDF ，如果一开始就拿到的是 CDF 当然更好），然后求其反函数就可以得到需要的变换 $\phi(\cdot)$ 了。实际上，这种方法有一个听起来稍微专业一点的名字：Inverse Transform Sampling Method 。不过，虽然看起来很简单，但是实际操作起来却比较困难，因为对于许多函数来说，求逆是比较困难的，求积分就更困难了，如果写不出解析解，不得已只能用数值方法来逼近的话，计算效率就很让人担心了。可事实上也是如此，就连我们最常见的一维标准正态分布，也很难用这样的方法来抽样，因为它的概率密度函数
\[ \displaystyle g(y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{1}{2}y^2} \]的不定积分没有一个解析形式。这可真是一点也不好玩，费了这么大劲，结果好像什么都干不了。看来这个看似简单的问题似乎还是比较复杂的，不过也不要灰心，至少对于高斯分布来说，我们还有一个叫做 Box Muller 的方法可以专门来做这个事情。因为高斯分布比较奇怪，虽然一维的时候概率分布函数无法写出解析式，但是二维的情况却可以通过一些技巧得出一个解析式来。
首先我们来考虑一个二维的且两个维度相互独立的高斯分布，它的概率密度函数为
\[ \displaystyle f(x,y) = \frac{1}{\sqrt{2\pi}}e^{-\frac{x^2}{2}}\cdot\frac{1}{\sqrt{2\pi}}e^{-\frac{y^2}{2}} = \frac{1}{2\pi}e^{-\frac{x^2+y^2}{2}} \]这个分布是关于原点对称的，如果考虑使用极坐标 $(\theta,r)$ （其中 $\theta\in[0,2\pi), r\in[0,\infty)$ ）的话，我们有 $x = r\cos\theta$，$y=r\sin\theta$ 这样的变换。这样，概率密度函数是写成：
\[ \displaystyle f(\theta,r) = \frac{1}{2\pi}e^{-\frac{r^2}{2}}r \]注意到在给定 $r$ 的情况下其概率密度是不依赖于 $\theta$ 的，也就是说对于 $\theta$ 来说是一个均匀分布，这和我们所了解的标准正态分布也是符合的：在一个圆上的点的概率是相等的。确定了 $\theta$ 的分布，让我们再来看 $r$，用类似于前面的方法：
\[ \displaystyle \begin{aligned} P(r<R) &= \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/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() |
画出来的图像这个样子：
不太好看，但是大概的形状是可以看出来的。其实有了二维的高斯分布，再注意到两个维度在我们这里是相互独立的，那么直接取其中任意一个维度，就是一个一维高斯分布了。如下：
如果 $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 要好一点啦！
然后，到这里为止，我们算是把高斯分布弄清楚了，不过这只是给一个介绍性的东西，里面的数学推导也并不严格，而 Box Muller 也并不是最高效的高斯采样的算法，不过，就算我们不打算再深入讨论高斯采样，采样这个问题本身也还有许多不尽人意的地方，我们推导出来的结论可以说只能用于一小部分简单的分布，连高斯分布都要通过 trick 来解决，另一些本身连概率密度函数都写不出来或者有各种奇怪数学特性的分布就更难处理了。所以本文的标题里也说了，这是上篇，如果什么时候有机会抽出时间来写下篇的话，我将会介绍一些更加通用和强大的方法，诸如 Rejection Sampling 、Gibbs Sampling 以及 Markov Chain Monte Carlo (MCMC) 等方法。如果你比较感兴趣，可以先自行 Google 一下解馋！ 😀
末了，还得废话两句，虽然是冒了一个泡，但是也是冒得很艰难，发了这篇日志也并不代表 blog 会恢复往日的发文频率，实际上这篇 blog 也差不多凑了半个多月的时间才折腾完，文章开头的“三个月”我想是不是也该更新为“四个月”了。 :p 不过我想，如果有时间和有趣的 idea 的话，我还是会写的，毕竟这也是我的乐趣之一呢！ 😀
仰慕！
看到一半发现随机数原来有这么多的讲究，看到后面发现原来这么这么多讲究。
不愿意面对棘手的问题，yy 一种情况：如果事先确切知道要生成 n 个随机数 (n 比较大) ，和随机数的精度 (比如 0.01。或者说，最终的结果要有 m 个不同的数， m << n，可能有时满足不了) 的话，先按照值域和 p(x) 从小到大直接生成这 m 种共 n 个数。然后再使用常见的平均分布随机数给每个元素弄一个位置，达到乱序。
thx a lot
最近在精算课上刚好讲到随机数的问题，pluskid的这篇博文可以当做补充讲义了，不过多元分布的情况还没弄清楚，明日继续研究……
@quark
嗯，你说的方法相当于把连续型的随机变量用离散型随机变量去近似，就是说你用 m 个 bin 的 histogram 去替代原来的连续曲线了。如果数值取得好的话，应该也会在某些情况下可用的。
@layla
哦？“精算”是一门讲什么的课？以前倒是没有听说过。是金融经济这个方面的吗？
邮件提醒失效了吗？
随机数总是整数个，数值计算总有精度下界，我更喜欢觉得通用的办法。
昨天 itsuhane 说直接用普通的随机数 count bit 就是正态分布了。想一下就知道是正确的，MSTC 水楼里贴的 VB 代码就是一个比较弱的实现。prada 上次讲 Beautiful Code 某一章的时候，介绍了高效的 bit counting 算法
0 x 10 + 0 != 0 ?! -,-
啊？普通随机数做 bit counting 吗？这个倒是没有听说过。你讲一下？我看你那个 VB 知道是怎么做的了，不过不明白，为什么是正态分布？如果把每个 bit 看出独立的 0/1 的话，数 bit 感觉应该是二项分布才对啊？
ps: 这个 blog 里的邮件提醒似乎从来没有生效过，我不知道为什么，明明和另一边用的是同一个插件同一个 WP 版本还有同一个 Host ，我想不同的只有主题了。
@pluskid
二项分布的极限情况就是正态分布了。
邮件提醒有两种插件，一个是 Mail to commenter，这个是检查评论中是否有”@user”+空格，或者是”@user:”这样子，它是直接用冒号和空格切掉评论文字的，然后再去找”@”，所以你的 Reply 链接要改一改，加一个空格或冒号，并在 Mail to commenter 设置好就行了。
另外一种是针对多级嵌套评论的，那个用的是另外一个插件，要开启嵌套评论才能用，不过这个和你目前的 Reply 链接肯定冲突了，最里面一层的评论也没有 reply 按钮，这点不太好。
偶尔看到那么高质量的文章！
尤其是图片的解释非常直观阿
精算的确是金融related的一门课，主要是通过数学的方法对金融、保险产品的定价、投资、风险分析进行计算……所以概率论、数学分析、实变泛函之类的知识都会用到……
@quark
原来我用的那个 plugin 是根据 threaded comment 来回复的啊，我还以为一股脑之前的用户全部邮件提醒呢。这个 theme 好像不支持嵌套主题，那就算了，把你说的那个 plugin 拿过来 hack 了一下，现在应该有邮件提醒了吧？
@layla
好复杂的样子…… @_@
@pluskid
其实也还好，就像是保险公司承保寿险，它就需要根据生命表来估算风险，在根据赔付额和赔付概率给保险产品定价，以保证它们不会亏本~
我也刚刚开始学，真正的精算的确很麻烦，要不然精算师的薪水就不会那么高了……
@layla
唔，我还以为这些是可以由程序来计算的，也许情况太复杂了，或者“精算”这个词翻译得不够贴切了，猜测应该是更多分析和评估的成分吧。
@pluskid
恩，精算师的确需要掌握一定编程能力和相关软件使用能力的，比如matlab，mathematica，spss之类的~但是数学知识是基础，因为精算很多时候也要用到计量的东西……
不过我私人认为最重要的其实还是建模和分析能力……
@pluskid
嗯，现在好了。
不过 Mail to commenter 的英文邮件模板很山寨，邮件标题就有语法错误的样子，建议改一下~
@layla, kid怕金融了，哈哈
不会吧~数理基础和计算机好的人很适合学金融工程，我们也常被这类精英抢饭碗:)
好文章,LZ威武…
我现在有个类似的问题,给定一个随机数X,已知它满足离散分布P(X=xi)=pi,如何像P的分布有解析式那样O(1)地生成满足这样分布的X呢?或者根本不能?谢谢LZ!
你把问题想得太复杂了，离散型的随即变量的话，是最简单的了。比如，你有三个可能的取值，概率分别是 0.1、0.5、0.4 ，那么你只要用均匀分布的随机数生成器生成一个随机数，如果它在 [0, 0.1) 里就取 x1，在 [0.1, 0.6) 里就取 x2 ，在 [0.6, 1] 里就取 x3 ，这样就好了呀。
恩,和我想的一样,但n比较大的时候,得二分查找,所以不是O(1)的…
看来这种情况还是不如文中的解析情况来的和谐…
谢谢LZ
这个倒似乎是这样，不知道有没有标准的解决办法。不过我想，如果 n 不是非常非常大的话，二分其实是相当快的，如果 n 真是很大，那么也许可以考虑用连续的情况去近似离散的情况，用一个多项式函数来拟合出一个概率密度函数，得到结果之后直接 truncate 成整数。
恩,这个思路很好,因为可以先对xi按照pi排序为x’i,再用多项式拟合就是单调函数的情况了,由x’i到xi可以hash搞一搞…
这个问题原本是这样的,比如我有google一年接受的所有query及每个query的频数,我想在上面按原分布抽样10000个,除了按上面的思路重复10000次,LZ有否更好的思路?谢谢LZ…
我好像没有理解你的问题，不知道“重复 10000 次”是指什么意思？如果是构造多项式拟合的话，构造一次就够了，后面直接用构造好的函数来生成随机数就可以了呀，如果要生成 10000 个 sample 的话，只要生成 10000 个 [0,1] 上的随机数，再按照计算出来的映射求出对应的值，然后全部 truncate 一下就可以了呀。
恩，是你说的意思。。。。我现在就是二分10000次做的。。。算了，估计也没更好的方法了。。。
吹毛求疵一下，CDF(Cumulative Distribution Function)翻译为累积分布函数是不是更合适一点？
呃，似乎确实是这样看起来比较自然一些！
非常好的教程，
逆累积分布函数方法还是比较常用的，我们经常用
result = mu + delta*tan(pi*( rand – 0.5 ));
来产生柯西分布的随机数，
期待LZ的后续教程，对MCMC方法比较感兴趣
厉害啊！终于有原创的高质量中文技术blog.
kid啊 我面twitter有问到这道题 让我写多维的正太随机数生成器 但是如果各个维度不相互独立怎么办？
不独立的 Gaussian 可以通过旋转得到独立的 Gaussian，你可以生成独立的 Gaussian 然后再逆向旋转回去。
ps: 你要去实习？
嗯 打算暑假去体验下码农。。。
有一个地方没看明白，还请博主解答！为什么求r的概率分布函数的时候是对1/2pi*e^(r^2/2)*r求积分呢？应该是漏了一个负号吧，还有最后的那个r是怎么来的呢？
多谢指出，确实少了一个负号，那个“多余”的 r 是由于概率密度变换的时候后面要乘以导数的绝对值得到的（这里是多维的情况所以实际上是 Jacob 矩阵的行列式的绝对值）
博主，我最近遇到一个问题。我有一个二维的密度函数f(x,y)=C*(x-0.2)*(0.5-x)*(y-0.5)*(0.8-y) 边缘密度函数为f(x)=根号C*(x-0.2)*(0.5-x) f(y)=根号C*(y-0.5)*(0.8-y). 现在我应该怎么求出一个二维随机数呢
先算条件密度函数，f(y|x)，然后生成x ~ f(x)，再根据x和f(y|x) 生成 y ~ f(y| x ~ f(x))
楼主讲解的Inverse Transform Sampling Method 和正态分布的Box Muller原理非常好。看了很多相关的博客，基本上这两块都是按照pattern recognition and machine machine书里写的，没有讲背后的原理。
赞一个！