编写: @_散沙_Python玩家_

Concrete Math with Python II

Python带你游览具体世界的数学基础

一般来说,高深的数学总会让人有种抗拒心理,但跟具体世界紧密联系的简单数学知识会则不容易让人有如下的疑问:

"我们知道这些知识到底有什么用处?"

这些数学知识第一次出现在高中、大学课本时可能会让人摸不到头脑,而他们往往是基础的计算机科学、统计学、经济学知识的基石。

这部分内容包括:

  • 数列与金融基础
  • 排列组合与概率基础
  • 方程组/矩阵
  • 线性规划
  • 质数与加密
  • 随机数与概率分布

注意:具备扎实的高中数学、理工科本科数学基础的朋友可以简单阅览或者跳过这一部分。

1 排列与组合基础

(Permutations & Combinations)


In [16]:
import numpy as np
import scipy as sp
import scipy.special as spsp
import math
import itertools as itt
import bokeh as bk
import bokeh.plotting as bplt
import bokeh.palettes as bpal

In [13]:
# 工具配置
TOOLS = "pan,wheel_zoom,box_zoom,box_select"
# 作图配置
plot_config = dict(plot_width=800, plot_height=500, tools=TOOLS)
bplt.output_notebook()


BokehJS successfully loaded.

1.0 接受挑战!

“四个人打扑克(52张牌),每人随机派发13张,如果不考虑花色,其中某个人的手牌一共有多少种不同的可能性呢?”

1.1 从排列问题开始

permutations返回一个iterable,迭代中遍历所有的排列可能性。

把它放入内存:


In [2]:
a = [i for i in itt.permutations(xrange(3),3)]

In [3]:
print a,len(a)


[(0, 1, 2), (0, 2, 1), (1, 0, 2), (1, 2, 0), (2, 0, 1), (2, 1, 0)] 6

重要概念:阶乘 —— 全排列计数

$N! = N * (N-1) * ... * 1 = \prod_{i=1}^N i$


In [4]:
a = [i for i in itt.permutations(xrange(5),5)]

print len(a)


120

In [5]:
print math.factorial(5)


120

再次遇到迭代问题:


In [6]:
def myFactorial(N):
    try:
        intN = int(N)
        if intN in [0,1]:
            return 1
        elif intN<0:
            return None
        else:
            return intN*myFactorial(intN-1)
    except:
        return None

In [7]:
print myFactorial(5)
print myFactorial('a')
print myFactorial(-9999)


120
None
None

从N里挑K个排列?

先挑选排头兵,再挑选第二个……挑选k次。

$$P_N^k = N * (N-1) * ... * (N-k+1) = \prod_{i = N-k+1}^{N} i$$

In [8]:
a = [i for i in itt.permutations(xrange(4),2)]

print a,len(a)


[(0, 1), (0, 2), (0, 3), (1, 0), (1, 2), (1, 3), (2, 0), (2, 1), (2, 3), (3, 0), (3, 1), (3, 2)] 12

重新定义$P_N^k$:

$$P_N^k = N * (N-1) * ... * (N-k+1) = \prod_{i = N-k+1}^{N} i = \frac{\prod_{i = 1}^N i}{\prod_{i = 1}^{N-k} i}$$

试验一下$P_5^3$:


In [9]:
print len([i for i in itt.permutations(xrange(5),3)])
print math.factorial(5)/math.factorial(2)
print 5*4*3
print spsp.perm(5,3)


60
60
60
60.0

1.2 组合问题

当不关心被选择对象的顺序特征,而只关心选择模式的时候,这就变成了一个组合问题。

5个人中选择2个人一起完成一项任务(选择A和B与选择B和A视为同一种选择),有几种方法呢?

使用combinations:


In [10]:
a = itt.combinations(xrange(5),2)

In [11]:
la = [i for i in a]
for item in la:
    print item,

print '\n',len(la)


