084. 大富翁概率(Monopoly odds)

大富翁游戏的标准棋盘大致如下图所示:

img

玩家从标记有\(GO\)的方格出发,掷两个六面的骰子并将点数和相加,作为本轮他们前进的步数。如果没有其它规则,那么落在每一格上的概率应该是\(2.5\%\)。但是,由于有方格 \(G2J\)(入狱)、\(CC\)(宝箱卡)和\(CH\)(机会卡)的存在,这个分布会有所改变。除了落在\(G2J\)方格上,以及在\(CC\)\(CH\)上抽到入狱卡以外,如果玩家连续三次都掷出两个相同的点数,则在第三次时也会直接入狱。

游戏开始时,\(CC\)\(CH\)对应的卡片将被洗牌打乱。当一个玩家落在\(CC\)\(CH\)上时,他们从宝箱卡和机会卡的牌堆最上方取一张卡并遵循指令行事,并将该卡再放回牌堆的最下方。宝箱卡和机会卡都各有\(16\)张,但我们只关心会影响到移动的卡片,其它的卡片我们都将无视它们的效果:

  • 宝箱卡(2/16张卡)
  • 回到起点\(GO\)
  • 进入监狱\(Jail\)
  • 机会卡(10/16张卡)
  • 回到起点\(GO\)
  • 进入监狱\(Jail\)
  • 移动到\(C1\)
  • 移动到\(E3\)
  • 移动到\(H2\)
  • 移动到\(R1\)
  • 移动到下一个\(R\)(铁路公司)
  • 移动到下一个\(R\)
  • 移动到下一个\(U\)
  • 后退三个方格

这道题主要考察掷出骰子后停在某个特定方格上的概率。显然,除了停在\(G2J\)上的可能性为\(0\)之外,停在\(CH\)格的可能性最小,因为有\(5/8\)的情况下玩家会移动到另一格。我们不区分是被送进监狱还是恰好落在监狱\(Jail\)这一格,而且不考虑需要掷出两个相同的点数才能出狱的要求,而是假定进入监狱的第二轮就会自动出狱。

从起点\(GO\)出发,并将方格依次标记\(00\)\(39\),我们可以将这些两位数连接起来表示方格的序列。从统计学角度,三个最有可能停下的方格分别是\(Jail(6.24\%)\)或方格\(10\),接下来是\(E3(3.18\%)\)或方格\(24\)以及\(GO(3.09\%)\)或方格\(00\)。这三个方格可以用一个六位数字串表示:\(102400\)。如果我们不用两个六面的骰子而是用两个四面的骰子,求出三个最有可能停下的方格构成的数字串。

分析:解决这道题的思路大概可以分为两种,一种是使用蒙特卡洛模拟,根据游戏的设定去模拟花游戏的过程,经过足够多轮的游戏,统计落到各个方格的概率,考虑题目并不要求给出各个方格具体的概率值,而只要概率的排序准确就可以了,所以模拟可能造成的误差并不会有太大影响,只要在写代码时能够充分准确的表达游戏设定,一般都能得到正确的结果。

另一种方法则是根据概率论知识直接计算出落在每个方格的概率,这需要用到一种叫做马尔克夫链的数学工具,我虽然在之前听说过很多次马尔科夫链,但从来没有深入了解过,所以为了加深对这种数学工具的理解,我选择了第二种方法。解决这道题大约花了我一周时间,其中一两天用来学习马尔科夫链,剩下的时间都在做无尽的debug工作,事实证明对大富翁游戏的熟悉还是很重要的,我犯的很多错误都是因为对游戏设定理解有误造成的。

接下来我会简要介绍一下马尔科夫链的基础知识,这里只会提到一些和题目相关的结论,至于证明与推导的过程大家可以参见相关的概率论书籍。在此基础上,我会介绍如何用马尔科夫链来分析大富翁游戏,如何计算出落在各个方格的概率。最后我会讲一下编程实现,尤其是如何生成正确的转移概率矩阵。

一、马尔科夫链简介

马尔科夫链是一种随机过程,和高斯过程、泊松过程等其它随机过程的主要区别在于它具有马尔科夫性质。所谓马尔科夫性质是指系统未来时刻的概率分布只与当前时刻有关,而与过去所有时刻都无关,这种性质也称作无记忆性。从数学定义的角度说,假设我们的序列状态是\(X_1,X_2,X_3,\cdots,X_n\),则:

