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()
permutations返回一个iterable，迭代中遍历所有的排列可能性。
把它放入内存：
In [2]:
a = [i for i in itt.permutations(xrange(3),3)]
In [3]:
print a,len(a)
重要概念：阶乘 —— 全排列计数
$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)
In [5]:
print math.factorial(5)
再次遇到迭代问题：
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)
从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)
重新定义$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)
当不关心被选择对象的顺序特征，而只关心选择模式的时候，这就变成了一个组合问题。
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)
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)
给定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)
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}}$$
In [21]:
N = 15
sum = 0
for i in xrange(N+1):
sum += spsp.comb(N,i)
print sum, 2 ** N
现实意义： $$ {N \choose 0} + {N \choose 1}+...+{N \choose {N-k}}+{N\choose N} = 2^N$$
等价与取子集问题：N个元素组成的集合，应该有多少子集？（包括本身的全集，以及没有元素的空集）
当关心每一个元素在子集中的状态时：每个元素都可以有2种状态：在或者不在 —— $2^N$ 种子集
关心每个子集包含的元素数：按元素数$0,1,...,N$分类，对每一类分别计数
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)
现实意义：
N+1个人要挑选一个k+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)
现在验证公式：$${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)
现实意义：“最大编号”问题
注意：如果用该公式同时除以等式右手边，就得到了随机分布中负二项分布的形式。
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)
每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)
上式表明：
假设现在有一个袋子中放置有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
现实意义：
将上式两侧同时除以左手边，就能得到随机分布中的“超几何分布”的形式。
In [33]:
# 一种可能的组合
x = [1,1,2,4,5,6,6,7,10,10,10,11,12]
print len(x)
In [33]:
xdict = {i:0 for i in xrange(13)}
for card in x:
xdict[card]+=1
print xdict.items()
字典变列表：
In [34]:
print xdict.values()
进一步，组合成一个字符串：
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
一组手牌、一组打印出来的变形后的字符串和最后那个元组形成一一对应。
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
刚才的代码用动态规划方法求解的话，明显可以保持较少的代码与更高效的运行：
注：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]
排列组合在算法面试里面是相当常见的问题了，来自测一下以下问题。如果都能够很快得出问题的答案说明过关啦！
对于排列组合有更深入兴趣的，建议阅读以下附加材料：
图论里也有很多排列组合问题，进一步深入的读者可以自由阅读以下内容
概率空间三要素： $(\Omega,\mathcal{F},P)$
$\Omega,\omega$的理解：骰子
事件集$\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)$
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
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
In [4]:
#Q3 Independence
p = (3/6.)**3 * spsp.comb(3,2)
print p
In [29]:
for item in itt.permutations([0,1,2]):
print item
这6种情况应该是等概率的，因此每一种情况占据1/6的概率。
后4种情况下，主持人一定会帮忙去掉2或者1而剩下0；前2种情况下，只要换门就是错的。因此总体而言，换门更合适。
如果现在问你，一局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.
这个结果和下面的结果哪个对？
In [36]:
# 使用超几何分布
print spsp.comb(4,4)*spsp.comb(48,9)/spsp.comb(52,13)
后者对，为什么？特定牌型之间的概率并不相等，观测到不同牌型是观测到不同事件集，在其上分配的概率不等。
比如当牌型为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
给定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)
In [ ]: