代码中的概率学（从理论到仿真的实践）

某日，我们几个同事在玩一个叫Splendor（璀璨宝石）的游戏。在翻贵族牌的时候，出现了一个极致的情况，连续翻出5块3x3宝石组合的牌。在众人惊呼这个奇迹的时候，突然有人提出，“这个概率有多少啊？”。于是，一帮理工科的死宅，开始了脑区的激活。

这个贵族牌的情况大致是这样的：

有两类牌，一类是2x4宝石组合，另一类是3x3宝石组合。

每类牌各有5张，一共10张。

一局游戏一次随机翻出其中5张牌。

可以看到，总共只有5张3x3宝石牌，如果需要翻出5张3x3宝石牌，那就意味着，每次随机翻牌时必须全部准确命中3x3宝石牌。也就是说，在所有可能的情况中，能够成功的情况只有一种，那就是5张翻出的牌即是全部5张已有的3x3宝石牌。

数学比较好的同学，马上想到概率为：510 ×49 ×38 ×27 ×16 \frac{5}{10}\times\frac{4}{9}\times\frac{3}{8}\times\frac{2}{7}\times\frac{1}{6}

数学更比较好的同学，直接想到了概率为：1C 5 10 \frac 1{C_{10}^5}

这两种方法计算的结果都为：0.003968 0.003968

也就是说，玩1000局游戏，大概只出现3局这样的情况。如果一天玩两局，半年左右才出现一局。

应该可以说，这确实是个极小概率事件了。

作为一名搞工程应用的同学，弄到这样，差不多就行了。

而作为一名搞理论研究的同学，只满足于计算结果是不够的，还必须努力发掘现象背后的普遍规律。

特别是作为一名有优秀数学内涵和编程背景的死宅，如果不能理论和实践两方面严格论证的话，那生活多没有色彩啊。

于是，有个死宅开始了自虐（自high）之路。

这个死宅发现，游戏时常常会出现，5张牌中有一张牌功亏一篑，众人每每叹息，又差一点点。这个现象引起了死宅的思考，这个差一点点的概率又是多少。或者说，5张牌里有4张3x3宝石牌的概率是多少？

最直接的反应就是，可能是“完美”概率的5倍吧？你看，5张里缺一张。

好吧，拿前面那个概率X5看看：0.003968×5=0.01984 0.003968\times5 = 0.01984

真的是这样的吗？

你这个死宅是怎么通过高数考试的？

你以为你说什么我们就信啊？你证明给我看啊！

靠………………

那就证明一下吧。

设卡牌总数为N，符合条件的卡牌数为M，一次随机选出K张牌。计算K张牌中全部为符合条件的牌的概率。

搞个图示说明一下：

O 可以是 X ，也可以是 Y。

+—————————–+

|| O O O O … O O O O || <– N个O，M个X

+—————————–+

取出K个O。

判断是否有K个X？

（这个图示好难看啊！）

（看懂了？不管了~）

从概率论的角度来讨论这个问题。

每一次选牌都是一种排列，也是一个样本点。

所有不同排列的数量，就是所有样本点的数量。

每一个样本点的几率是相等的。

这样一来，要计算概率，只要解决以下两个问题。

样本总数有多少？

符合条件的样本有多少？

思路1：

样本总数，就是洗牌的排列总数，就是从N张牌中依次取出N张牌的排列总数。

即：P n n P_n^n

符合条件的样本数，就是前面K张牌有符合条件的排列数，乘上后面N-K张牌任意的排列数。

即：P k m ×P n−k n−k P_m^k \times P_{n-k}^{n-k}

概率为：符合条件样本数样本总数 =P k m ×P n−k n−k P n n \frac{符合条件样本数}{样本总数} = \frac{P_m^k \times P_{n-k}^{n-k}}{P_n^n}

代入N=10，M=5，K=5，验算一下概率为：0.003968 0.003968

思路2：