(0, 1) (0, 2) (0, 3) (0, 4) (1, 2) (1, 3) (1, 4) (2, 3) (2, 4) (3, 4) 
10

5选2问题可以用$C_5^2$表示也可以用${ 5 \choose 2 }$表示。

如果先不管顺序从N个元素中选取出k个,选出之后突然又关心起k个元素彼此的顺序(此时又有1组元素对应k!种可能的关系)

那么组合问题和排列问题的联系:$C_N^k * k! = P_N^k$,进一步有

$$C_N^k = \frac{P_N^k}{k!} = \frac{N!}{k!(N-k)!} = C_N^{N-k}$$

5选2的数目等于5选3的数目。


In [12]:
print myFactorial(5)/myFactorial(3)/myFactorial(2)
print spsp.comb(5,2),spsp.comb(5,3)


10
10.0 10.0

给定N,到底什么k使得${N}\choose{k}$最大呢?

提示:参考${N}\choose{k+1}$ $= \frac{N-k}{k+1}$ ${N}\choose{k}$,能得到什么结论?


In [17]:
p1 = bplt.figure(title="给定N时,组合数的取值",
       background_fill="#FFFFFF",**plot_config)

colors = bpal.brewer['RdBu'][4]

N = range(6,10)

def binomialmax(N):
    x = np.arange(0,N+1)
    def Nchoose(k):
        return spsp.comb(N,k)
    y = np.frompyfunc(Nchoose,1,1)(x)
    p1.line(x, y, line_width=3, color=colors[N-6], legend='N={0}'.format(N))

for n in N:
    binomialmax(n)
    

p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = 'k'
p1.yaxis.axis_label = 'N选k'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1

bplt.show(p1)


组合数一般有多大

阶乘很大,组合数一般没有阶乘那么大:


In [18]:
print math.factorial(25)
print 2**25
print spsp.comb(25,12)


15511210043330985984000000
33554432
5200300.0

新问题:阶乘一般有多大

Stirling公式指出:$$\lim_{n\rightarrow\inf} {n!}/{(\frac{n}{e})^n} \sqrt{2\pi n}= 1$$

试验一下:


In [20]:
N = np.array(xrange(30))+1

def stirling_convergency(n):
    return math.factorial(n)/(n/math.e)**n/math.sqrt(2*n*math.pi)

y = np.frompyfunc(stirling_convergency,1,1)(N)

p1 = bplt.figure(title="给定N时,组合数的取值",
       background_fill="#FFFFFF",**plot_config)

p1.line(N, y, line_width=3, color='steelblue', legend='Stirling公式的收敛性')

p1.grid.grid_line_alpha=0.3
p1.xaxis.axis_label = 'Stirling公式左手侧'
p1.yaxis.axis_label = 'N'
p1.ygrid.band_fill_color="steelblue"
p1.ygrid.band_fill_alpha = 0.1

bplt.show(p1)


因此 $${{N}\choose{N/2}}\approx 2^n \sqrt{\frac{2}{\pi n}}$$

组合数就是二项式系数

我们进行普通的代数式连乘展开计算,并且合并同类项,此时:

$$(a+b)^N = (a+b)(a+b)...(a+b) = \sum_{k=0}^N C_{N}^{k} a^{k}b^{N-k}$$

当$a = b = 1$时……简单验证一下,实际上$a,b$取任何值都满足


In [21]:
N = 15
sum = 0
for i in xrange(N+1):
    sum += spsp.comb(N,i)
print sum, 2 ** N


32768.0 32768

现实意义: $$ {N \choose 0} + {N \choose 1}+...+{N \choose {N-k}}+{N\choose N} = 2^N$$

等价与取子集问题:N个元素组成的集合,应该有多少子集?(包括本身的全集,以及没有元素的空集)

  • 当关心每一个元素在子集中的状态时:每个元素都可以有2种状态:在或者不在 —— $2^N$ 种子集

  • 关心每个子集包含的元素数:按元素数$0,1,...,N$分类,对每一类分别计数

