# The Three Little Pigs

The game of pig has simple rules but complex strategies. It was described for the first time in 1945  by a magician called John Scarne. Playing the pig game is easy: each turn, a player repeatedly rolls a die until either a 1 is rolled or the player decides to hold:

• If the player rolls a 1, they score nothing and it becomes the next player’s turn
• If the player rolls any other number, it is added to their turn total and the player’s turn continues
• If a player chooses to hold, their turn total is added to their score, and it becomes the next player’s turn

The first player who reach at least 100 points is the winner. For example: you obtain a 3 and then decide to roll again, obtaining a 1. Your score is zero in this turn. Next player gets the sequence 3-4-6 and decides to hold, obtaining a score of 13 points in this turn.

Despite its simplicity, the pig game has a very complex and against-intuition optimal strategy. It was calculated in 2004 by Todd W. Neller and Clifton Presser from Gettysburg College of Pennsilvania with the help of computers.

To illustrate the game, I simulated three players (pigs) playing the pig game with three different strategies:

• The Coward pig, who only rolls the die a small number of times in every turn
• The Risky pig, who rolls the die a more times than the coward one
• The Ambitious pig, who tries to obtain in every turn more points than two others

I simulated several scenarios.

• Some favorable scenarios for Coward pig:

In first scenario, the Coward pig rolls the die between 1 and 5 times each round and wins if the Risky pig asumes an excessive level of risk (rolling each time between 10 and 15 times). Trying to obtain more than the Coward is a bad option for the Ambitious pig. Simulating this scenario 100 times gives victory to Coward a 51% of times (25% to Risky and 24% to Ambitious).

Second scenario puts closer Coward and Risky pigs (first one rolls the die between 4 and 7 times  each round and second one between 6 and 9 times). Coward wins 54% of times (34% Risky and only 12% Ambitious).

Being coward seems to be a good strategy when you play against a reckless or when you are just a bit more conservative than a Risky one.

• Some favorable scenarios for Risky pig:

Rolling the die between 4 and 6 times each round seems to be a good option, even more when you are playing against a extremely conservative player who rolls no more than 3 times each time. Simulating 100 times these previous scenarios gives victory to Risky pig a 58% of times in first the case in which Coward rolls allways 1 and Risky 6 times each round (0% for Coward and only 42% form Ambitious) and 66% of times in the second one (only 5% to Coward and 29% to Ambitious).

Being Risky is a good strategy when you play against a chicken.

• Some favorable scenarios for Ambitious pig:

The Ambitious pig wins when two others turn into extremely coward and risky pigs as can be seen in the first scenario in which Ambitious wins 65% of the times (31% for Coward and 4% for Risky). Ambitious pig also wins when two others get closer and hit the die a small number of times (2 rolls the Coward and 4 rolls the Risky). In this scenario the Ambitious wins 58% of times (5% for Coward and 37% for Risky). By the way, these two scenarios sound very unreal.

Being ambitious seems to be dangerous but works well when you play against a crazy and a chicken or against very conservative players.

From my point of view, this is a good example to experiment with simulations, game strategies and xkcd style graphics.

The code:

```require(ggplot2)
require(extrafont)
#Number of hits for Coward
CowardLower=2
CowardUpper=2
#Number of hits for Risky
RiskyLower=4
RiskyUpper=4
game=data.frame(ROUND=0, part.p1=0, part.p2=0, part.p3=0, Coward=0, Risky=0, Ambitious=0)
while(max(game\$Coward)<100 & max(game\$Risky)<100 & max(game\$Ambitious)<100)
{
#Coward Little Pig
p1=sample(1:6,sample(CowardLower:CowardUpper,1), replace=TRUE)
s1=min(min(p1-1),1)*sum(p1)
#Risky Little Pig
p2=sample(1:6,sample(RiskyLower:RiskyUpper,1), replace=TRUE)
s2=min(min(p2-1),1)*sum(p2)
#Ambitious Little Pig
s3=0
repeat {
p3=sample(1:6,1)
s3=(p3+s3)*min(min(p3-1),1)
if (p3==1|s3>max(s1,s2)) break
}
game[nrow(game)+1,]=c(max(game\$ROUND)+1,s1,s2,s3,max(game\$Coward)+s1,max(game\$Risky)+s2,max(game\$Ambitious)+s3)
}
opts=theme(
panel.background = element_rect(fill="darkolivegreen1"),
panel.border = element_rect(colour="black", fill=NA),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour="black"),
panel.grid.major = element_line(colour="white", linetype = 1),
panel.grid.minor = element_blank(),
axis.text.y = element_text(colour="black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=25, family="xkcd"),
legend.key = element_blank(),
legend.position = c(.2,.75),
legend.background = element_blank(),
plot.title = element_text(size = 50)
)
ggplot(game, mapping=aes(x=game\$ROUND, y=game\$Coward)) +
geom_line(color="red", size=1.5) +
geom_line(aes(x=game\$ROUND, y=game\$Risky), color="blue", size=1.5) +
geom_line(aes(x=game\$ROUND, y=game\$Ambitious), color="green4", size=1.5) +
geom_point(aes(x=game\$ROUND, y=game\$Coward, colour="c1"), size=5.5) +
geom_point(aes(x=game\$ROUND, y=game\$Risky, colour="c2"), size=5.5) +
geom_point(aes(x=game\$ROUND, y=game\$Ambitious, colour="c3"), size=5.5) +
ggtitle("THE THREE LITTLE PIGS") +
xlab("ROUND") + ylab("SCORING") +
geom_text(aes(max(game\$ROUND), max(max(game\$Coward, game\$Risky, game\$Ambitious)), hjust=1.2, family="xkcd", label="WINNER!"), size=10)+
geom_hline(yintercept=100, linetype=2, size=1)+
scale_y_continuous(breaks=seq(0, max(max(game\$Coward, game\$Risky, game\$Ambitious))+10, 10))+
scale_x_continuous(breaks=seq(0, max(game\$ROUND), 1))+
scale_colour_manual("",
labels = c(paste("Coward: ", CowardLower, "-", CowardUpper, " hits", sep = ""), paste("Risky: ", RiskyLower, "-", RiskyUpper, " hits", sep = ""), "Ambitious"),
breaks = c("c1", "c2", "c3"),
values = c("red", "blue", "green4"))+ opts
```

# Do Not Play With Mr. Penney

Facts do not speak (Henry Poincare)

Mr. Penney is my best friend. He is maths teacher and loves playing. Yesterday we were in his office at the university when he suggested me a game:

When you toss a coin three times, you can obtain eight different sequences of tails and heads: TTT, TTH, THT, HTT, THH, HTH, HHT and HHH. Using a fair coin, all sequences have the same chances to appear. Choose one sequence and I will then choose another one. I will toss a coin until either your or my sequence appears as a consecutive subsequence of the coin toss outcomes. The player whose sequence appears first wins. I will repeat this procedure 100 times. The one with more games won is the winner of the game.  Don’t worry: I will not toss the coin manually. I will simulate using my computer. What’s your bet?

Ok, my bet is THT, I said. After some seconds, Mr. Penney said: My bet is TTH.

This was the result of the first round:

Another chance? told me Mr. Penney. Of course! Now my bet is TTH! I said. In fact, I was thinking Take that! Now I chose your previous bet. Do you think I am foolish?. After some seconds, Mr. Penney said: My bet now is HTT.

This was the result of the second round:

Another chance? told me Mr. Penney. At this point, I was very suspicious but I wanted the last chance so I told him Of course! Now my bet is HTT! I wanted to try my strategy one more time. After some seconds, Mr. Penney said: My bet now is HHT.

This was the result of the third round:

Ok, I give it up! What’s the trick? I said. And Mr. Penney explained it to me. You can find the explanation here. This is the last time I play with you! I told him once he finished the explanation.

Here you have the code. Feel free to play:

```library(gridExtra)
library(gridExtra)
Me     <- "TTH"
Penney <- "HTT"
results <- data.frame(play= numeric(0), Penney = integer(0), Me = character(0))
for (i in 1:100) {
play <- c()
repeat {play <- do.call(paste, c(play, as.list(sample(c("H","T"), 1)), sep=""))
if (grepl(Penney, play)|grepl(Me, play)) {
results <- rbind(results, data.frame(play= i, Penney = as.numeric(grepl(Penney, play)), Me = as.numeric(grepl(Me, play))))
break}}}
grid.newpage()
table <- rbind(
c("Me", Me, sum(results\$Me), if(sum(results\$Penney) > sum(results\$Me)) "Loser" else "Winner"),
c("Penney", Penney, sum(results\$Penney), if(sum(results\$Penney) > sum(results\$Me)) "Winner" else "Loser"))
grid.table(table,
cols = c("Player", "Bet", "Games Won", "Result"),
gpar.colfill = gpar(fill="palegreen3",col="White"),
gpar.corefill =  gpar(fill="palegreen",col="White"),
gpar.rowfill = gpar(fill=NA, col=NA))
```

# Random Love

Anyone who considers arithmetical methods of producing random digits is, of course, in a state of sin (John von Newman)

Ms. Positive and Mr. Negative live in a one-dimensional world and are falling in love. But beginnings are not always easy. They have a big problem: none of them like the other’s neighborhood. Ms. Positive only wants to walk around Positive Integer Numbers Neighborhood and Mr. Negative around Negative Integers Numbers one. This is a prickly problem they need to deal with as soon as possible. But they have a good idea. They will start their walks from Zero, an impartial place between both neighborhoods and will let fate to guide their feet. They will toss a coin to decide every step: if result is head, they will advance 1 step toward positive numbers neighborhood; if tail, they will advance 1 step toward negative numbers one. For example, if the first 5 tosses are face, face, tail, tail and tail, the their first 5 steps will be +1, +2, +1, 0 and -1. It seems to be a fair agreement for both. Maybe is not the most pleasant way to take a walk but It is well known that lovers use to do silly things constantly, especially at the beginnings. They always walk for two hours, so they toss the coin 7.200 times every walk (these lovers are absolutely crazy as you can see). This was their first walk:

After this first walk, Mr Negative was really upset. Ms. Positive, watching his face fell, ask him: What’s the matter, honey? and Mr. Negative replied: What’s the matter? What’s the matter? The matter is that we spent almost all the time walking around your horrible neighborhood! What comes next is too hard to be reproduced here. Anyway, they agreed to give a chance to the method they designed. How can one imagine that a coin can produce such a strange walk! There must be an error! After 90 walks, the situation of our lovers was extremely delicate. A 57% of the walks were absolutely awful for one of them since more than 80% of the steps were around the same neighborhood. Another 32% were a bit uncomfortable for one of them since between 60% and 80% of the steps were around the same neighborhood. Only 11% of the walks were gratifying. How is it possible?, said Mr. Negative. How is it possible?, said Ms. Positive.

But here comes Ms. Positive, who always looks on the brigth side of life: Don’t worry, darling. In fact, we don’t have to be sad. We get angry the same amount of times! For me is enough. What about you?, said her. For me is perfect as well!, said Mr. Negative. In that moment, they realise they were made for each other and started another random walk with a big smile on their faces.

This is the code:

```library(ggplot2)
steps   <- 2*60*60 #Number of steps
results <- data.frame()
walks<-90 #Number of walks
for (i in 1:walks)
{
state <- cumsum(sample(c(-1,1), steps, replace = TRUE))
results <- rbind(results, c(sum(state<0), sum(state>0), sum(state==0),
if (sum(state<0) >= sum(state>0)) 1 else 0))
}
colnames(results) <- c("neg.steps", "pos.steps", "zero.steps", "ind.neg")
results\$max.steps <- apply(results, 1, max)/apply(results, 1, sum)
#Plot of one of these walks
mfar=max(abs(max(state)),abs(min(state)))
plot1 <- qplot(seq_along(state),
state,
geom="path")+
xlab("Step") +
ylab("Location") +
labs(title = "The First Walk Of Ms. Positive And Mr. Negative")+
theme(plot.title = element_text(size = 35))+
theme(axis.title.y = element_text(size = 20))+
theme(axis.title.x = element_text(size = 20))+
scale_x_continuous(limits=c(0, length(state)),breaks=c(1,steps/4,steps/2,3*steps/4,steps))+
scale_y_continuous(limits=c(-mfar, mfar), breaks=c(-mfar,-mfar/2, 0, mfar/2,mfar))+
geom_hline(yintercept=0)
ggsave(plot1, file="plot1.png", width = 12, height = 10)
#Summary of all walks
hist1 <- ggplot(results, aes(x = max.steps))+
geom_histogram(colour = "white",breaks=seq(.4,1,by=.2),fill=c("blue", "orange", "red"))+
theme_bw()+
labs(title = paste("What Happened After ", toString(walks), " Walks?",sep = ""))+
scale_y_continuous(breaks=seq(0,(nrow(results[results\$max.steps>.8,])+10),by=10))+
theme(plot.title = element_text(size = 40))+
xlab("Maximum Steps In The Same Location (%)") +
ylab("Number of Walks")
ggsave(hist1, file="hist1.png", width = 10, height = 8)
#Data for waterfall chart
waterfall <- as.data.frame(cbind(
c("Total Walks", "Satisfactory Walks", "Uncomfortable Walks", "Awful Walks for Mr. +", "Awful Walks for Ms. -"),
c("a", "b", "c", "d", "d"),
c(0,
nrow(results),
nrow(results)-nrow(results[results\$max.steps<.6,]),
nrow(results)-nrow(results[results\$max.steps<.6,])-nrow(results[results\$max.steps>=.6 & results\$max.steps<.8,]),
nrow(results)-nrow(results[results\$max.steps<.6,])-nrow(results[results\$max.steps>=.6 & results\$max.steps<.8,])-nrow(results[results\$max.steps>=.8 & results\$ind.neg==1,])
),
c(nrow(results),
nrow(results)-nrow(results[results\$max.steps<.6,]),
nrow(results)-nrow(results[results\$max.steps<.6,])-nrow(results[results\$max.steps>=.6 & results\$max.steps<.8,]),
nrow(results)-nrow(results[results\$max.steps<.6,])-nrow(results[results\$max.steps>=.6 & results\$max.steps<.8,])-nrow(results[results\$max.steps>=.8 & results\$ind.neg==1,]),
0
),
c(nrow(results),
nrow(results[results\$max.steps<.6,]),
nrow(results[results\$max.steps>=.6 & results\$max.steps<.8,]),
nrow(results[results\$max.steps>=.8 & results\$ind.neg==1,]),
nrow(results[results\$max.steps>=.8 & results\$ind.neg==0,]))
))
colnames(waterfall) <-c("desc", "type", "start", "end", "amount")
waterfall\$id <- seq_along(waterfall\$amount)
waterfall\$desc <- factor(waterfall\$desc, levels = waterfall\$desc)
#Waterfall chart
water1 <- ggplot(waterfall, aes(desc, fill = type)) +
geom_rect(aes(x = desc, xmin = id-0.45, xmax = id+0.45, ymin = end, ymax = start))+
xlab("Kind of Walk") +
ylab("Number of Walks") +
labs(title = "The Ultimate Proof (After 90 Walks)")+
theme(plot.title = element_text(size = 35))+
theme(axis.title.y = element_text(size = 20))+
theme(axis.title.x = element_text(size = 20))+
theme(legend.position = "none")
```

# Dora’s Choice

Arithmetic is being able to count up to twenty without taking off your shoes (Mickey Mouse)

On her last mission, Dora The Explorer sails down the Amazon river to save her friend Isa The Iguana from Swiper The Fox claws. After some hours of navigation, Dora sees how the river divides into 3 branches and has to choose which one to follow. Before leaving, her friend Map told her that just one of these branches is safe. Two others end in terrible waterfalls, both impossible to escape alive. Although Dora does not know which one is the good one, she decides to take the branch number 1. Suddenly, her friend Boots The Monkey yells from the top of a palm tree:

Dora, do not take branch number 3! I can see from here that it ends in a horrible waterfall!

After listening to Boots, Dora changes her mind and decides to take branch number 2. Why Dora switches? Because she knows that this change has significantly increased her probability of ending the mission alive.

There are several ways to convince yourself of this. One is to simulate the situation that has faced Dora and compare results of switching and not switching . Switching, Dora saves her life 2 of each 3 simulations while if she does not, Dora only saves 1 of each 3 simulations. Changing her mind, Dora doubles her chances of survival!

Carefully considering what happens, you can see that switching Dora saves herself when her first choice is erroneus, which occurs with probability 2/3. On the other hand, if Dora remains faithful to her first choice, obviously only saves herself with probability 1/3.

This is an example on my own of the famous Monty Hall Problem. You can see a nice explanation of it in a chapter of Numb3rs or in the film 21 Black Jack. Not long ago I exposed the problem in a family meeting. Only my mum said she would switch (we were 6 people in the meeting). It is fun to share this experiment and ask what people would do. Do it with your friends and family. First time I knew the problem I thought there were no difference between switching and not since I gave both possibilities 1/2 of probability. If I had been Dora, pretty sure I would tumbled over a terrible waterfall. What about yo?

Note: this is an update of the post, which was not a correct formulation of Monty Hall Problem. Thanks to David Robinson and Scott Kostyshak for showing me my error. A correct formulation of the problem may be this:

On her last mission, Dora The Explorer sails down the Amazon river to meet her cousin Diego. After some hours of navigation, Dora sees how the river divides into 3 branches and has to choose which one to follow. Before leaving, her friend Map told her that just one of these branches is safe. Two others end in terrible waterfalls, both impossible to escape alive. Although Dora does not know which one is the good one, she decides to take the branch number 1. After putting the bow towards branch number one, Dora sees Swiper The Fox smiling from the shore, in a high place where obviously can see the end of all three branches. Dora yells him:

– Help me Swiper! Which one should I take?

Swiper replies:

– I am the villain of this story so I will give you only an advice: do not take branch number 3. It ends into a terrible waterfall.

Dora, who has a sixth sense to notice when Swiper is lying, knows he is telling the truth and immediately changes her mind and decides to take branch number 2. Why Dora switches? Because she knows that this change has significantly increased her probability of ending the mission alive.

Here you have the code:

```library(ggplot2)
library(extrafont)
nchoices <- 3
nsims <- 500
choices <- seq(from=1, to=nchoices, by=1)
good.choice <- sample(choices, nsims, replace=TRUE)
choice1 <- sample(choices, nsims, replace=TRUE)
dfsims <- as.data.frame(cbind(good.choice, choice1))
dfsims\$advice <- apply(dfsims, 1, function(x) choices[!choices %in% as.vector(x)][sample(1:length(choices[!choices %in% as.vector(x)]), 1)])
dfsims\$choice2 <- apply(dfsims, 1, function(x) choices[!choices %in% as.vector(c(x[2], x[3]))][sample(1:length(choices[!choices %in% as.vector(c(x[2], x[3]))]), 1)])
dfsims\$win1 <- apply(dfsims, 1, function(x) (x[1]==x[2])*1)
dfsims\$win2 <- apply(dfsims, 1, function(x) (x[1]==x[4])*1)
dfsims\$csumwin1 <- cumsum(dfsims\$win1)/as.numeric(rownames(dfsims))
dfsims\$csumwin2 <- cumsum(dfsims\$win2)/as.numeric(rownames(dfsims))
dfsims\$nsims <- as.numeric(rownames(dfsims))
dfsims\$xaxis <- 0
### XKCD theme
theme_xkcd <- theme(
panel.background = element_rect(fill="darkolivegreen1"),
panel.border = element_rect(colour="black", fill=NA),
axis.line = element_line(size = 0.5, colour = "black"),
axis.ticks = element_line(colour="black"),
panel.grid = element_line(colour="white", linetype = 2),
axis.text.y = element_text(colour="black"),
axis.text.x = element_text(colour="black"),
text = element_text(size=18, family="Humor Sans"),
plot.title = element_text(size = 50)
)
### Plot the chart
p <- ggplot(data=dfsims, aes(x=nsims, y=csumwin1))+
geom_line(aes(y=csumwin2), colour="green4", size=1.5, fill=NA)+
geom_line(colour="green4", size=1.5, fill=NA)+
geom_text(data=dfsims[400, ], family="Humor Sans", aes(x=nsims), colour="green4", y=0.7, label="if Dora switches ...", size=5.5, adjust=1)+
geom_text(data=dfsims[400, ], family="Humor Sans", aes(x=nsims), colour="green4", y=0.3, label="if Dora does not switch ...", size=5.5, adjust=1)+
coord_cartesian(ylim=c(0, 1), xlim=c(1, nsims))+
scale_y_continuous(breaks = c(0,round(1/3, digits = 2),round(2/3, digits = 2),1), minor_breaks = c(round(1/3, digits = 2),round(2/3, digits = 2)))+
scale_x_continuous(minor_breaks = seq(100, 400, 100))+
labs(x="Number Of Simulations", y="Rate Of Survival", title="Dora's Choice")+
theme_xkcd
ggsave("doras_choice.jpg", plot=p, width=8, height=5)
```

