Rを用いてビンゴの確率について

ビンゴ（bingo）する確率について調べてみたけどあまりよくわからなかった。
とりあえず5×5でフリーマスありで、10000回（コードでは人数になってるけど）やってみた結果を添付します。

細かい結果を見てみると、mizutaの結果と少し異なり、
dinとほとんど同じ結果が得られた。
つまり、mizutaさんの式がどこか間違っているのだろうか。考えれば考えるほどわからなくなった。
論文では他に、"A New Look at the Probabilities in Bingo"、Some Probability Problems Concerning the Game of Bingoなど読んだ。

コード中では3×3や、4×4などを任意に変えられるようにしました。freeマスは4×4などでは上手く行かないかも。

以下コード

players<-10000
numrows <- numcols <- 5
maxNum <- 75
isfree <- TRUE # free : center

all_bingo<-list()
numBingoMatrix <- matrix(0,nrow=maxNum,ncol=players)

countBingo <- function(bingo,numrows){
  #bingo bingo matrix
  count <- 0
  for(i in 1:numrows){
    if(sum(bingo[,i])==0){
      count <- count + 1
    }
  }
  for(j in 1:numrows){
    if(sum(bingo[j,])==0){
      count <- count + 1
    }
  }
  
  #antidiagonal
  if(sum(diag(apply(bingo,2,rev)))==0){
    count <- count + 1
  }
  if(sum(diag(bingo))==0){
    count <- count + 1
  }
  return(count)
}

for(i in 1:players){
  bingo<-matrix( (sample(1:maxNum,size = numrows*numcols,replace=FALSE)),ncol=numcols,nrow=numrows)
  if(isfree){
    bingo[(numcols+1)/2,(numrows+1)/2]<-0
  }
  all_bingo<-c(all_bingo,list(bingo))
}

skeleton <- all_bingo
calls <- sample(1:maxNum,size = maxNum,replace=FALSE)
bingoNum <- numeric(players)

for(j in 1:maxNum){
  unlisted <- unlist(all_bingo)
  unlisted[unlisted==calls[j]]<-0
  all_bingo<-relist(flesh=unlisted,skeleton = skeleton)
  
  for(k in 1:players){
    bingoNum[k] <- countBingo(all_bingo[[k]],numrows)
  }
  numBingoMatrix[j,]<-bingoNum
}

temp <- numBingoMatrix>0
numberOfHit<-apply(temp,1,sum)
plot(numberOfHit/players,xlab="calls",ylab="cumulative probability")
numberOfHit/players