1.3 其他有趣、有用的性质

1.3.1 帕斯卡等式

Pascal Identity:

$${{N}\choose{k}}+{{N}\choose{k+1}}= {{N+1}\choose{k+1}}$$

随机找几个数字试验一下:


In [24]:
for i in xrange(10):
    N = np.random.randint(20)+2
    k = np.random.randint(N-1)
    print "{0:>3} {1:>3} ".format(N,k),spsp.comb(N,k)+spsp.comb(N,k+1) == spsp.comb(N+1,k+1)


  7   1  True
 19   4  True
  6   1  True
 12   1  True
 19  16  True
  9   2  True
 18  11  True
 13   8  True
 17  15  True
 12   9  True

现实意义:

N+1个人要挑选一个k+1个人的小组去执行任务,而你恰好是其中一个。

  • 如果这个小组带着你去:$C_N^k$
  • 如果这个小组不带你去:$C_N^{k+1}$

1.3.2 杨辉三角形

然而杨辉三角形早已看穿了一切:第N排的第k个数字(N与k都从0开始)代表$N \choose k$

我们可以注意到三角形的两侧全都是1,即${N \choose 0} = {N \choose N} = 1$

帕斯卡等式:中间的某个数字来自于它左上方与右上方元素的相加。


In [30]:
# 为了排版美观,请让N<=10
def YangHui(N):
    for i in range(N+1):
        print '{0: ^55}'.format(' '.join(['{0: >4}'.format(np.int64(spsp.comb(i,j))) for j in range(i+1)]))
        print '\n',
# Pascal Identity的递推版本
YangHui(10)


                            1                          

                          1    1                       

                       1    2    1                     

                     1    3    3    1                  

                  1    4    6    4    1                

                1    5   10   10    5    1             

             1    6   15   20   15    6    1           

           1    7   21   35   35   21    7    1        

        1    8   28   56   70   56   28    8    1      

      1    9   36   84  126  126   84   36    9    1   

   1   10   45  120  210  252  210  120   45   10    1 

现在验证公式:$${k \choose k}+{k+1 \choose k}+...+{N \choose k} = {N+1 \choose k+1 }$$

首先把${k \choose k}$ 替换为 ${k+1 \choose k+1}$,实际上这就是杨辉三角形斜向一列相加,反复使用帕斯卡等式。

再次用代码验证:


In [31]:
for i in xrange(10):
    N = np.random.randint(20)+2
    k = np.random.randint(N-1)
    s = 0
    for j in xrange(k,N+1):
        s += spsp.comb(j,k)
    print "{0:>2}  {1:>2}  ".format(N,k), s == spsp.comb(N+1,k+1)


 4   1   True
13   7   True
18   3   True
11   4   True
 2   0   True
20   8   True
 7   3   True
11   8   True
14   6   True
 3   0   True

现实意义:“最大编号”问题

  • 等式右侧代表N+1个人中选出k+1个人:$N+1 \choose k+1$
  • 等式左侧:先假设选出k+1个人后,所有人中的最大编号是x,x可以取值$k+1,k+2,...N+1$,分别对x的取值分类。剩下的k人...

注意:如果用该公式同时除以等式右手边,就得到了随机分布中负二项分布的形式。

1.3.3 蛙跳

每2项取数加和:

$$I = {N \choose 0}+{N \choose 2}+{N \choose 4} + {N \choose 6}+...$$

考虑$(1+x)^n$,当$x=-1,1$的时候展开相加……

$$(1+1)^N = \sum_{k=0}^N C_{N}^{k} 1^{k}1^{N-k}$$$$(1-1)^N = \sum_{k=0}^N C_{N}^{k} (-1)^{k}1^{N-k}$$

k为奇数的时候被抵消了,而剩下的部分:

$$2^N = 2I$$