$$ \Pr(X_{{n+1}}=x\mid X_{1}=x_{1},X_{2}=x_{2},\ldots ,X_{n}=x_{n})=\Pr(X_{{n+1}}=x\mid X_{n}=x_{n}) $$

如果一个随机过程具有马尔科夫性质,我们就可以更加容易的对其进行处理。特别的,如果状态之间的转移概率保持不变,即:

$$ \Pr(X_{{n+1}}=x\mid X_{n}=y)=\Pr(X_{n}=x\mid X_{{n-1}}=y)\, $$

这种马尔科夫链我们通常称之为时齐马尔可夫链与静态马尔科夫链,如果同时具有以上两个性质,既然某一时刻状态转移的概率只依赖于它的前一个状态,那么我们只要能求出系统中任意两个状态之间的转移概率,这个马尔科夫链的模型就确定了。幸运的是,我们在实际中经常接到的马尔科夫链都是静态的。因此我们在已知初始概率分布\(\pi_0\)和转移概率\(p\)情况下,就可以推导出任意有限步之后的概率分布\(\pi_n\)。我们看一个例子:

img

上图通过有向图的形式展示了一个与股市有关的马尔科夫链,股市有三个状态,分别是牛市、熊市和横盘。图中的箭头和旁边的数字表示从一个状态转移到另一个状态的概率,比如我们可以看到从牛市转移到熊市的概率为\(0.075\),转移到横盘的概率为\(0.025\),而保持在牛市状态的概率是\(0.9\),可以看到三个概率加起来正好为一。因此假设现在处于牛市的概率为0.5,熊市的概率为0.25,横盘的概率为0.25,则初始的概率分布可以表示为\(\pi_0=(0.5,0.25,0.25)\)。我们把牛市、熊市以及横盘状态分别标记为,则如果我们想知道下一时刻为牛市的概率,根据条件概率与全概率公式,则有:

$$ p\left( x_1=0 \right) =\sum_{i=0}^2{p\left( x_0=i \right) p\left( x_1=0|x_0=i \right)}=0.5\times 0.9+0.25\times 0.15+0.25\times 0.25=0.55 $$

同理我们可以得到下一时刻为熊市的概率:

$$ p\left( x_1=1 \right) =\sum_{i=0}^2{p\left( x_0=i \right) p\left( x_1=1|x_0=i \right)}=0.5\times 0.075+0.25\times 0.8+0.25\times 0.25=0.3 $$

以及下一时刻为横盘的概率:

$$ p\left( x_1=2 \right) =\sum_{i=0}^2{p\left( x_0=i \right) p\left( x_1=2|x_0=i \right)}=0.5\times 0.025+0.25\times 0.05+0.25\times 0.5=0.15 $$

如果我们把初始分布表示为一个行向量\(\pi_0\),同时把转移概率表示为一个矩阵,矩阵中元素\(p_{ij}\)表示当前状态为\(i\)且下一个状态为\(j\)的条件概率,即\(p_{ij}=p(x_{n+1}=j|x_n=i)\),则上述股市状态的转移矩阵为:

$$ P=\left( \begin{array}{ccc} 0.9&0.075&0.025 \\ 0.15&0.8& 0.05 \\ 0.25&0.25&0.5 \end{array} \right) $$

则我们只需要用表示初始分布的行向量右乘这个矩阵,就可以得到下一个状态的分布:

$$ \pi _1=\pi _0P=\left( 0.5,0.25,0.25 \right) \cdot \left( \begin{matrix} 0.9& 0.075& 0.025\\ 0.15& 0.8& 0.05\\ 0.25& 0.25& 0.5\\ \end{matrix} \right) =\left( 0.55,0.3,0.15 \right) $$

可以看到这里的计算结果与上面我们分别计算的结果是一致的。如果我们想知道再下一状态的概率分布\(\pi_2\),只需要再右乘一个转移矩阵就可以了,也即\(\pi_2=\pi_1P=\pi_0P^2\)。同理我们可以得到第\(n\)步之后的概率分布\(\pi_n=\pi_0P^n\),因此我们只需要计算出转移矩阵的乘方就可以计算出任意步之后的概率分布。

进一步来看,对于满足某些性质的转移矩阵\(P\),当步数趋于无穷时,其分布会逐渐趋于稳定,最终会收敛到一个极限分布,假设记其为\(\pi^*\),即\(\lim _{{n\rightarrow \infty }}\pi{\mathbf {P}}^{n}=\pi ^{*}\)。需要注意的是,极限分布的定义与初始分布无关,即对任意的初始分布,当时间步趋于无穷时,随机变量的概率分布都趋于极限分布。根据极限分布的定义易知\(\pi^*P=\pi^*\),也即极限分布在转移矩阵的作用下不会发生变化。对于上面介绍的股市的例子,下表展示了它三十步内的概率分布,可以发现分布在第\(24\)步以后就基本维持不变了,这个分布就是极限分布,同时也是稳定分布。

img

通过上面的例子我们可以发现,可以通过在初始分布上不断右乘转移矩阵来得到极限分布,如果发现概率分布向量不再变化,那么这个分布就是极限分布。除此之处,我们也可以通过极限分布的性质\(\pi^*P=\pi^*\)发现\(\pi^*\)为矩阵\(P\)的特征值为一左特征向量,因此我们只需要求出这个特征向量就知道极限分布了。实际中我们可以列方程组求解,假设\(\pi^*=(x_1,x_2,x_3)\),则:

$$ \pi ^*P=\left( x_1,x_2,x_3 \right) \cdot \left( \begin{matrix} 0.9& 0.075& 0.025\\ 0.15& 0.8& 0.05\\ 0.25& 0.25& 0.5\\ \end{matrix} \right) \Rightarrow \begin{cases} 0.9x_1+0.15x_2+0.25x_3=x_1\\ 0.075x_1+0.8x_2+0.25x_3=x_2\\ 0.025x_1+0.05x_2+0.5x_3=x_3\\ \end{cases} $$

解以上方程组,我们可也可以得到完全一致的极限分布。

二、使用马尔科夫链分析大富翁游戏

题目要求最有可能停下的三个方格,实际上就是要求在游戏轮数趋于无穷时落在各个方格概率的极限分布,因此我们可以用马尔科夫链来进行分析。在此之前,我们需要对大富翁的游戏设定做一下梳理。

有两个设定为大富翁游戏引入了随机因素,一是通过抛掷两个骰子并将其点数相加决定移动的方格数;二是当移动到宝箱卡和机会卡方格后,通过抽卡来决定移动的方格。这两个设定显然是相互独立的,则我们可以分别考虑两个事件发生的概率。我们还需要对大富翁的部分游戏规则做一些简化,这些规则要不违反了马尔科夫性质,要不会给我们的分析带来不必要的麻烦,主要如下:

在以上假设前提下,我们可以开始计算大富翁游戏的转移矩阵了。记这个矩阵为\(T\),则其为一个\(40\times40\)的方阵,其中每个元素\(t_{ij}\)表示从方格\(i\)移动到方格\(j\)的概率。前面已经提到移动的方格数量受到两个随机事件的影响,分别是掷骰子的点数和以及抽到的宝箱与机会卡。因为每一轮都是先掷骰子后抽卡,且两个事件为相互独立事件,则设\(d_{ij}\)表示因为掷骰子从\(i\)移动\(j\)的概率,\(c_{ij}\)表示因为抽卡从\(i\)移动\(j\)的概率,则有:

$$ t_{ij}=\sum_{k=0}^{39}{d_{ik}\cdot c_{kj}} $$

上式的意思是从方格\(i\)到方格\(j\)的转移概率等于因为掷骰子从方格\(i\)到所有方格的概率乘以因为抽卡从所有方格到方格\(j\)的概率之和。如果我们把因为掷骰子从\(i\)移动\(j\)的概率表示成矩阵\(D\),把因为抽卡从\(i\)移动\(j\)的概率记为矩阵\(C\),两个矩阵都是\(40\times40\)的方阵,则转移矩阵为两个矩阵的乘积,即\(T=D\cdot C\)

