Plain vanilla BlackJack simulation with R

[This article was first published on The Beginner Programmer, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Please before continue reading, make sure to read the disclaimer at the bottom of this article.

Here is a simulation I run with R in the same period I created the poker one.

I have just decided to call it plain vanilla since neither double down or split pairs are allowed. Rules are as basic as they can be.

The code looks like messy, I know, If I had to do it now, I guess I would do many things differently, however, the results looks fine, and the code runs fine as well, therefore it is not entirely to be thrown away I believe.

The simulation is divided in two parts. In the first one I looked for the probabilities for each possible event (that’s to say: “player wins”, “tie”, “dealer wins”). Since the code looks pretty hairy, many explanations are provided.

A small but important note: if you want to run the simulation, since it is a bit demanding computationally speaking, in order not to crash RStudio you may want to run a line a time and not the entire code all at once, or just maybe use only some functions, just be careful because since we run the simulation a big number of times and the functions are not the fastest to be run, your computer might complain a little.

Here is the code:

# In this example we simulate a plain vanilla Blackjack game
# and we calculate the probabilites of win and tie for both the player and the dealer
# A deck of 52 cards (Jokers are excluded), A is an ace.
# Please note: as in BlackJack figures worth 10 each, we will deal with each picture as if it were a ten
deck = c("A","A","A","A",10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,10,9,9,9,9,8,8,8,8,7,7,7,7,6,6,6,6,5,5,5,5,4,4,4,4,3,3,3,3,2,2,2,2)
# Usually, in casinos, Blackjack games are played
# with 4 to 6 decks. The following function will take in
# account this matter and return the appropriate deck
# according to the number the casino uses (n)
deck_in_use <- function(n)
{
deck_to_return = rep(deck,n)
return(deck_to_return)
}
# The function hand takes a deck and returns a hand of Blackjack
hand <- function(deck_to_use, player_hand)
{
# We fix player's hand
p_hand = player_hand
# We fix the deck
deck = deck_to_use
# Player's hand gets removed from the deck,
# since it is no longer available to be dealt
for(i in p_hand)
{
if(is.element(i,deck))
{
deck<-deck[-match(i,deck)]
}
}
# In case player's hand contains aces, the following
# code trasforms aces in 11 or 1 according to BJ rules
if((p_hand[1]) == "A" &amp;&amp; (p_hand[2] == "A"))
{
p_hand[1] <- 1
p_hand[2] <- 11
}else if((p_hand[1] == "A") && (p_hand[2] != "A"))
{
p_hand[1] <- 11
}else if((p_hand[2] == "A") &amp;&amp; (p_hand[1] != "A"))
{
p_hand[2] <- 11
}
p_hand <- as.numeric(p_hand)
# Dealer's drawing function: if dealer's hand is
# lower than 16, then another card is dealt
d_draw <- function(hand)
{
if(sum(hand) > 16)
{
return(hand)
}else
{
new_card = sample(deck,1,F)
hand = append(hand,new_card)
if(is.element(new_card, deck))
{
deck<-deck[-match(new_card, deck)]
}
# In case the newly dealt card is an ace,
# the following code deals with it
if(new_card == "A")
{
place_holder_ace_d2 <- match(new_card, hand)
vect_partial_sum_hand = hand[-place_holder_ace_d2]
vect_partial_sum_hand <- as.numeric(vect_partial_sum_hand)
if(sum(vect_partial_sum_hand) <= 10)
{
hand[place_holder_ace_d2] <- 11
}else
{
hand[place_holder_ace_d2] <- 1
}
}
hand <- as.numeric(hand)
# In case something went missing from before,
# the following code checks if every ace has
# the appropriate value
hand_minus_new_card &lt;- hand[-match(new_card, deck)]
if((sum(hand) <= 21) && (is.element(11,hand_minus_new_card)))
{
hand[match(11, hand)] <- 1
}
hand <- as.numeric(hand)
# Redraw is needed in case dealer's hand is lower than 16
return(d_draw(hand))
}
}
# Dealer's hand is eventually sampled and implemented with the
# d_draw function to simulate a real draw.
dealer_hand <- function()
{
hand = sample(deck,2,F)
for(i in hand)
{
if(is.element(i,deck))
{
deck<-deck[-match(i,deck)]
}
}
# If the dealer has two Aces they are converted into a 11 and a 1
if((hand[1] == "A") && (hand[2] == "A"))
{
hand[1] <- 11
hand[2] <- 1
hand <- as.numeric(hand)
}else
{
# If an A is in the dealer hand, it is converted into an 11
for(i in hand)
{
place_holder_ace_d <- match(i, hand)
if(i == "A")
{
hand[place_holder_ace_d] <- 11
}
}
hand <- as.numeric(hand)
}
hand <- d_draw(hand)
return(hand)
}
# Results are eventually collected into a list and returned by the function
results = list(p_hand, dealer_hand())
return(results)
}
# The following function simulates n Blackjack hands and store the
# results in a matrix. It returns some probabilities (please check below)
# n is the number of hands to be played in the simulation
simulate_games<;- function(deck_to_use, player_hand, n)
{
names_ <- c("Dealer","Player", "Test_win", "Test_draw")
result_matrix <- matrix(ncol = 4)
colnames(result_matrix) <- names_
for(i in 1:n)
{
game <- hand(deck_to_use,player_hand)
vec_result <- c(sum(game[[2]]), sum(game[[1]]))
draw <- (vec_result[2] == vec_result[1])
# Tie
if(draw)
{
vec_result[3] <- 0
vec_result[4] <- 1
result_matrix <- rbind(result_matrix, vec_result)
# Dealer hits Blackjack
}else if((vec_result[1] == 21) && (length(game[[1]]) == 2) &&; (vec_result[2] != 21) )
{
vec_result[3] <- 0
vec_result[4] <- 0
result_matrix <- rbind(result_matrix, vec_result)
# Player wins on points
}else if((vec_result[2] < vec_result[1]) | (vec_result[1] < 21))
{
vec_result[3] <- 1
vec_result[4] <- 0
result_matrix <- rbind(result_matrix, vec_result)
# Dealer wins
}else
{
vec_result[3] <- 0
vec_result[4] <- 0
result_matrix <- rbind(result_matrix, vec_result)
}
}
# Delete first row of the results matrix
result_matrix <- result_matrix[-1,]
# Show results matrix
#View(result_matrix)
probab_win = sum(result_matrix[,3])/n
probab_tie = sum(result_matrix[,4])/n
probab_dealer_win = 1 - probab_win - probab_tie
#name_data <- c("Prob of win", "Prob of a tie", "Prob of dealer win")
vect_data <- c(probab_win, probab_tie, probab_dealer_win)
#names(vect_data) <- name_data
return(vect_data)
}
# Simulation with 3 decks a hand equal to 19 and 100 games
# Sidenote: this game is as if the dealer was shuffling the
# cards at each game
simulate_games(deck_in_use(3), c(10,9), 100)
# Simulation of n_games*n_hand_game with n_deck decks and a hand equal to hand
simulate_n_games &lt;- function(n_deck, hand, n_hand_game, n_games)
{
res_matrix <- matrix(ncol = 3)
names = c("prob of win", "prob of tie", "prob of dealer win")
colnames(res_matrix) &lt;- names
for(i in 1:n_games)
{
data_vec <- simulate_games(deck_in_use(n_deck), hand, n_hand_game)
res_matrix <- rbind(res_matrix, data_vec)
}
res_matrix <- res_matrix[-1,]
View(res_matrix)
probab_win <- sum(res_matrix[,1])/nrow(res_matrix)
probab_tie <- sum(res_matrix[,2])/nrow(res_matrix)
probab_dealer_win &lt;- 1 - probab_win - probab_tie
probabilities <- c(probab_win, probab_tie, probab_dealer_win)
name_data <- c("Prob of win", "Prob of a tie", "Prob of dealer win")
names(probabilities) <- name_data
return(probabilities)
}
# Simulate 2000*100 Blackjack games with 3 decks and a hand equal to 18
simulate_n_games(3,c(10,8),2000,100)
# You may want to run the following lines one a time since they
#are a little "heavy" to compute.
# The simulation is run for hands from 16 to 21
first_row_16 <- simulate_n_games(3,c(10,6),2000,100)
second_row_17 <- simulate_n_games(3,c(10,7),2000,100)
third_row_18 <- simulate_n_games(3,c(10,8),2000,100)
fourth_row_19 <- simulate_n_games(3,c(10,9),2000,100)
fifth_row_20 <- simulate_n_games(3,c(10,10),2000,100)
sixth_row_21 <- simulate_n_games(3,c(10,11),2000,100)
BJ_probabilities <- matrix(ncol = 3)
BJ_probabilities <- rbind(BJ_probabilities,first_row_16)
BJ_probabilities <- rbind(BJ_probabilities,second_row_17)
BJ_probabilities <- rbind(BJ_probabilities,third_row_18)
BJ_probabilities <- rbind(BJ_probabilities,fourth_row_19)
BJ_probabilities <- rbind(BJ_probabilities,fifth_row_20)
BJ_probabilities <- rbind(BJ_probabilities,sixth_row_21)
BJ_probabilities <- BJ_probabilities[-1,]
View(BJ_probabilities)
hands_points <- c(16,17,18,19,20,21)
plot(hands_points, BJ_probabilities[,1], xlab = "Hands full points", main = "Probability of player win", type = "b", col="blue")
plot(hands_points, BJ_probabilities[,2], xlab = "Hands full points", main = "Probability of tie", type = "b", col="blue")
plot(hands_points, BJ_probabilities[,3], xlab = "Hands full points", main = "Probability of dealer win", type = "b", col="blue")

Let’s have a look at the results that we have got:


im1


By simply standing with a hand equal to 19, you can expect to win 59% of the times, and loose 28%. A tie is likely to occur 13% of the times.


Let’s now simulate 2000*100 games with a hand equal to 18 and 3 decks and then average the probabilities:


im5


im4


In this case the odds seem to be a lot worse than before, by just subtracting 1 point to the player’s hand we raised the probability of win for the dealer from 28% to 43.4%.


Why not run the simulation for each hand from 16 to 21 and then plot the results? Warning: your computer may yell at you after this demanding step!


im6


As we could expect, as the player’s hand gets higher, the probabilities of winning increases. For some reason which I ignore, R did not print in plain the probabilities of winning for the dealer, however you can easily get them from the table above or by calculating for each row 1 – Prob of Win – Prob of a tie, since the three events should (must) sum up to one. It is interesting to note that the probability of a tie is more or less near 12-14% and that, as we expected, it is equal to 0 when the hand is equal to 16 (since the dealer must draw on 16 and therefore a tie is not possible).


Here are the plots:


Rplot


Rplot01


Rplot02


This is the end of part 1 of the simulation.


 




And now the second part. Here we want to know what are the probabilities of not bust when we draw 1 or 2 cards given a certain hand.

# In this example we calculate the probabilities of bust
# and survival by asking for 1 or 2 cards given a hand
# Auxiliary functions are set first
# This function gives a player hand after n_draw given a hand,
# n decks, and a fixed number of draw (hits). Note: hand is a vector
hand_to_try <- function(n_decks, hand, n_draw)
{
deck_to_use <- deck_in_use(n_decks)
# Player's hand gets removed from the deck, since it is no longer available
for(i in hand)
{
# Remove aces
if((i == 11) | (i == 1))
{
deck_to_use <- deck_to_use[-match("A", deck_to_use)]
}else
{
deck_to_use <- deck_to_use[-match(i, deck_to_use)]
}
}
# Dealer draws
dealer_hand <- hand(deck_to_use, hand)[[2]]
# Dealer's hand gets removed from the deck, since it is no longer available
for(i in dealer_hand)
{
# remove aces
if((i == 11) | (i == 1))
{
deck_to_use <- deck_to_use[-match("A", deck_to_use)]
}else
{
deck_to_use <- deck_to_use[-match(i, deck_to_use)]
}
}
# Assuming hand is <= 11, n_draw cards are dealt and added to player's hand
for(i in 1:n_draw)
{
new_card <- sample(deck, 1, F)
if(new_card == "A")
{
new_card <- 1
new_card <- as.numeric(new_card)
}
hand <- append(hand, new_card)
}
# the function returns the final hand for the player
return(hand)
}
# Here we try the function
hand_to_try(1,c(10,"A"),2)
# This function simulates n_games trial at drawing n_draw
# cards given a player's hand with n_decks BJ game
# and returns the matrix with all the hands
m_pbusted_afterndraw <- function(n_games, n_decks, hand, n_draw)
{
hand_matrix <- matrix(ncol = length(hand) + n_draw)
for(i in 1:n_games)
{
p_hand <- hand_to_try(n_decks, hand, n_draw)
hand_matrix <- rbind(hand_matrix, p_hand)
}
#View(hand_matrix)
hand_matrix <- hand_matrix[-1,]
#View(hand_matrix)
return(hand_matrix)
}
# This code generates a matrix of 200 hands with two draws
# with a hand of 15 and 2 decks
matrix_hands <- m_pbusted_afterndraw(200, 2, c(10,5),2)
# This function, given a matrix of hands as matrice_hands,
# simply calculates the frequency (probability) of bust and not bust.
probab_bust <- function(matrix)
{
result_matrix <- matrix(ncol = 2)
colnames(result_matrix) <- c("Busted", "Not busted")
#In the matrix, 1 = True, 0 = False
for(i in 1:nrow(matrix))
{
if(sum(as.numeric(matrix[i,])) < 21)
{
result_matrix <- rbind(result_matrix, c(1,0))
}else
{
result_matrix <- rbind(result_matrix, c(0,1))
}
}
result_matrix <- result_matrix[-1,]
#View(result_matrix)
probab_vect <- c(sum(result_matrix[,1])/nrow(result_matrix), sum(result_matrix[,2])/nrow(result_matrix))
names(probab_vect) <- c("Probability of Bust", "Probability of survival")
return(probab_vect)
}
probab_bust(matrix_hands)
# The following code simulates 1000 hands each, with 3 decks,
# a hand (from 10 to 17) and one hit (draw), they return a matrix each
matrice_hands_hand_11_1draw <- m_pbusted_afterndraw(1000, 3, c(9,2),1)
matrice_hands_hand_12_1draw <- m_pbusted_afterndraw(1000, 3, c(9,3),1)
matrice_hands_hand_13_1draw <- m_pbusted_afterndraw(1000, 3, c(9,4),1)
matrice_hands_hand_14_1draw <- m_pbusted_afterndraw(1000, 3, c(9,5),1)
matrice_hands_hand_15_1draw <- m_pbusted_afterndraw(1000, 3, c(9,6),1)
matrice_hands_hand_16_1draw <- m_pbusted_afterndraw(1000, 3, c(9,7),1)
matrice_hands_hand_17_1draw <- m_pbusted_afterndraw(1000, 3, c(9,8),1)
# The following code simulates 1000 hands each, with 3 decks,
# a hand (from 10 to 17) and two hits (2 new cards), they return a
# matrix each
matrice_hands_hand_10_2draw <- m_pbusted_afterndraw(1000, 3, c(9,1),2)
matrice_hands_hand_11_2draw <- m_pbusted_afterndraw(1000, 3, c(9,2),2)
matrice_hands_hand_12_2draw <- m_pbusted_afterndraw(1000, 3, c(9,3),2)
matrice_hands_hand_13_2draw <- m_pbusted_afterndraw(1000, 3, c(9,4),2)
matrice_hands_hand_14_2draw <- m_pbusted_afterndraw(1000, 3, c(9,5),2)
matrice_hands_hand_15_2draw <- m_pbusted_afterndraw(1000, 3, c(9,6),2)
matrice_hands_hand_16_2draw <- m_pbusted_afterndraw(1000, 3, c(9,7),2)
matrice_hands_hand_17_2draw <- m_pbusted_afterndraw(1000, 3, c(9,8),2)
# The following code calculates the probabilities of bust for
# the two sets of arrays
prob_bust_hand_11_1draw <- probab_bust(matrice_hands_hand_11_1draw)
prob_bust_hand_12_1draw <- probab_bust(matrice_hands_hand_12_1draw)
prob_bust_hand_13_1draw <- probab_bust(matrice_hands_hand_13_1draw)
prob_bust_hand_14_1draw <- probab_bust(matrice_hands_hand_14_1draw)
prob_bust_hand_15_1draw <- probab_bust(matrice_hands_hand_15_1draw)
prob_bust_hand_16_1draw <- probab_bust(matrice_hands_hand_16_1draw)
prob_bust_hand_17_1draw <- probab_bust(matrice_hands_hand_17_1draw)
prob_bust_hand_10_2draw <- probab_bust(matrice_hands_hand_10_2draw)
prob_bust_hand_11_2draw <- probab_bust(matrice_hands_hand_11_2draw)
prob_bust_hand_12_2draw <- probab_bust(matrice_hands_hand_12_2draw)
prob_bust_hand_13_2draw <- probab_bust(matrice_hands_hand_13_2draw)
prob_bust_hand_14_2draw <- probab_bust(matrice_hands_hand_14_2draw)
prob_bust_hand_15_2draw <- probab_bust(matrice_hands_hand_15_2draw)
prob_bust_hand_16_2draw <- probab_bust(matrice_hands_hand_16_2draw)
prob_bust_hand_17_2draw <- probab_bust(matrice_hands_hand_17_2draw)
# The final array for 1 draw is built
bust_p_array_1draw <- matrix(ncol = 2)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_11_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_12_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_13_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_14_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_15_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_16_1draw)
bust_p_array_1draw <- rbind(bust_p_array_1draw, prob_bust_hand_17_1draw)
bust_p_array_1draw <- bust_p_array_1draw[-1,]
View(bust_p_array_1draw)
# The final array for 2 draw is built
bust_p_array_2draw <- matrix(ncol = 2)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_11_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_12_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_13_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_14_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_15_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_16_2draw)
bust_p_array_2draw <- rbind(bust_p_array_2draw, prob_bust_hand_17_2draw)
bust_p_array_2draw <;- bust_p_array_2draw[-1,]
View(bust_p_array_2draw)
hands_points <- c(11,12,13,14,15,16,17)
plot(hands_points, bust_p_array_1draw[,1], xlab = "Starting hand", main = "Probability of bust \ndrawing 1 card given \na BJ hand", type = "b", col="red")
plot(hands_points, bust_p_array_1draw[,2], xlab = "Starting hand", main = "Probability of survival \ndrawing 1 card given \na BJ hand", type = "b", col="blue")
plot(hands_points, bust_p_array_2draw[,1], xlab = "Starting hand", main = "Probability of bust \ndrawing 2 cards given \na BJ hand", type = "b", col="red")
plot(hands_points, bust_p_array_2draw[,2], xlab = "Starting hand", main = "Probability of survival \ndrawing 2 cards given \na BJ hand", type = "b", col="blue")

Let’s have a look at the new results we have got:
This is a hand we have got by asking 2 times card from a single deck.


im7


with a starting hand equal to (10,”A”) if we were to ask 2 cards, here are the probabilities that we should face:


im8


Eventually, we can compute the probabilities of “survival” and bust for each hand from 16 to 17 and plot the results.


im9


im10


Rplot03


Rplot04


Rplot05


Rplot06


The results looks pretty interesting and it is nice to see that the probability of bust if we ask for 1 card and our starting hand is 11 is 0 as we could easily predict, however, if we were to ask 2 cards, our chances of survival would suddenly drop to 25%. A big leap!


Hope this was interesting.





Disclaimer: This article is for educational purpose ONLY. Odds generated by this code are calculated by a random simulation. As such the odds will represent an approximation of the true odds. They might even be completely wrong or misleading. This code must NOT be used for anything other than educational purpose. The provider of this code does not guarantee the accuracy of the results and accepts no liability for any loss or damage that may occur as a result of the use of this code. Understanding and agreeing to the terms of this disclaimer is a condition of use of this code. By reading the article you confirm you have understood and will comply with this disclaimer.

To leave a comment for the author, please follow the link and comment on their blog: The Beginner Programmer.

R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

Never miss an update!
Subscribe to R-bloggers to receive
e-mails with the latest R posts.
(You will not see this message again.)

Click here to close (This popup will not appear again)