样本总数，就是选出的M张牌的排列总数，就是从N张牌中依次取出M张牌的排列总数。

即：P m n P_n^m

符合条件的样本数，就是前面K张牌有符合条件的排列数，乘上后面M-K张牌任意的排列数。

即：P k m ×P m−k n−k P_m^k \times P_{n-k}^{m-k}

概率为：符合条件样本数样本总数 =P k m ×P m−k n−k P m n \frac{符合条件样本数}{样本总数} = \frac{P_m^k \times P_{n-k}^{m-k}}{P_n^m}

代入N=10，M=5，K=5，验算一下概率为：0.003968 0.003968

写到这里，我们回头去检验那个4张牌为3x3宝石牌的概率。

N=10，M=5，K=4，代入计算

概率为：0.023810 0.023810

尼玛，和我的计划有出入啊！

所以说，要学好数学。真理和胡说八道之间还是有很大差距的。

（……）靠！到底谁对啊？

（……）貌似公式君好像有点道理~

……

理论是这样写完了，规律似乎也发现了，任意N，M，K貌似也能代入计算了。

但是，何以证明？你说我们就信啊！

概率这种东西怎么证明啊？

具有优秀数学内涵的死宅告诉你：有一门数学课，叫做“统计学”。

那么统计学里，概率应该怎么证明啊？很简单，我们用的方法是，抽样->调查->统计。简单来说，就是把样本“泼洒”出去，然后“捡”一些回来看看，对比一下正确/错误的比例，再和理论结果验证一下。

如果说得服你，就信我吧，说不服你，就去信上帝吧~

古典的统计学家要怎么说服大众呢？很简单，我们真的去捡样本啊。深入基层，了解群众，手拿纸本，奋笔疾书。要说计算钢珠落入圆圈中的概率，我们就真的抓一把钢珠丢下去，看有多少钢珠掉圆圈里。

不管你信不信，反正我是信了~

现代的统计学家怎么办？也去捡样本？

好吧，我们的时代已经很先进了，我们有计算机这货了，统计学家也进化成高级码农了。

那么，统计学码农应该做什么呢？

很明显好不好，埋头写代码呗！

写代码，这么难听！我们有名字的好不好。我们写的叫做“仿真”程序。

我们模仿现实发生的事情写程序，从而不用再苦逼的去鸟不拉屎的地方搞抽样，然后一字一句的写在古老发黄的记事本上，再拿个算盘或者计算器敲上一个通宵。

（……）虽然说不用再敲算盘或计算器了，但是作为一个码农，还是会在电脑键盘上敲上一个通宵的。

废话了那么多，来看看仿真程序吧。

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>

#define N 10
#define M 5
#define K 5

#define TEST_ROUND 100000

///------------------------------------------------------------------------------
double fact(int n)
{
#define FACT_MAX_NUMBER 170

if (n < 0 || n > FACT_MAX_NUMBER)
return 0;

double dRet = 1;
for (int i = n ; i > 1 ; --i)
dRet *= i;

return dRet;
}

void test_analytical()
{
double dHits = (fact(M)/fact(M-K))*fact(N-K);
double dAll = fact(N);
printf("Analytical : %f\n", dHits/dAll);
}

///------------------------------------------------------------------------------
bool is_hit1(int k, int n, int m, bool bBreak)
{
bool bHit = true;
while ( k-- )
{
int nIndex = rand() % n;

if ( nIndex + 1 > m )
{
bHit = false;

if ( bBreak )
break;
}
else
m--;

n--;
}

return bHit;
}

void test_probability1(int nSeed, bool bBreak)
{
srand(nSeed);

int nHitTimes = 0;
int nCount = TEST_ROUND;
while ( nCount-- )
{
if ( is_hit1(K, N, M, bBreak) )
nHitTimes++;
}

double dProbability = (double)nHitTimes/TEST_ROUND;
printf("Simulation 1 (%s): %f\n", bBreak ? "Break":"Loop All", dProbability);
}