# Shoot The Heart With Monte Carlo

The heart has its reasons which reason knows not (Blaise Pascal)

You only need two functions to draw a heart mathematically. The upper part is generated by `(1-(|x|-1)2)1/2` and the lower one by `acos(1-|x|)-PI`. Here is how this heart is:

Whats the area of this heart? It’s easy: integrating `heart.up(x)-heart.dw(x)` between -2 and 2 and you will obtain that heart measures 9.424778, but there is a simple and nice way to approximate to this value: shoot the heart.

The idea is very simple. Heart is delimited by a square with vertex in `(-2, heart.dw(0))`, `(-2, 1)`, `(2, heart.dw(0))` and `(2, 1)`. Generating a set of points uniformly distributed inside the square and counting how many of them fall into the heart in relation to the area of the square gives a very good approximation of the exact area of the heart. This is a plot representing a simulation of 2.000 shots (hits in red, fails in blue):

Given a simulation of n points, the estimated area of the heart is the area of the square by percentage of points that falls inside the heart. And of course, precision increases with the number of shots you make, as you can see in the following plot, where exact area is represented by the red horizontal line:

Here you have the code:

```library("ggplot2")
heart.up <- function(x) {sqrt(1-(abs(x)-1)^2)} #Upper part of the heart
heart.dw <- function(x) {acos(1-abs(x))-pi}    #Lower part of the heart
#Plot of the heart
ggplot(data.frame(x=c(-2,2)), aes(x)) +
stat_function(fun=heart.up, geom="line", aes(colour="heart.up")) +
stat_function(fun=heart.dw, geom="line", aes(colour="heart.dw")) +
scale_colour_manual("Function", values=c("blue","red"), breaks=c("heart.up","heart.dw"))+
labs(x = "", y = "")+
theme(legend.position = c(.85, .15))
sims <- 2000 #Number of simulations
rlts <- data.frame()
for (i in 1:sims)
{
msm <- cbind(as.matrix(runif(i, min=-2, max=2)), as.matrix(runif(i, min=heart.dw(0), max=1)))
nin <- 0
for (j in 1:nrow(msm)) {if (msm[j,2]<=heart.up(msm[j,1]) & msm[j,2]>=heart.dw(msm[j,1])) nin=nin+1}
rlts <- rbind(c(i, 4*(1-heart.dw(0))*nin/i), rlts)
}
colnames(rlts) <- c("no.simulations","heart.area")
exact.area <- integrate(function(x) {heart.up(x)-heart.dw(x)},-2,2)\$value
mean.area <- mean(rlts\$heart.area) #Mean of All Estimated Areas
ggplot(data = rlts, aes(x = no.simulations, y = heart.area))+
geom_point(size = 0.5, colour = "black", alpha=0.4)+
geom_abline(intercept = exact.area, slope = 0, size = 1, linetype=1, colour = "red", aes(color="My Line"), alpha=0.8, show_guide = TRUE)+
labs(list(x = "Number of Shots", y = "Estimated Area"))+
ggtitle("Shot The Heart With Monte Carlo") +
theme(plot.title = element_text(size=20, face="bold"))+
scale_x_continuous(limits = c(0, sims), expand = c(0, 0))+
expand_limits(x = 0, y = 0)+
scale_y_continuous(limits = c(0, 2*exact.area), expand = c(0, 0), breaks=c(0, exact.area/4, exact.area/2, 3*exact.area/4, exact.area, 5*exact.area/4, 3*exact.area/2, 7*exact.area/4, 2*exact.area))+
geom_text(x = 1000, y = exact.area/2, label=paste("Exact Area =", sprintf("%7.6f", exact.area)), vjust=-1, colour="red", size=5)+
geom_text(x = 1000, y = exact.area/2, label=paste("Mean of All Estimated Areas=", sprintf("%7.6f", mean.area)), vjust=+1, colour="red", size=5)
```