• 0x00 前言
• 0x01 从$\pi$的奇葩算法谈起
  • 内切圆
  • 蒲丰投针
• 0x02 基数计数(Cardinality Counting)
  • bitmap与Hash函数
  • Linear Counting
  • LogLog Counting
• 0x03 总结
• 0x04 参考文献

0x00 前言

一直以为，概率与统计是一门非常有传奇色彩的学科。它由赌博业催生而来，但你确不能靠它在赌场发家致富。它不能告诉你抛一次硬币是正是反，但它能告诉你抛1万次硬币正反面大概各5000次。而现实生活中的概率并非简单的抛硬币，我们所看到的并非真的是随机，而那些真正的随机又常常难以理解。

最近知乎上有一篇很不错的科普，《伪随机的上位和真随机的逆袭》，推荐一看~

0x01 从$\pi$的奇葩算法谈起

圆周率$\pi$我们每个人都很熟悉，说起圆周率的计算，我们最熟悉的版本无非是祖冲之的“割圆法”，当然随着科学的不断进步，数学界的牛人们不断刷新着$\pi$的精度，什么贝利－波尔温－普劳夫公式（BBP公式）、高斯-勒让德算法（GLA）、快速傅里叶变换、牛顿迭代法……

你如果对自己的数学充分自信的话可以看看《一起来算圆周率》，如果你和我一样已经完全抓瞎了，那么还是继续往下看吧╮(￣▽￣)╭

有没有正常人可以理解同时也非常简单粗暴的方法来算$\pi$呢？当然有，这就是传说中的蒙特卡罗方法(Monte Carlo method)，俗称统计模拟，也可以称之为熔火之心(MC)╮(￣▽￣)╭

说白了，这是一种用随机数或伪随机数来解决计算问题的方法，这话听上去非常矛盾，我们认知中的计算是非常精确的，怎么可以通过随机数来计算呢？别急，我们先来看看怎么用MC方法来计算$\pi$

内切圆

在平面上画一个正方形和一个内切圆，然后随机向正方形区域内撒$n$个点，统计落在圆内点的个数$m$，如下图所示

可以认为我们随机撒的点服从均匀分布，因此点的个数与平面区域的面积有如下关系

最终可以得到

你看，我们并没有用到什么高深的数学知识，只需要统计出区域内点的个数，利用上面的式子即可计算出$\pi$，而且我们撒出的点越多，计算出的$\pi$就越精确。用程序实现起来也非常方便

from random import random
n=10**6
print sum(1 if random()**2 + random()**2 < 1 else 0 for i in range(n))*4.0/n

蒲丰投针

再来看第二种利用统计模拟计算$\pi$的方法，这就是大名鼎鼎的蒲丰投针问题(Buffon’s needle)

在平面上画上无数条平行线，间距为$d$， 取一根长度为$l(l < d)$的针，随机向该平面投掷$n$次

然后统计与平行线相交的次数$m$，如下图所示

那么$\pi$可以由下式算出

方法是很简单，但是为什么呢？下面我们来证明一下。

首先我们把问题简化一下，只考察两条平行线之间的半区域（其他区域都是等价的）。容易观察出针与平行线是否相交取决于两个因素，针的位置和旋转的角度，我们用$x(0\le x \le \frac{d}{2})$表示针中心点到平行线的距离，用$\theta(0 < \theta < \pi)$表示针与平行线的夹角，如下图所示

显然，当$x \le \frac{l}{2}\cdot sin\theta$时，针与平行线相交。如下图所示

从图中可以看出，当$x$与$\theta$的取值落在区域$g$内时，上述条件即可满足，所以针与平行线相交的概率为

而$p=\frac{m}{n}$，所以

python实现

from random import uniform
from math import pi, sin
#注意:为了便于计算,这里我们假定l=d=1
def direct_needle(N):
   N_hits=0
   for i in range(N):
      x_center=uniform(0.,0.5)
      phi=uniform(0,pi/2)
      x_tip=x_center - sin(phi)/2.
      if x_tip < 0: N_hits += 1
   return N_hits