In [20]:
def step2binom(N):
    s = 0
    for i in xrange(0,N//2+1):
        s += spsp.comb(N,i*2)
    print "{0:>2}".format(N), s == 2**(N-1)

for i in xrange(5):
    N = np.random.randint(50)+2
    step2binom(N)


28 True
42 True
32 True
44 True
19 True

每4项取数加和:

$${N \choose 0} + {N \choose 4}+{N \choose 8}+ {N \choose 12}+...$$

考虑$(1+x)^n$,当$x=1,i,-1,-i$的时候……


In [21]:
def step4binom(N):
    s = 0
    for i in xrange(0,N//4+1):
        s += spsp.comb(N,i*4)
    print "{0:>2}".format(N) , np.int(s)  == np.int((2**N + (1+1j)**N + (1-1j)**N).real/4.0)

for i in xrange(10):
    N = np.random.randint(50)+4
    step4binom(N)


18 True
 9 True
35 True
 6 True
20 True
30 True
39 True
29 True
10 True
40 True

1.3.4 盲人摸球

$${M+N \choose k} = \prod_{i=0}^k {M \choose i} {N \choose k-i}$$

上式表明:

假设现在有一个袋子中放置有M颗红色的球,N颗绿色的球,而我们一把抓出来k颗球。

那么我们可以根据k颗球中红色球的数目x分类,x可以取值$0,1,2,...k$,从而得到如上公式:


In [30]:
for i in xrange(5):
    k = np.random.randint(10)+1
    M = np.random.randint(10)+1+k
    N = np.random.randint(10)+1+k
    s0 = spsp.comb(M+N,k)
    s1 = 0
    for j in xrange(k+1):
        s1 += spsp.comb(M,j)*spsp.comb(N,k-j)
    print "{0:>2}  {1:>2}  {2:>2}".format(M,N,k),s0 == s1


 5  13   4 True
14  17   8 True
13   8   7 True
 8  11   2 True
14  15   9 True

现实意义:

将上式两侧同时除以左手边,就能得到随机分布中的“超几何分布”的形式。

1.4 实战:卡牌大师

真正的卡牌大师,都是计算排列组合的高手。

Challenge:“四个人打扑克(52张牌,不含Joker),每人随机派发13张,如果不考虑花色,其中某个人的手牌一共有多少种不同的可能性呢?”

1.4.1 具体牌型


In [33]:
# 一种可能的组合
x = [1,1,2,4,5,6,6,7,10,10,10,11,12]
print len(x)


13

1.4.2 手牌库的变形

字典这个数据结构是很适合做统计的。

列表变字典:


In [33]:
xdict = {i:0 for i in xrange(13)}
for card in x:
    xdict[card]+=1
print xdict.items()


[(0, 0), (1, 2), (2, 1), (3, 0), (4, 1), (5, 1), (6, 2), (7, 1), (8, 0), (9, 0), (10, 3), (11, 1), (12, 1)]

字典变列表:


In [34]:
print xdict.values()


[0, 2, 1, 0, 1, 1, 2, 1, 0, 0, 3, 1, 1]

进一步,组合成一个字符串:


In [25]:
cardpool = []
translist = [cardcount*"*" for cardcount in xdict.values()]
transstr = "|".join(translist)
print transstr

finalstr = tuple([idx for idx,i in enumerate(transstr) if i == '|'])
print finalstr


|**|*||*|*|**|*|||***|*|*
(0, 3, 5, 6, 8, 10, 13, 15, 16, 17, 21, 23)

一组手牌、一组打印出来的变形后的字符串和最后那个元组形成一一对应。

  • 字符串永远有12个竖线,代表0-12这13种牌被12个竖线隔成13个区域
  • 第k个区域里有多少个*,就代表着这个人有多少张k号牌
  • 一个人共13张牌,所以字符串总长永远是25
  • 12个竖线的位置,唯一决定了牌型

1.4.3 卡池有限深

也许你已经发现了两件事:

  • 问题转换为:这不就是25选12吗?元组和combinations里的每一项很像!
  • 哪里不对?卡池有限深。每一张牌在整个卡组里毕竟最多只有4张!

问题转化为:只要存在任何一个区域中大于4星,就无法和真实情况对应了,应该排除。

  • 边界上超过4星
  • 任意两个竖线之间超过4星

In [2]:
import numpy as np
import itertools as itt
s = 0
for item in itt.combinations(xrange(25),12):
    if np.any(np.diff(np.array([-1]+list(item)+[25])) > 5):
        continue
    else:
        s += 1
print s


3598180

1.4.4 对动态规划有兴趣?

刚才的代码用动态规划方法求解的话,明显可以保持较少的代码与更高效的运行:

注:res[i,j]表示:在只能选择j张不同的点数的牌,且相同点数的牌不能超过4张的时候,一共抽i张牌所能获得的不同牌型(不考虑不同花色)


In [9]:
import numpy as np
res = np.zeros((14,14))
for i in xrange(14):
    res[0,i] = 1
for j in xrange(1,14):
    for i in xrange(1,14):
        if i>4*j:
            res[i,j] = 0
        else:
            for k in xrange(min(i+1,5)):
                res[i,j] += res[i-k,j-1]
print  res[13,13]


3598180.0

1.5 There's more...

排列组合在算法面试里面是相当常见的问题了,来自测一下以下问题。如果都能够很快得出问题的答案说明过关啦!

  • 一个人在无限长X轴横轴上随机游走,n步以后回到原点,有多少种不同走法?
  • M*N的格子从左下角走到右上角,只允许向右走或者向上走,有多少种不同的走法?
  • M*N的格子,里面有多少个矩形?(正方形和长方形都算!)
  • 有k个不同的盒子,其中放入N个没有差异的小球,问不同的放法有多少。

对于排列组合有更深入兴趣的,建议阅读以下附加材料:

图论里也有很多排列组合问题,进一步深入的读者可以自由阅读以下内容

2. 真$\cdot$概率论初步

2.0 概念

概率空间三要素: $(\Omega,\mathcal{F},P)$

  • 样本空间$\Omega$,由样本点$\omega$组成
  • 在其上观测到的事件集 $\mathcal{F}$
  • 以及为事件集分配的一个从0到1的概率值 $P$

$\Omega,\omega$的理解:骰子

事件集$\mathcal{F}$的理解和性质:

  • 如果$\mathcal{F}$存在事件$A$,那么$\bar{A}$也就是“非A”也存在
  • 如果$\mathcal{F}$存在事件$A,B$,那么$A\cap B$和$A\cup B$也存在
  • $\mathcal{F}$存在全集和空集

$P$是为了$\mathcal{F}$中的元素服务的而不是为$\Omega$中的元素$\omega$服务的

独立:$P(A\cap B) = P(A)*P(B)$

互斥:$P(A\cup B) = P(A)+P(B), P(A\cap B) = 0$

包含:$P(A \cup B) = P(A), P(A \cap B) = P(B)$

2.1 无非就是计数问题+排列组合

赌博问题:独立、公平的骰子点数

  • 2颗骰子>7?
  • 3颗骰子>13?
  • 3颗骰子2奇数1偶数?

In [27]:
#Q1 List Them
l = [ (x, y) for x in xrange(6) for y in xrange(6) if x+1+y+1 >7 ]
print len(l) ,len(l) / 6.0 ** 2


15 0.416666666667

In [28]:
#Q2 List Them ... Lazily
g = ( (x, y ,z) for x in xrange(6) for y in xrange(6) for z in xrange(6) if x+1+y+1+z+1 > 13 )
cnt = 0
for x in g:
    cnt += 1
print cnt / 6.0 ** 3


0.162037037037

In [4]:
#Q3 Independence
p = (3/6.)**3 * spsp.comb(3,2)
print p


0.375

2.2 经典的宝藏问题

如果一个概率问题让你很头痛,那么:

搞清楚划分事件的集合、以及其上的概率分布是如何均匀分布的十分重要。

问题:3道门后有且仅有1个门有宝藏。随机选定某个门之后,主持人帮忙去掉一个无宝藏的门。换还是不换?

我们可以让编号0代表有宝藏,1号门后有垃圾,2号门后有废铁(怎样都好啦...)

换一个思路:现在3个门是随机打乱的,而我们选中最前面的门。(这和门固定下来我们随机选是一样的效果。)


In [29]:
for item in itt.permutations([0,1,2]):
    print item


(0, 1, 2)
(0, 2, 1)
(1, 0, 2)
(1, 2, 0)
(2, 0, 1)
(2, 1, 0)

这6种情况应该是等概率的,因此每一种情况占据1/6的概率。

后4种情况下,主持人一定会帮忙去掉2或者1而剩下0;前2种情况下,只要换门就是错的。因此总体而言,换门更合适。

2.3 回到卡牌上来

如果现在问你,一局4人卡牌游戏中,某个人抽到4张A的概率?


In [35]:
# 含有4张A的牌型比上所有牌型的总数
import numpy as np
import itertools as itt
s1 = 0
for item in itt.combinations(xrange(20),11):
    if np.any(np.diff(np.array([-1]+list(item)+[20])) > 5):
        continue
    else:
        s1 += 1
print s1/3598180.


0.0421268530201

这个结果和下面的结果哪个对?


In [36]:
# 使用超几何分布
print spsp.comb(4,4)*spsp.comb(48,9)/spsp.comb(52,13)


0.00264105642257

后者对,为什么?特定牌型之间的概率并不相等,观测到不同牌型是观测到不同事件集,在其上分配的概率不等。

比如当牌型为x = [0, 2, 1, 0, 1, 1, 2, 1, 0, 0, 3, 1, 1]时,

$P(x) = \prod_{i=0}^{12}$ ${4}\choose{x[i]}$ / ${52}\choose{13}$

并不等于1/3598180


In [37]:
x = [0, 2, 1, 0, 1, 1, 2, 1, 0, 0, 3, 1, 1]

px = 1

for i in x:
    px *= spsp.comb(4,i)
    
print px/spsp.comb(52,13)
print 1/3598180.0


9.28836858809e-07
2.77918280909e-07

条件概率

给定A事件时发生B事件的概率,就是在A条件下B的条件概率,记为$P(B|A)$

事件A:抽13张牌(52张牌每张都不同,样本点而非牌型$\omega$构成了概率空间$\Omega$)

事件B:抽13张牌时,抽到4张K

事件C:抽13张牌时,抽到3张K

事件C':抽13张牌时,抽到3张K或以上

$P(B\cap C) = 0, P(C') = P(B \cup C) = P(B) + P(C)$

$P(A) = 1, P(B) = P(AB) = \frac{card(B)}{card(A)}$

更多的,如果告诉你抽13张牌时,至少抽到3张K了(C'),那么此时抽到4张的概率?

$$P(B|C') = \frac{P(B C')}{P(C')} = \frac{P(B)}{P(B)+P(C)}$$

In [2]:
pb = spsp.comb(48,9)*spsp.comb(4,4)
pc = spsp.comb(48,10)*spsp.comb(4,3)
print pb/(pb+pc)


0.0602409638554

3 小结

本次我们深入了解了排列组合中的问题。

将日常中的问题翻译为排列组合问题之后,我们对于解决简单的概率问题(古典概率问题)就会更加得心应手。

以及前面的面试题答案:

  • ${N}\choose{N/2}$
  • ${N+M}\choose{M}$
  • ${M+1}\choose{2}$ $\cdot$ ${N+1}\choose{2}$
  • ${N+k-1}\choose{k-1}$

In [ ]: