I often want to simulate outcomes varying across a set of

`parameters.`* *In order to accomplish this in an efficient

manner I have coded up a little function that takes parameter

vectors and produces results. First I will show how to set it

up with some dummy examples and next I will show how it can

be used to select the optimal blackjack strategy.

SimpleSim <- function(..., fun, pairwise=F) { # SimpleSim allows for the calling of a function varying

# multiple parameters entered as vectors. In pairwise form

# it acts much like apply. In non-paiwise form it makes a

` # combination of each possible parameter mix`

# in a manner identical to block of nested loops. returner <- NULL L <- list(...) # Construct a vector that holds the lengths of each object vlength <- unlist(lapply(L, length)) npar <- length(vlength) CL <- lapply(L, "[", 1) # Current list is equal to the first element # Pairwise looping if (pairwise) { # If pairwise is selected than all elements greater than 1 must be equal. # Checks if all of the elements of a vector are equal if (!(function(x) all(x[1]==x))(vlength[vlength>1])) { print(unlist(lapply(L, length))) stop("Pairwise: all input vectors must be of equal length", call. =F) } for (i in 1:max(vlength)) { # Loop through calling the function CL[vlength>1] <- lapply(L, "[", i)[vlength>1] # Current list returner <- rbind(returner,c(do.call(fun, CL),pars="", CL)) } } # End Pairwise # Non-pairwise looping if (!pairwise) { ncomb <- prod(vlength) # Calculate the number of combinations print(paste(ncomb, "combinations to loop through")) comb <- matrix(NA, nrow=prod(vlength), ncol=npar+1) comb[,1] <- 1:prod(vlength) # Create an index value comb <- as.data.frame(comb) # Converto to data.frame colnames(comb) <- c("ID", names(CL)) for (i in (npar:1)) { # Construct a matrix of parameter combinations comb[,i+1] <- L[[i]] # Replace one column with values comb<-comb[order(comb[,(i+1)]),] # Reorder rows } comb<-comb[order(comb[,1]),] for (i in 1:ncomb) { for (ii in 1:npar) CL[ii] <- comb[i,ii+1] returner <- rbind(returner,c(do.call(fun, CL),pars="", CL)) } } # End Non-Pairwise return(returner) } # END FUNCTION DEFINITION # Let's first define a simple function for demonstration minmax <- function(...) c(min=min(...),max=max(...)) # Pairwise acts similar to that of a multidimensional apply across columns SimpleSim(a=1:20,b=-1:-20,c=21:40, pairwise=T, fun="minmax") # The first set of columns are those of returns from the function "fun" called. # The second set divided by "par" are the parameters fed into the function. SimpleSim(a=1:20,b=-1:-20,c=10, pairwise=T, fun="minmax") # Non-pairwise creates combinations of parameter sets. # This form is much more resource demanding. SimpleSim(a=1:5,b=-1:-5,c=1:2, pairwise=F, fun="minmax") # Let's try something a little more interesting. # Let's simulate a game of black jack strategies assuming no card counting is possible. blackjack <- function(points=18, points.h=NULL, points.ace=NULL, cards=10, cards.h=NULL, cards.ace=NULL, sims=100, cutoff=10) { # This function simulates a blackjack table in which the player # has a strategy of standing (not asking for any more cards) # once he has either recieved a specific number of points or # a specific number of cards. This function repeates itself sims # of times. # This function allows for up to three different strategies to be played. # 1. If the dealer's hole card is less than the cuttoff # 2. If the dealer's hole card is greater than or equal to the cuttoff # 3. If the dealer's hole card is an ace # In order to use 3 level strategies input parameters as .h and .ace # It returns # of wins, # of losses, # of pushes (both player and dealer gets 21) # and the number of blackjacks. # This simulation assumes the number of decks used is large thus # the game is like drawing with replacement. if (is.null(points.h)) points.h <- points if (is.null(points.ace)) points.ace <- points.h if (is.null(cards.h)) cards.h <- cards if (is.null(cards.ace)) cards.ace <- cards.h bdeck <- c(11,2:9,10,10,10,10) # 11 is the ace bdresult <- c(ppoints=NULL, pcards=NULL, dpoints=NULL, dcards=NULL) for (s in 1:sims) { dhand <- sample(bdeck,1) # First draw the deal's revealed card phand <- sample(bdeck,2, replace=T) # Specify target's based on dealer's card if (dhand<cutoff) { pcuttoff <- points ccuttoff <- cards } if (dhand>=cutoff) { pcuttoff <- points.h ccuttoff <- cards.h } if (dhand==11) { pcuttoff <- points.ace ccuttoff <- cards.ace } # player draws until getting above points or card count while ((sum(phand)<pcuttoff)&(length(phand)<ccuttoff)){ phand <- c(phand, sample(bdeck,1)) # If player goes over then player may change aces to 1s if (sum(phand)>21) phand[phand==11] <- 1 } # Dealer must always hit 17 so hand is predetermined while (sum(dhand)<17) { dhand <- c(dhand, sample(bdeck,1)) # If dealer goes over then dearler may change aces to 1s if (sum(dhand)>21) dhand[dhand==11] <- 1 } bdresult <- rbind(bdresult, c(ppoints=sum(phand), pcards=length(phand), dpoints=sum(dhand), dcards=length(dhand))) } # Calculate the times that the player wins, pushes (ties), and loses pbj <- (bdresult[,1]==21) & (bdresult[,2]==2) dbj <- (bdresult[,3]==21) & (bdresult[,4]==2) pwins <- ((bdresult[,1] > bdresult[,3]) & (bdresult[,1] < 22)) | (pbj & !dbj) push <- (bdresult[,1] == bdresult[,3]) | (pbj & dbj) dwins <- !(pwins | push) # Specify the return. c(odds=sum(pwins)/sum(dwins), pwins=sum(pwins), dwins=sum(dwins), push=sum(push), pcards=mean(bdresult[,2]), dcards=mean(bdresult[,4]), pblackjack=sum(pbj), dblackjack=sum(dbj)) } blackjack(points=18, sims=4000) # We can see unsurprisingly, that the player is not doing well. blackjack(points=18, points.h=19, sims=4000) # We can see that by adopting a more aggressive strategy for when # the dealer has a 10 point card or higher, we can do slightly better. # But overall, the dealer is still winning about 3x more than us. # We could search through different parameter combinations manually to # find the best option. Or we could use our new command SimpleSim! MCresults <- SimpleSim(fun=blackjack, points=15:21, points.h=18:21, points.ace=18:21, cutoff=9:10, cards=10, sims=100) # Let's now order our results from the most promising. MCresults[order(-unlist(MCresults[,1])),] # By the simulation it looks like we have as high as a 50% ratio of loses to wins. # Which means for every win there are 2 loses. # However, I don't trust it since we only drew 100 simulations. # In addition, this is the best random draw from all 224 combinations which each # have different probabilities. # Let's do the same simulation but with 2000 draws per. # This might take a little while. MCresults <- SimpleSim(fun=blackjack, points=15:21, points.h=18:21, points.ace=18:21, cutoff=9:10, cards=10, sims=5000) # Let's now order our results from the most promising. MCresults[order(-unlist(MCresults[,1])),]

hist(unlist(MCresults[,1]), main="Across all combinations\nN(Win)/N(Loss)",

xlab = "Ratio", ylab = "Frequency")

# The best case scenario 38% win to loss ratio appears around were we started, # playing to hit 18 always and doing almost as well when the dealer is high # (having a 10 or ace) then playing for 19. # Overall, the odds are not in our favor. For every win we expect 1/.38 (2.63) loses.

Along similar lines....

ReplyDeletehttp://www.unt.edu/rss/class/Jon/R_SC/Module10/CardsBeta.R

~Jon

woooow! that script is great

ReplyDelete