n=10**6
print 2.0*n/direct_needle(n)

这里还有一个关于蒲丰投针证明的视频，讲的很生动形象（台湾腔的大叔讲话还是很有感觉的╮(￣▽￣)╭）

与MC方法相关的还有许多理论与算法，如马尔可夫链蒙特卡洛方法(Markov Chain Monte Carlo，简称MCMC)，吉布斯采样(Gibbs Sampling)……这些方法在物理、经济等诸多领域发挥着重要作用~当然本文的重点并不是介绍MC方法，只是想通过这个例子说明，在理论和事实之间有着某种神秘的联系（就像我们可以从大量统计数据中估算出$\pi$一样），而概率与统计就是我们认识和研究这种联系的利器。

ps:后面的内容较长且涉及的数学知识较多，理解起来可能不太轻松，感兴趣的同学可以把本文存为书签，有时间的时候再看~

0x02 基数计数(Cardinality Counting)

基数(cardinality，也译作势)，是指一个集合中不重复元素的个数。注意这里的集合和我们学过的严格定义的集合不同，允许存在重复元素，另外，本文所讨论的均为有限集。如给定这样的一个集合$\left \{ 1,2,3,1,2 \right \}$，它有$5$个元素，但它的基数为$3$。

基数计数在很多领域都有应用，最为的典型同时也最为大家所熟知的就是网站统计了。如果把PV（Page View， 即页面浏览量或点击量，用户每次刷新即被计算一次）看作一个集合，那么UV（Unique Visitor，访问该页面的一台电脑客户端为一个访客，也就是我们常说的独立访客）就是这个集合基数。显然UV是一个页面价值和影响力的一个重要评判指标，那么如何对其计数呢？

有人说这还不简单，定义一个包含独立访客标识，访问计数等信息的数据结构，最后用数组或者链表保存起来不就行了。每当有新的访问行为发生，则依据其访客标识进行查找，命中则说明该访客已访问过该链接，访问计数加一；若未命中，说明该访客是首次访问，则分配一个新空间保存该访客信息。

如果只是一个小型网站这么做是完全可以的，但是如果是一个每日PV在上亿甚至更高数量级的大型网站中，线性表的查询效率可想而知，而进行二分查找、B树等优化也会遇到各种各样的问题。这里我不想炒大数据这个概念，但在现实场景中，数据到达一定数量级之后这都是不得不考虑的问题。

bitmap与Hash函数

bitmap是非常经典的海量数据处理数据结构，其本质是用bit数组的某一位表示某一数据，从而一个bit数组可以表示海量数据。用0表示某一元素不在集合中，用1表示某一元素在集合中，如0100000011000000可以用来表示集合$\left \{ 1,8,9 \right \}$。

Hash函数，也称散列函数，原本是密码学领域的概念。说白了就是将一长串输入转换为固定长度的输出，即$H(M)=h$，其中$H$为Hash函数，$M$为密文，$h$为定长的hash值。同时要满足以下几个条件：

• 单向性
• 抗冲突性
• 映射分布均匀性和差分分布均匀性

更多细节内容就与本文无关了，这里不再赘述。但是在数据结构中我们用到的Hash函数并没有那么严格的要求，一般满足第二点就可以了，有的时候甚至简单的如求余函数也可作为Hash函数。那么为什么要使用Hash函数呢？很简单，单凭其时间复杂度为$O(1)$便足以傲视群雄了，在海量数据处理中，这个优点显得尤为突出。

怎样将两者结合起来进行基数计数呢？比如我们选取md5作为哈希函数，只取其结果的前8位，那么整个Hash函数的映射空间就有$16^{8}$（40多亿，几乎不会产生冲突）个值，取一个长度为$16^{8}$的bitmap（大约536MB），每一位对应Hash函数映射空间中的一个值，初始值全为0。每当有新访问产生，对该访客标识进行Hash，并映射到bitmap中的某一位上，若该位置为0，则置1；若为1，则不作处理。最后统计整个bitmap中1的个数即为基数。

整个过程没有特别耗时的查找操作，大大提高了效率。但这种方法占用了大量空间，我们这里只是统计某一个链接的数据，如果要统计全站的链接的话就不行了，即便换一个映射空间较小的Hash函数也不能从本质上解决这个问题。

当然这是属于比较精确的计数，相应的代价自然会比较高。在实际应用中我们可能并不需要知道太过精确的数据，10000000的访问量和10000100的访问量并没有太大差别，有没有一种方便快捷的估算方法呢？

又到了概率与统计出场的时候了~

Linear Counting

这里我们先从简单的LC算法(Linear Counting)讲起，仔细分析上面的例子不难发现其空间占用较多是因为其过于追求Hash函数的抗冲突性，进而导致映射空间过大。LC算法正是大大降低了Hash函数的要求，并利用概率与统计的相关知识，最终给出基数的一个估计。

LC算法描述如下，设有一Hash函数H，其Hash映射空间有m个值（最小值0，最大值m-1），且服从均匀分布。取一个长度为m的bitmap，每一位与Hash函数映射空间中每个值一一对应，初始值全为0。设一个集合的基数为n，此集合所有元素通过Hash函数映射到bitmap中，如果某一个元素被Hash到第k个比特并且第k个比特为0，则将其置为1。当集合所有元素Hash完成后，设bitmap中还有u个bit为0。则：

为n的一个估计，且为最大似然估计（MLE）。

乍一看，这个算法的描述和前面的差不多，其实最大的区别在于这里的Hash函数是允许冲突的，也就是说允许$H(M)=H(M’)$的情况出现。

下面来看看如何进行证明。首先我们要明确一点，这里bitmap的最终值只与集合的基数有关，比如$A=\left \{ k_1,k_2,\cdots,k_n \right \}$和$B=\left \{ k_1,k_2,\cdots,k_n ,x_1,x_2,\cdots,x_i \right \}$两个集合，其中$k_1,k_2,\cdots,k_n$为$n$个不同元素，$x_i \in \left \{ k_1,k_2,\cdots ,k_n \right \}$。显然，$A$和$B$的基数是相同的，经过Hash函数映射得到的bitmap也是一样的。

所以下面的证明我们就用集合$A$来进行辅助思考~

LC算法定义的Hash函数是满足映射均匀分布的，也就是说集合$A$中的每一个元素映射到bitmap中的每一位都是等可能的（这里有的同学可能会有点困惑，Hash函数不是确定的么，为什么又和等可能扯上关系了呢？可以这样理解，在未进行Hash运算时，你是不知道输入值将映射到哪个位置的，而由于Hash函数的均匀分布性，所以其映射到bitmap任何位置都是等可能的，但是相同的输入值一定会映射到相同的位置上）。显然，在一次映射中，bitmap上的某一位被映射到的概率为$\frac{1}{m}$，不被映射到的概率为$1-\frac{1}{m}$，设事件$C_j$为“经过$n$次元素Hash后，bitmap上第$j$位为0的概率”，则有

由于bitmap上每一位都是独立的，所以$u$的期望为

如果你高数还没忘光的话，一定还记得下面这个式子

自然数$e$就是这么定义的~

所以，当$m$和$n$都趋于无穷时有，

bitmap上每一位的值服从参数相同0-1分布，因此u服从二项分布。由概率论知识可知，当n很大时，可以用正态分布逼近二项分布，因此可以认为当n和m趋于无穷大时u服从正态分布。

u的概率密度函数

由于我们观察到的0的个数u是从正态分布中随机抽取的一个样本，因此它就是$\mu$的最大似然估计（正态分布的期望的最大似然估计是样本均值）。

又由如下定理：

$f(x)$是可逆函数，$\hat{x}$是$x$的最大似然估计，则$f(\hat{x})$是$f(x)$的一个最大似然估计。

$f(x)=-m\ln\frac{x}{m}$是可逆函数，$u$为$E(u)$的最大似然估计，所以$\hat{n}=-m\log\frac{u}{m}$为$n=-m\log\frac{E(u)}{m}$的一个最大似然估计。

