没有看过上一篇的同学请看正态分布随机数的生成 (1)。
求反变换固然还可行，但是碰到无法解析求逆的函数，用数值方法总归比较慢。下面我们就来说说另一个能够适合任何概率密度分布的方法——接受—拒绝法 (Acceptance-Rejection Method)，国内也有翻译成叫做舍选法的。接受—拒绝法的思路其实很简单——比如说你想要正态分布，我们就弄个方框框把它框起来，然后均匀地往里面扔飞镖。扔到曲线以下我就留着，扔到曲线以上就不要了。这样搞好以后来看，曲线之下的点就是（二维）均匀分布的。那这些点的横坐标就正好满足我们要的分布——高的地方的点就多，低的地方的点就少嘛。
正态分布——听起来非常耳熟，用起来也很顺手——因为很多语言都已经内置了有关正态分布的各种工具。但其实，在这个最普遍、最广泛的正态分布背后，要生成它还有很多学问呢。
$$f(x \; | \; \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi} } \; e^{ -\frac{(x-\mu)^2}{2\sigma^2} }$$
难道教科书上没有讲吗？看看概率书上是怎么说的……比如我手头这本浙大版的《概率论与数理统计》（第四版）第 378 页上说……“标准正态变量的分布函数 $\Phi(x)$ 的反函数不存在显式，故不能用逆变换法产生标准正态变量……”
Continue reading “正态分布随机数的生成 (1)”
前面说的是均匀抽样，要是想加个权怎么办呢？先说加权有什么用呢？比如我们已经统计好了搜索的关键字和词频，那么有了加权就可以直接用这个数据来抽样而无需把关键字重复好多遍了。
我们先来看看这个抽样应该是怎么样的。假设这总共 $n$ 项都摆在你面前了，设第 $i$ 项的权值为 $w_i$，那么我们抽到这一项概率应该是
$$P_i = \frac{w_i}{w_1 + w_2 + \cdots + w_n}$$
我们现在来看看，如果把这 $n$ 项按照第 $n$ 项、第 $n-1$ 项、…… 第 $2$ 项、第 $1$ 项的顺序抽出来，这个概率是多少。因为这 $n$ 项的顺序是可以随便摆的，所以它求出来是有一定普遍意义的。
第一次抽出第 $n$ 项的概率是
$$P_n(1) = \frac{w_n}{w_1 + w_2 + \cdots + w_n}$$
第 $n$ 项已经被抽走了，那第二次抽到第 $n-1$ 项的概率是
$$P_{n-1}(2) = \frac{w_{n-1}}{w_1 + w_2 + \cdots + w_{n-1}}$$
以此类推，按照这个顺序抽出所有元素的概率就是
$$P(S) = \prod_{i=1}^n \frac{w_i}{w_1+w_2+\cdots w_i} $$
这有什么用呢？别着急，先记下这个结果，我们下面来说说算法怎么做。
没有看过前一篇的同学请看这里：蓄水池抽样浅说 (1)
前面介绍的方法都很好，但是要计算 $n-k$ 个随机数实在是有点浪费时间…… 有兴趣的同学可以把那个调用函数换成比如 reservoirSample(range(10**8), 10**5)，就知道这东西还是要算上一会儿的了。
再来看看这个过程吧——把前 $k$ 个元素放进去之后，随机决定接下来的元素要不要放进去，但是每次决定时都需要产生一个随机数。要是我们能够确定应该跳过多少个元素，不就可以省掉很多生成随机数的工夫了吗？每一步过程就变成了
我们有时候会有进行数据抽样的需要，比如要从文件中随机取出若干行，或从数据集中随机取出若干数据进行分析。通常情况下这并不是什么难事，比如 Python 中直接提供了 random.sample() 来做这件事，Numpy 中更提供了功能更为强大的 numpy.random.choice()。然而这些东西都有一个问题，就是你必须把整个数据集读到内存里。如果数据集超出了内存的限制，或者要对一个持续的输入流做抽样，即从包含 $n$ 个项目的集合中（等概率）选取 $k$ 个样本，其中 $n$ 为很大或未知，又该怎么做呢？
怎么求出从 $2$ 到 $N$ 之间的所有素数之和？
素数求和是数论中一个很典型的问题。这事情乍听起来并不难，只要把素数都求出来再加起来就好了嘛。可能你已经了解了一些求素数的方法，比如最简单的试除。写过这个程序的同学都知道这样搞有多慢。或者常用的是埃氏筛法，简单来说，就是列好从 $2$ 到 $N$ 的范围，先用 $2$ 去筛，把从 $2 \times 2$ 开始的 2 的倍数剔掉；再用下一个素数 3 筛，把 $3 \times 3$ 开始的 3 的倍数剔除掉…… 一直筛完不大于 $\sqrt{N}$ 所有素数为止。我直接从 wiki 上借了一个动画来演示这个过程，有兴趣的同学可以自己写一个试试看。
本文中的方法来自于解决了 ProjectEuler 全部问题的大神 Lucy_Hedgehog。
Continue reading “素数求和的动态规划方法”
前面的方法已经非常快了，应该说很满意。那还有没有再优化的余地呢？
我们再来仔细看一下这个过程：深度优先搜索，逐级回溯，典型的二叉树遍历方式。为了方便说明，我标出了结点遍历的顺序：
Continue reading “0-1 背包问题详解 (4)”
我们先把动态规划的想法放一下，退一步来看看这个问题。每一个物品都可能放或者不放，如果用一个决策树来表示的话，就得到这样一棵满二叉树：
这对于 \(n\) 个物品就有 \(2^n\) 种放法，总共要计算 \(2^{n+1} – 2\) 个结点，当然是太多了。不过，有没有什么办法能让遍历的结点少一点呢？
不知道有没有读者试过了这个大包包，我试了一下，在我的 2.30GHz i5 5300U 上，这个函数用 cPython 跑了…… 1632.07 秒，将近半个小时。照这个搞法，每次那个包有多大，就得搞出多少项。那个包的大小可是个任意数啊？这 $ W$ 比 $ 2^n$ 大都说不定呢？
Continue reading “0-1 背包问题详解 (2)”
所有代码可见于：https://github.com/jameslao/Algorithmic-Pearls/tree/master/0-1-Knapsack
这些程序大部分都是 2013 年跟着 Tim Roughgarden 看算法课的时候写的，还算有些心得，尤其是优化的部分，似乎中文书中讲到的不是很多。而且也很久没有写博客了，回顾一下当初学过的东西也还有些收获。
不知道你还记不记得我们小时候看过的《太阳山》的故事：
从前，有弟兄俩，老大是个贪心的富人，老二是个穷人。
老二没有地，种老大的地。他一年到头在地里做活；可是老是吃不上，穿不上，租子交不够。
有一天，老二上山打柴。天都快黑了，他还坐在山上；他发愁，加上打柴，租子还是交不够。这时候，忽然刮起一阵黑风，飞来一只大鸟。大鸟落在他跟前，对他说：“你不用发愁。太阳山上有很多金子、银子，还有宝石。我背你上太阳山去吧！你从那儿可以拿些金子回来。”
老二爬到大鸟的背上。大鸟叫他闭上眼睛。他刚闭上眼睛，觉得飕一下子，就听见大鸟说：“到了。你看，有多少金子！去拿吧！你可不要贪心。这是太阳山。如果我们在这儿待得时间太长，赶上太阳回来，就把你烧死了。那时候，我也没办法救你。”
老二一边答应着，一边从大鸟的背上跳下来。哎呀！遍地是金子、银子、宝石，金光闪闪，晃得人睁不开眼睛。老二想了想，就拿一小块金子，装进口袋，要大鸟背他回去。大鸟问他为什么只拿一小块儿。他说够了。
老二又爬到大鸟的背上，闭上眼睛。飕一下子，大鸟落到他家门口了。他谢过大鸟，大鸟就飞走了。
老二有了金子，买了地，盖了房子。他早晨起来去种地，晚上回到新盖的房子里休息，日子过得很好。
老大见老二有金子，很眼红。他问老二的金子是哪儿来的。老二把大鸟背自己到太阳山的事告诉他了。
第二天，老大学着老二到山上去打柴。天快黑了，他故意不回家，坐在山上假装发愁。大鸟飞来了，也答应背他上太阳山去拿些金子。
大鸟背老大到了太阳山，嘱咐他不要贪心，像嘱咐老二一样，可是他连哼也不哼。他看见那么多金子、银子、宝石，忙着张开大布袋子，捡最大块儿的金子往里装。
老大拼命的往大布袋子里装金子，装个没完。大鸟催他回去，他说再让他拿几块。大鸟再催他回去，他还是这么说。最后，大鸟说：“时候到了！太阳马上就回来了！”他这才站起来，朝着大鸟走，摇摇晃晃的，大布袋子压得他走不稳了。没走几步，他又弯下腰来，说：“等一会儿。我再拿几颗宝石吧！”
正在这时候，火红的太阳回来了。它把烈火似的阳光喷到太阳山上。大鸟飞走了。老大被烧死了。
要是换你去拿，你会怎么拿呢？你是只拿一小块，还是拿到自己走不动？有没有什么办法，能够拿到最多的金子呢？
Continue reading “0-1 背包问题详解 (1)”