///------------------------------------------------------------------------------
void get_indexs(int vIndexs[], int n, int k)
{
for ( int ii = 0; ii < n; ii++ )
vIndexs[ii] = 0;

while ( k-- )
{
int nNumeric = rand() % n;
nNumeric++;

int nCount = 0;
for ( int nPos = 0; nPos < N; nPos++ )
{
if ( 0 == vIndexs[nPos] )
{
nCount++;
if ( nCount == nNumeric )
{
vIndexs[nPos] = 1;
break;
}
}
}

n--;
}
}

bool is_hit2(int vCards[], int n, int k, bool bFirstK)
{
if ( bFirstK )
{
while ( k-- )
{
if ( 0 == vCards[k] )
return false;
}
}
else // random K position
{
int vIndexs
;
get_indexs(vIndexs, n, k);

for ( int nPos = 0; nPos < n; nPos++ )
{
if ( vIndexs[nPos] )
{
if ( 0 == vCards[nPos] )
return false;
}
}
}

return true;
}

void test_probability2(int nSeed, bool bFirstK)
{
srand(nSeed);

int vCards
;
int nHitTimes = 0;
int nCount = TEST_ROUND;
while ( nCount-- )
{
get_indexs(vCards, N, M);
if ( is_hit2(vCards, N, K, bFirstK) )
nHitTimes++;
}

double dProbability = (double)nHitTimes/TEST_ROUND;
printf("Simulation 2 (%s K): %f\n", bFirstK ? "First":"Random", dProbability);

}

///------------------------------------------------------------------------------
void main()
{
test_analytical();
test_probability1(time(NULL), true);
test_probability1(time(NULL), false);
test_probability2(time(NULL), true);
test_probability2(time(NULL), false);
}

输出结果：

N 10

M 5

K 5

Analytical : 0.003968

Simulation 1 (Break): 0.003950

Simulation 1 (Loop All): 0.003940

Simulation 2 (First K): 0.003940

Simulation 2 (Random K): 0.004080

N 10

M 5

K 4

Analytical : 0.023810

Simulation 1 (Break): 0.023850

Simulation 1 (Loop All): 0.023800

Simulation 2 (First K): 0.023780

Simulation 2 (Random K): 0.023870

N 10

M 5

K 3

Analytical : 0.083333

Simulation 1 (Break): 0.083320

Simulation 1 (Loop All): 0.083690

Simulation 2 (First K): 0.083000

Simulation 2 (Random K): 0.083910

N 10

M 4

K 4

Analytical : 0.004762

Simulation 1 (Break): 0.004860

Simulation 1 (Loop All): 0.004890

Simulation 2 (First K): 0.004890

Simulation 2 (Random K): 0.004600

N 10

M 4

K 3

Analytical : 0.033333

Simulation 1 (Break): 0.033970

Simulation 1 (Loop All): 0.032990

Simulation 2 (First K): 0.033500

Simulation 2 (Random K): 0.034150

N 10

M 4

K 2

Analytical : 0.133333

Simulation 1 (Break): 0.133760

Simulation 1 (Loop All): 0.132460

Simulation 2 (First K): 0.132750

Simulation 2 (Random K): 0.133990

可以看到仿真方法有两种。

挑选法。每一步都挑选符合条件的，遇到不符合条件的，挑选就可以结束了。这种方法适合条件比较简单，并适合每一个单步可以直接判断的。

模拟法。模仿现实环境，散布样本点，并抽样判断符合条件。这种方法仿真性比较强，适合条件比较复杂，尤其是组合条件比较多的。缺点是，程序运行速度比较慢。

从结果来看，无论N，M，K的值如何变化，理论计算概率值和仿真计算概率值都非常接近。

所以，可以有足够的理由认为：

你应该相信我们，相信科学。

当然，赞美一下上帝也不错哟~