马尔可夫过程和马尔可夫链_R：蛇形和梯形马尔可夫链-CSDN博客

马尔可夫过程和马尔可夫链

几天前，我读了一篇非常酷的博客文章，解释了如何在蛇和梯子游戏中使用马尔可夫链来建模可能的状态转换，这是我什至没有想到的马尔可夫链的用法！

尽管该示例对于理解该概念非常有帮助，但我对代码的理解是，它假定一个掷骰子大于100的掷骰都是成功的假设。

在我所知道的游戏版本中，您必须准确地赢得100。例如，如果您在98号方格上掷6，则可以向前移动2个空格至100，然后向后反弹4个空格至96。

我认为调整代码以适应此问题将是一个不错的练习：

n=100
 
# We have 6 extra columns because we want to represent throwing of the dice which results in a final square > 100
M=matrix(0,n+1,n+1+6)
rownames(M)=0:n
colnames(M)=0:(n+6)
 
# set probabilities of landing on each square assuming that there aren't any snakes or ladders
for(i in 1:6){
  diag(M[,(i+1):(i+1+n)])=1/6
}
 
# account for 'bounce back' if a dice roll leads to a final score > 100
for(i in 96:100) {
  for(c in 102:107) {
    idx = 101 - (c - 101)  
    M[i, idx] = M[i, idx] + M[i, c]
  }  
}

我们可以检查最后几行，以检查转换矩阵是否正确：

> M[95:100,95:101]
 
   94        95        96        97        98        99       100
94  0 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667
95  0 0.0000000 0.1666667 0.1666667 0.1666667 0.3333333 0.1666667
96  0 0.0000000 0.0000000 0.1666667 0.3333333 0.3333333 0.1666667
97  0 0.0000000 0.0000000 0.1666667 0.3333333 0.3333333 0.1666667
98  0 0.0000000 0.1666667 0.1666667 0.1666667 0.3333333 0.1666667
99  0 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667 0.1666667

如果我们在第99平方（最后一行）上，我们可以掷1并以100结束，以2并以99（前进1，后退1）结束，以3并以98（前进1， 2后卫），4和97（后卫1，前卫3），5和96（前卫1，后卫4）或6和95（1前卫，后卫5）。即我们可以以1/6的概率降落在95、96、97、98、99或100上。

如果我们在第96平方（第3行）上，我们可以掷1并以97结束，2并以98结束，3并以99结束，4并以100结束，5并以99（前进4，后退1）或6结束，以98（前进4，后退2）结束。即我们可以以1/6的概率降落到97，以2/6的概率降落到98，以2/6的概率降落在99或以1/6的概率降落在100。

我们可以对其他平方做类似的分析，但似乎正确计算了概率。

接下来，我们可以使用蛇和梯子更新矩阵。该代码保持不变：

# get rid of the extra columns, we don't need them anymore
M=M[,1:(n+1)]
 
# add in the snakes and ladders
starting = c(4,9,17,20,28,40,51,54,62,64,63,71,93,95,92)
ending   = c(14,31,7,38,84,59,67,34,19,60,81,91,73,75,78)
 
for(i in 1:length(starting)) {
  # Retrieve current probabilities of landing on the starting square
  v=M[,starting[i]+1]  
  ind=which(v>0)
 
  # Set no probability of falling on the starting squares
  M[ind,starting[i]+1]=0
 
  # Move all existing probabilities to the ending squares
  M[ind,ending[i]+1]=M[ind,ending[i]+1]+v[ind]
}

我们还可以简化powermat函数，该函数用于模拟一定数量的掷骰后板子的外观：

# original
powermat=function(P,h){
  Ph=P
  if(h>1) {
    for(k in 2:h) {
      Ph=Ph%*%P
    }
  }
  return(Ph)
}
 
#new 
library(expm)
powermat = function(P,h) {
  return (P %^% h)
}

initial=c(1,rep(0,n))
h = 1
> (initial%*%powermat(M,h))[1:15]
     0         1         2         3 4         5         6 7 8 9 10 11 12 13        14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45
[1,] 0 0.1666667 0.1666667 0.1666667 0 0.1666667 0.1666667 0 0 0  0  0  0  0 0.1666667  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0
     46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100
[1,]  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0  0   0

我注意到的一件有趣的事情是，与不需要完全得分100来获胜的情况相比，现在完成游戏平均需要更多的回合：

> sum(1 - game)
[1] 999

distrib=initial%*%M
game=rep(NA,1000)
for(h in 1:length(game)){
game[h]=distrib[n+1]
distrib=distrib%*%M}
plot(1-game[1:200],type="l",lwd=2,col="red",
ylab="Probability to be still playing")

我希望花更长的时间才能完成游戏，但是不会太久！我想我可能犯了一个错误，但我不确定在哪里…

更新资料

安东尼奥斯发现了我犯的错误–在第100平方时，到达100平方的概率应该为1。即我们需要像这样更新M：

M[101,101] = 1

现在，如果我们可视化他仍在玩游戏的可能性，我们将获得更准确的曲线：

distrib=initial%*%M
game=rep(NA,1000)
for(h in 1:length(game)){
game[h]=distrib[n+1]
distrib=distrib%*%M}
plot(1-game[1:200],type="l",lwd=2,col="red",
ylab="Probability to be still playing")