在得到转移矩阵后,我们只需要使用上一节介绍的方法计算出各个方格概率的极限分布,对其从大到小排序,找到排名前三概率对应的方格编号,就可以给出题目的答案了。下面我们具体来看一看如何构造矩阵以及计算极限分布。

三、构造矩阵

首先来看一下如何构造矩阵\(D\),该矩阵中每个元素\(d_{ij}\)表示因为掷骰子从方格\(i\)移动到方格\(j\)的概率。因为假设四,我们只需要考虑一次掷骰子中两个骰子点数之和的概率即可。对于正六面体和正四面体的骰子,其点数之和的分布概率如下表:

img

如对于正四面体骰子,两个骰子的点数组合共有\(4^2=16\)种,其中点数之和为二的有一种,为三的有两种,依次类推,而出现特定点数之和的概率等于它的种数除以总的组合数,如点数之和为五的概率等于\(4/16=0.25\)。对于正六面体骰子也是一样的,只不过其组合数更多,同时各个点数之和的分布更加均匀。

知道骰子点数之和的概率分布后,就可以构造矩阵了。比如\(d_{02}\)表示从方格0移动到方格2的概率,也就是对应的点数之和为2,所以对于正四面体,这个概率应该为\(6.25\%\);对于正六面体,这个概率为\(2.78\%\)。与之类似,\(d_{36}\)应为从方格3移动到方格6总共移动3格的概率,所以应该是\(12.5\%\)\(5.56\%\)。需要注意的是,当出发方格编号加上移动的方格数量大于四十的时候,最终到达的方格编号应该等于移动的方格数量除以四十的余数。所以从方格\(i\)移动到方格\(j\)的概率\(d_{ij}=p(scores=(j-i)\ mod\ 40)\)。使用这个公式以及\(numpy\)模块,我们就可以构造出矩阵\(D\)了,代码如下:

import numpy as np

def get_dice_mat():
    n,s = 40,4
    freq = {2:1,3:2,4:3,5:4,6:3,7:2,8:1}
    prob = {k:v/s**2 for k,v in freq.items()}
    mat = np.zeros((n,n))
    for i in range(n):
        for j in range(n):
            mat[i,j] = prob.get((j-i)%40,0)
    return mat

下一步是构造矩阵\(C\),其中每个元素\(c_{ij}\)表示因为抽卡从方格\(i\)移动到方格\(j\)的概率。需要注意的 是,四十个方格中有三个机会卡方格,三个宝箱卡方格以及一个入狱方格,虽然入狱方格不涉及抽卡,但他们都会改变玩家移动的路径,所以这里我们放在一起来考虑:入狱方格比较简单,它的编号是30,一旦落在这个方格就马上入狱,也就是说进入方格10,所以我们只需要设定\(c_{30,10}=1\)就可以了。

接下来是宝箱卡方格,编号分别为2,17和33,在总共16张卡中有2张涉及到方格移动,分别是回到起点以及入狱,根据假设一,抽到这两张的卡概率均为\(2/16/2=1/16\)。因为16张卡中只有2张涉及到方格移动,剩下14张卡只会让玩家停在原地,所以玩家留在原地的概率应该为\(14/16=7/8\)。因此,我们可以推出\(c_{2,0}=0.0625,c_{2,10}=0.0625,c_{2,2}=0.875\),另外两张宝箱卡也是一样的。

机会卡方格是最复杂的,编号分别为7,22和36,总共16张卡中有10张涉及到方格移动,因此每个移动方向的概率应该等于\(10/16/10=1/16\),而停留在原地的概率为\(6/16=3/8\)。这十张卡片中前六张与所处的具体方格无关,直接移动到相应方格就可以了,而后四张与所处方格位置有关。比如当处于方格7时,下一个火车站为方格15,而处于方格22时,下一个火车站则为方格25,因此对这四张卡片我们应该根据所处方格来处理。结果如下:

img

需要注意的是,对于方格36,它的下一个火车站也是R1,也就是方格5,所以方格5的概率实际上是\(0.0625+0.125=0.1875\)。另外对于这七个方格以外的其它方格,并不会因为抽卡而移动到别的方格,所以他们留在当前方格的概率应该为一,也即:

$$ c_{ii}=1,for\,\,i\,\,\ne \,\,2,7,17,22,30,33,36 $$

