Mathematical Formulation of Bingo Probability and MonteCarlo Simulation(R)

This article is poor translation of my japanese articles.

There are a few papers related probability of Bingo games.

"A New Look at the Probabilities in Bingo"Some Probability Problems Concerning the Game of Bingo

However, I cound not find a mathematical formulation of Bingo probability.

At first, I suppose there are  n \times n halls, max number  m (often 75), and number of calls  k.

For example, if we have  5 \times 5 halls, there are 12 lines to bingo. I suppose probability of bingo in a line  p_{1},

 p_{1}=\frac{\dbinom{n}{n} \times \dbinom{m-n}{k-n}}{\dbinom{m}{k}}=\frac{\dbinom{m-n}{k-n}}{\dbinom{m}{k}}

when there is no "FREE" hall.
This probability formulation is related to hypergeometric distribution.

Bingo happens in at least one of 12 lines, bingo probability  p is,
 p = 1-(1-p_{1})^{2 n + 2}
.

In case there is free hall,
 p = 1-(1-p_{1})^{2 n - 2}(1-p_{2})^4
where
 p_{2}=\frac{\dbinom{m-n+1}{k+1-n}}{\dbinom{m}{k}}
.





Figure shows a comparison MonteCarlo simulation and derived mathematical formulation, dotted points are simulation value and red line is derived one.

MonteCarlo simulation in 50000 times,  5 \times 5 halls , and no "FREE" hall.

f:id:saikeisai:20170131013213p:plain


R simulation code below.

players<-50000
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)


bingoProbability <- function(calls,maxNum,numrows,isfree){
p1 <- choose(maxNum-numrows,calls-numrows)/choose(maxNum,calls)
p2 <- p1
if(isfree){
p2 <- choose(maxNum-numrows+1,calls+1-numrows)/choose(maxNum,calls)
}
p <- 1-((1-p1)^(numrows*2-2))*((1-p2)^4)
return(p)
}


plot(numberOfHit/players,xlab="calls",ylab="cumulative probability")
lines(bingoProbability(calls=1:75,maxNum = 75,numrows=5,isfree=isfree),col=2)