这里我们也要注意到，m的取值不能太小，不然很有可能bitmap上所有位都被映射到了，这样u就为0了，整个算法就失去了意义。原论文作者给出了一张表

可以看出m大约为n的十分之一左右。

LogLog Counting

从前文可以看出LC算法虽然有一些提升，但毕竟只是线性的减少了空间占用，如果面对上亿级别的基数计数，这个方法依然要占用大量空间。但这小节将要介绍的LogLog Counting的空间复杂度只有$O(\log_2(\log_2(N)))$，是的你没看错，取了两次$\log$。

在开始介绍LLC算法前，要提前说明的是这里的Hash函数又是要求比较严格的，要保证碰撞几率极小，同时映射空间是近似均匀分布的。同时，LLC算法中不使用bitmap这个数据结构。

设a为待估集合中的一个元素，$h=H(a)$，这里把h表示为长度为L的比特串（如h为098950fc，写为00001001100010010101000011111100），将这L个比特串从左至右依次编号为1、2、……、L。因为Hash函数是均匀分布的，所以这L个比特服从如下分布且相互独立

$ P(x=k) = \left \{_{0.5(k=1)}^{0.5(k=0)} \right. $

其实也就是说h中每一个比特位出现0或者1都是等可能且互相独立的，设$\rho(h)$为比特串h中第一个1出现的位置（例如h为098950fc，则$\rho(h)=5$）。现在对集合中的每个元素都作Hash运算，对每个生成的h都计算其$\rho(h)$，取$\rho_{max}$为所有$\rho(h)$中的最大值。此时，我们可以估计集合的基数为

为什么可以这样估计呢？

我们可以把上述寻找比特串中第一个1的过程看作一个投硬币试验：当硬币为反面时，记为0；当硬币为正面时，记为1，试验停止，记录投掷次数。设n次试验中，最大投掷数为k。

现在考虑如下两个事件：

A：进行n次试验，每次投掷次数都不大于k

B：进行n次试验，至少有一次投掷次数大于等于k

在一次试验中，投掷次数大于k的概率为$(\frac {1}{2})^k$，相当于只要连续投掷出k个反面，投掷次数肯定大于k，那么在一次试验中，投掷次数不大于k的概率为$1-\frac {1}{2^k}$，所以

B的对立事件为“进行n次试验，每次投掷次数都不大于k-1”，所以

注意到，当$n \gg 2^k$时，$P(A) \to 0$，当$n \ll 2^k$时，$P(B) \to 0$，转换为自然语言描述就是，当n远远大于$2^k$时，每次试验投掷次数都不大于k的概率几乎为0，当n远远小于$2^k$时，至少有一次试验投掷次数大于等于k的概率也为0。而这些均与我们观察到的试验结果不符，即存在至少一次试验的投掷次数等于k，且不存在比k更多的投掷次数。所以唯一合理的推断是$n \approx 2^k$。

现在回到我们的初始问题，自然有

当然在实际使用中，直接这样估计仍会存在较大误差，所以LLC算法还采用了分桶平均的方法（类似物理实验的多次测量求平均值），这里不再深入。

0x03 总结

当然这里介绍的基数计数方法都是一些基本的算法和思想，除此之外还有Adaptive Counting、HyperLogLog Counting等方法，相关链接会在参考文献中列出。

概率与统计的神奇之处还不仅限于这些，在统计学习中有更多让人着迷的应用。最后，引用《统计与真理》这本书的里的话作为结语~

All knowledge is, in final anlysis, history.

All sciences are, in the abstract, mathematics.

All judgements are, in their rationale, statistics.

如果你觉得本文对你有帮助，欢迎打赏我一杯咖啡钱，支持我写出更多好文章~

0x04 参考文献

• 解读Cardinality Estimation算法（第一部分：基本概念）
• 解读Cardinality Estimation算法（第二部分：Linear Counting）
• 解读Cardinality Estimation算法（第三部分：LogLog Counting）
• 解读Cardinality Estimation算法（第四部分：HyperLogLog Counting及Adaptive Counting）