基于以上逻辑,我们可以通过修改矩阵相应元素的方法来构造矩阵\(C\),代码如下:

def get_card_mat():
    n = 40
    card_mat = np.zeros((n,n))
    for i in range(n):
        if i in [2,17,33]:
            card_mat[i,0] = 0.0625
            card_mat[i,i] = 0.875
            card_mat[i,10] = 0.0625
        elif i == 30:
            card_mat[i,10] = 1
        elif i == 7:
            card_mat[i,i] = 0.375
            card_mat[i,15] = 0.125
            for j in [0,4,5,10,11,12,24,39]:
                card_mat[i,j] = 0.0625
        elif i == 22:
            card_mat[i,i] = 0.375
            card_mat[i,25] = 0.125
            for j in [0,5,10,11,19,24,28,39]:
                card_mat[i,j] = 0.0625  
        elif i == 36:
            card_mat[i,i] = 0.375
            card_mat[i,5] = 0.1875
            for j in [0,10,11,12,24,33,39]:
                card_mat[i,j] = 0.0625
        else:
            card_mat[i,i] = 1.0
    return card_mat 

通过以上方式,我们就把矩阵\(D\)与矩阵\(C\)构造完了,然后我们用矩阵\(D\)乘以矩阵\(C\)就可以得到转移矩阵\(T\)。对于这三个矩阵,我们需要检查一下它们每一列的和是否为一,如果不是则说明构造过程中有错误,需要再回头检查一下,直到每一列的和都为一为止。

在得到正确的转移矩阵之后,我们可以使用第一节介绍的几种方法来计算它的极限分布,我使用的是不断右乘转移矩阵的方法,发现经过两百多次有迭代,分布向量不再发生变化,这个分布就是极限分布了:

def get_steady_vector():
    C = get_card_mat()
    D = get_dice_mat()
    T = D @ C
    vec_i = np.array([1]+[0]*39)                                 
    while True:
        vec_s = vec_i @ T
        if np.array_equal(vec_s,vec_i):
            return vec_s
        else:
            vec_i = vec_s    

最后结果如下:

img

可以发现,对于正六面体骰子,概率前三的方格编号分别为10,24和0,和题目中给出的顺序一致。而对于正四面体的骰子,概率前三的方格编号分别为10,15和24。

四、一点延伸

根据欧拉计划网站这道题的讨论帖子,题目的灵感应该是来源于Ian Stewart在1996在《科学美国人》杂志上发表的文章,这篇文章从落在每个方格概率的角度讨论了大富翁游戏是否公平的问题,因为作者对在大富翁游戏的设定做了很多不符合现实的简化,所以得出结论落在每个的概率都是\(2.5\%\),所以这个游戏是公平的。文章发表后激起了很多讨论,很多人给作者发信表示在考虑大富翁游戏中各项现实设定后,落在每个方格的概率实际上会有所差异,因此作者在几个月后又发表了一篇文章,介绍了如何使用马尔科夫链的方法来分析大富翁游戏,并考虑了各项现实游戏中的规则,最终得到的落在各个方格的概率和我在上表中和大家展示的结果一致,概率排名前三的方格分别为10,24和0,对应的概率分别为5.89%,3.19%和3.11%。根据这个结果,作者显然没有考虑连续掷出相同点数骰子的规则,在考虑这个规则后,得出的结果应该和题目给出的结果是一致的。但因为连续两次和三次掷出相同点数骰子的概率实际上很小,所以对最后的概率数值影响不大,对概率的排序没有影响。

实际上,网上讨论如何用马尔科夫链分析大富翁游戏的文章很多,但是基本上都只介绍了原理,而对具体如何来构造和计算转移矩阵语焉不详,所以当自己具体来分析时就不知道怎么做了。为了展示这些矩阵以及具体的计算过程,我做了一个EXCEL表格,大家可以从这个表格中看出具体的构造与计算过程是怎么样的。

作为一个对大富翁游戏和马尔科夫链均不熟悉的门外汉,在结束这篇冗长的文章之后,也许应该尝试玩一下这个游戏,希望这里的分析不会毁掉我们玩这个游戏的乐趣。

完整代码参见:https://github.com/sorrowise/euler/blob/master/084.py