今天在查找马尔科夫链的过程中,在网上看到一个有意思的问题,于是用python将其做了实现和改进,题目如下:
如下图所示的迷宫共有9个格子,相邻格子有门相通,9号格子就是迷宫出口. 整个迷宫将会在5分钟后坍塌.
1号格子有一只老鼠,这只老鼠以每分钟一格的速度在迷宫里乱窜(它通过各扇门的机会均等)。求此老鼠在迷宫坍塌之前逃生的概率。如果这只老鼠速度提高一倍,则老鼠在迷宫坍塌之前逃生的概率能增加多少?
题目中可以看空间1为入口,空间9为出口,思路上既可以用迭代的方法,也可以用矩阵形式来编写 。
先是基于马尔科夫链的矩阵形式。马尔科夫链的本质是一个时间、状态离散的马尔可夫过程,本质特征就是,模型当前的状态,依赖也仅依赖于前面的那一个状态。
设S为状态分布矩阵,P为状态转移矩阵,那么老鼠的每一次移动,这个过程即是乘上P的过程,当前的状态分布矩阵,即是P与S的矩阵乘积。(此处S需以列向量的形式呈现)具体不在赘述,如果编程采用矩阵运算的形式,那么对于题中这样一个长度为9的状态分布矩阵,就需要用到9*9的状态转移矩阵来计算,如果迷宫变大了,矩阵将会变得非常大。
那么从迭代的思路来看呢:
在初始状态的情况下,老鼠所在的位置,空间状态概率被记为1,其余位置的概率为零。
现在老鼠只有两个选择,要么去2,要么去4,在机会均等的情况下,那么下一个状态(也就是老鼠动一步的情况下),老鼠到达2和4的概率都是0.5,而老鼠留在1的概率为0,所以2、4的当前状态的概率都变成了0.5,而1的当前状态的概率,变成了0。
现在要进行下一次移动,又会面临新的选择,老鼠每一次的移动,就以此类推。
我们抛开确定的状态,从整体任选一个状态。如果老鼠的当前位置处于四个角上,即周围存在两个可以去的空间,那么移动机会均等的情况下,去往每一个空间的可能性为1/2。如果老鼠当前位于边上,周围就有三个可以去的空间,那么这个概率就变成了1/3,以此类推,内部的话就会有四个空间,概率是1/4。
然后我们再考虑,从上一个状态到达当前状态的概率是多少呢?那就要看,当前位置,老鼠可能是从哪里来的。如果老鼠来源的那个空间,周围有三个可以去的空间,那么现在空间的概率,就是老鼠来源空间概率的1/3,当然,当前空间不仅仅只会有这一个老鼠来源,会有很多个,那么当前空间的概率,就是这些来源的累积。
除此以外,要考虑几个点:老鼠走到出口后,是不会再动的,这就意味着,出口相邻的点,是不会接受来自出口的老鼠的,也就是从出口出来的老鼠,概率应当统统记为0。既然老鼠在出口去不了其他地方,是在出口是不会动的,那么老鼠在出口的概率,就每次更新状态时候应当保留。(出口,是一个只进不出的点)
此外,这个程序可以自行设定迷宫长与宽,设定老鼠当前位置和出口位置,设置步数。
可以的改进:可以扩展到三维或是多维;设定某个移动方向的倾向(也就是概率,目前是概率均等)
import numpy as npdef connect(point,size,rat_exit):flag = 0m,n = sizex,y = pointif [x,y] == rat_exit:return 0for [xx,yy] in [[x-1,y],[x+1,y],[x,y-1],[x,y+1]]:if 0<=xx<m and 0<=yy<n:flag += 1if flag == 1:return 1elif flag == 2:return 0.5elif flag == 3:return 1/3elif flag ==4:return 0.25#mainsize = input('棋盘的长与宽:')size = list(int(x) for x in size.split(','))m,n = sizerat = input('入口所在坐标:')rat = list(int(x) for x in rat.split(','))a,b = ratrat_exit = input('出口所在坐标:')rat_exit = list(int(x) for x in rat_exit.split(','))aa,bb = rat_exit# l = int(input('逃离步长:'))t = int(input('逃离步数:'))A = np.zeros([m,n])A[a,b] = 1print('初始状态为:')print(A)for i in range(t):AA = A.copy()for x in range(m):for y in range(n):if [x,y] != rat_exit:AA[x,y] = 0for [xx,yy] in [[x-1,y],[x+1,y],[x,y-1],[x,y+1]]:if 0<=xx<m and 0<=yy<n:AA[x,y] += A[xx,yy]*connect([xx,yy],size,rat_exit)A = AAprint('第{}次行动后的状态为:'.format(str(i+1)))print(A)print('\n')P = A[aa,bb]print('老鼠到达出口的概率为:'+ str(P))
运行结果示例: