Exam 1
Nick Strayer
October 20, 2015
Question 1:
What is the span or support of S?
\[[3, 42]\]
Question 2:
If rolling of all three dice, what is the probability of getting a sum of 3 on that roll, i.e. P(S=3)?
\[P(3) = \frac{1}{10} \cdot \frac{1}{12} \cdot \frac{1}{20} = 0.0004166667\]
Question 3:
If rolling of all three dice, what is the probability of getting a sum of 21 on that roll given the twenty-sided landed on 18, i.e. P(S=21 | D20=18)?
\[P(A|B) = \frac{P(A,B)}{P(B)}\]
To get a 21 with the D20 as an 18 you can either have (D10 = 1, D12 = 2), or (D10 = 2 , D12 = 1).
Probability of this can be calculated as such:
\[P(A|B) = \frac{P(A,B)}{P(B)}\] \[P(A,B) = \frac{1}{10} \cdot \frac{1}{12}\cdot \frac{1}{20} + \frac{1}{10} \cdot \frac{1}{12}\cdot \frac{1}{20}\] \[P(B) = \frac{1}{20}\] \[\frac{P(A,B)}{P(B)} = 0.01666667\]
Question 4:
If holding an 18 on the twenty-sided die and rolling the other two dice, what is the probability of getting a sum of 21 on that roll, i.e. P(S=21 | D20 held at 18)?
Due to the rolls all being independent this probability is simply the same as the previous one.
\[\frac{1}{10} \cdot \frac{1}{12} + \frac{1}{10} \cdot \frac{1}{12} = 0.01666667\]
Question 5:
If rolling all three dice, what is the expectation of the dice sum, i.e. E[S]? Hint: If you think of S as D10+D12+D20, you don’t need the probability mass function of S.
First we write a quick function to calculate the expected value of a dice with d sides.
diceE = function(d) sum((1/d)* (1:d))We can then use the fact that expectations are aditive (e.g. \(E[X_1 + X_2] = E[X_1] + E[X_2]\)) to calculate the expected value of all the dice combined.
(sumOfRoll <- diceE(10) + diceE(12) + diceE(20) ) ## [1] 22.5So \(E[S] =\) 22.5.
Question 6:
If rolling of all three dice, what is the standard deviation of the dice sum, i.e. SD[S]?
Another convenient fact about expectations is that the variance (aka \(E[X^2] - E[X]^2\)) is also additive. As such we can write another function for calculating \(E[X^2]\) and then the variance and be on our way.
#Function for E[X^2]
diceE2 = function(d) sum((1/d)* (1:d)^2)
#Function for Variance
varD = function(d) diceE2(d) - (diceE(d))^2
#Calculate the variance.
varOfRoll <- varD(10) + varD(12) + varD(20)
#Take the square root to get the Standard Deviation.
(SdOfRoll <- sqrt(varOfRoll) )## [1] 7.308671From this we can see that \(SD[S] =\) 7.3086706
Question 7:
Call strategy A simply rolling of all three dice every roll until getting a Dragon Jack (S=21). I wrote a script in R to play strategy A 17,000 times and saved the number of rolls until getting a Dragon Jack. The simulation’s sample mean and sample standard deviation were 20.31224 and 19.74028, respectively. A histogram of the simulation is below. (I have ommited the histogram.)
Calculate a 95% confidence interval expressed as (LB, UB) for the true mean number of rolls to get a Dragon Jack using strategy A. In addition, justify the confidence interval method you used in a couple short sentences.
sampleMean <- 20.31224
sampleSD <- 19.74028
n <- 17000I will choose to use a T interval for this sample as we don’t know the true standard deviation. Just plugging in the sample standard deviation into a normal interval will give almost the same results though because as the T’s degrees of freedom get large (and 17,000 is definitely large in T terms) it approximates the normal very closely.
alpha = 0.05 #95% interval
tVal <- qt(1 - (alpha/2), n - 1)
width <- tVal * (sampleSD/ sqrt(n - 1))
(lowerBound = sampleMean - width)## [1] 20.01547(upperBound = sampleMean + width)## [1] 20.60901From this we can see that the final 95% confidence interval is: \[(20.0154696, 20.6090104)\]
Question 8:
Call strategy B rolling the dice as follows until getting a Dragon Jack. If D20 > 15 or < 5, reroll all three dice. If D20 is in [5, 15], hold that value of D20 and keep rerolling the D10 and D12. Rolling real dice, I simulated 17 games using strategy B which are described below.
mean(strategyB) = 15.41176
sd(strategyB) = 20.29796
sort(strategyB) = 1 4 5 5 6 7 8 8 9 10 13 13 15 18 20 31 89
Calculate a 95% confidence interval expressed as (LB, UB) for the true mean number of rolls to get a Dragon Jack using strategy B. In addition, justify or critique the confidence interval method you used in a couple short sentences.
sampleMeanB <- 15.41176
sampleSDB <- 20.29796
nB <- 17We proceed in the same manor as the last question. Again I am choosing to use the T confidence interval. This case is a more classic use case for the T as it has not only an unknown SD but also a small \(n\). The normal would be overly confident in this instance.
alpha = 0.05 #95% interval
tValB <- qt(1 - (alpha/2), nB - 1)
widthB <- tValB * (sampleSDB/ sqrt(nB - 1))
(lowerBoundB = sampleMeanB - widthB)## [1] 4.654322(upperBoundB = sampleMeanB + widthB)## [1] 26.1692From this we can see that the final 95% confidence interval is: \[(4.6543218, 26.1691982)\]
Question 9
Referring to the data presented for strategy B in the previous question, calculate a 95% confidence interval expressed as (LB, UB) for the true median number of rolls to get a Dragon Jack using strategy B. Hint: you have not done this in any of our quizzes.
Note: I completely forgot to come back and do this one on the in class. Oops.
A way to get confidence intervals for the median that we learned in Brian’s Probability class was to use bootstrapping to generate a ton of fake results, get their medians and then find the percentiles of those medians. Code inspiration for this comes from this stackoverflow question.
#Load in the results
results = c(1,4,5,5,6,7,8,8,9,10,13,13,15,18,20,31,89)
#Generate bootstraped medians
bootstraped <- replicate(1000, median(sample(results, length(results), replace=T)))
#Find the quantiles
quantile(bootstraped, c(0.025, 0.975))## 2.5% 97.5%
## 6 15Question 10:
When creating the data for strategy B, I played versus strategy A and found that B won 12 out of the 17 games. Calculate a 95% confidence interval expressed as (LB, UB) for the probability of strategy B beating A using a Wilson interval.
k <- 12 #k successes
n <- 17 #n trials
#Get 95% confidence intervals for a binomial sample.
binconf(k, n, method = "wilson")## PointEst Lower Upper
## 0.7058824 0.4686689 0.8672001Question 11:
Suppose I wrote code for strategy B and played B vs A for one million games. Since ties are possible, I’ll calculate the probability of B strictly winning (not losing or tieing). How would the 95% asymptotic Normal (Wald) interval for the probability of B winning compare to the Wilson and Exact intervals in terms of coverage and interval width?
In general the Wald interval results in lower coverage and wider intervals than the Wilson and Exact (see Quiz 5 and Brown et al., 2002 [1]). However, given that we are simulating one million games the intervals will be very close to eachother. We saw in the Quiz 5 that even at 4000 simulations we saw the intervals coverage rapidly converging to eachother, however, that being said the Wald interval was still the worst performing.
Question 12:
We usually think of confidence interval symmetry as an important property, but maybe we could create a better interval by sacrificing it. I honestly don’t know, but I have written this question to explore the idea. We’ve studied the symmetric 95% Exact CI for a proportion defined loosely as the (LB, UB) that satisfies P(X > x | theta = LB) = 0.025 and P(X < x | theta = UB) = 0.025. Consider an asymmetric 95% Exact CI for a proportion defined loosely as the (LB, UB) that satisfies P(X > x | theta = LB) = a1 and P(X < x | theta = UB) = a2, where (a1, a2) are defined as follows.
1. If the sample proportion = 0.5, then (a1, a2) = (0.025, 0.025).
2. If the sample proportion > 0.5, then (a1, a2) = (0.01, 0.04).
3. If the sample proportion < 0.5, then (a1, a2) = (0.04, 0.01).
For the data where I found that strategy B beat A 12 out of 17 games, calculate a 95% confidence interval expressed as (LB, UB) for the probability of strategy B winning using this asymmetric 95% Exact CI.
First we write a function to generate the asymptotic confidence interval as described:
#Function for computing the asymmetric exact interval
asymExact = function(n, k){
if((k/n) > 0.5){
lb = binconf(k, n, alpha = 0.02, method = "exact")[2]
ub = binconf(k, n, alpha = 0.08, method = "exact")[3]
} else if((k/n) < 0.5){
lb = binconf(k, n, alpha = 0.08, method = "exact")[2]
ub = binconf(k, n, alpha = 0.02, method = "exact")[3]
} else {
lb = binconf(k, n, alpha = 0.05, method = "exact")[2]
ub = binconf(k, n, alpha = 0.05, method = "exact")[3]
}
return(c(k/n, lb, ub))
}Note: On the in-class I messed up and wrote n/k as opposed to k/n. Oops.
We then call this for our given data.
(asymmetricInterval <- asymExact(17,12))## [1] 0.7058824 0.3974867 0.8834306From this we can get that the asymmetric interval for the outcome is: \[(0.3974867, 0.8834306)\]
Question 13:
Fix n = 17 and vary theta, i.e. let X ~ Bin(17, theta). Create a plot of coverage rate by theta for theta in (theta, 1) for the following three 95% CI methods (Wilson, Symmetric Exact, and the proposed Asymmetric Exact). Take a fine gradation of theta, e.g. delta <- 0.001; theta <-seq(delta,1-delta,delta);. Describe the relative performance of the methods. State which interval method would you pick for this setting of a small n and unknown theta and explain why you picked that method.
To answer this I am utilizing the method that we used in Quiz 5 to get coverage with a high accuracy. The code is the same with a few small modifications for using the asymmetric interval.
#Function to check coverage.
covered <- function(r,n,p, method = "wilson"){
if(method == "asymmetric"){ #Add a switch to use our custom asymmetric function
ci <- asymExact(n, r)
} else {
ci <- binconf(r,n, method = method) #generate ci
}
return(as.numeric(ci[2] < p & ci[3] > p)) #return 1 if p is in ci, 0 if not.
}
expectation <- function(n, p, method = "wilson"){
possible <- seq(0,n) #Make a vector with the n+1 possible results.
C <- sapply(possible, function(r) covered(r, n, p, method)) #grab coverage results
probs <- sapply(possible, function(r) dbinom(r, n, p)) #generate a vector with the prob of all of results.
return(sum(C * probs))
}
performance_plot <- function(n){
#Make our range of ps
delta <- 0.001;
ps <- seq(delta,1-delta,delta)
#Generate the coverage levels for each interval type.
coverages <- data.frame(p = ps,
wilson = sapply(ps, function(p) expectation(n, p, method = "wilson")),
exact = sapply(ps, function(p) expectation(n, p, method = "exact")),
asymmetric = sapply(ps, function(p) expectation(n, p, method = "asymmetric")))
#Manipulate and plot.
coverages_plot = melt(coverages, id = c("p"))
colnames(coverages_plot) = c("p", "interval", "coverage")
plt <- ggplot(coverages_plot, aes(x = p, y = coverage, group = interval, color = interval)) +
geom_line() + geom_hline(aes(yintercept = 0.95))
return(plt) #return the ggplot object for later customization.
}
performance_plot(17) + ylim(0.88, 1) + labs(title = "Interval comparison with n = 17")
From the plot we can see that the wilson method is generally the worst for coverage. The exact interval as seen before performs consistantly above 95% and is particularly strong at the edges. Our new asymmetric interval outperforms the exact in the middle range (from \(\theta = 0.25 \to 0.75\)). This seems to be because the exact interval was already so good at dealing with the edge cases that meddling with it just messed it up. Although the performance didn’t shift too much.
Personally, I would still stick with the exact interval for this case. I think that it’s higher performance on the edges is more important than the asymmetric’s good performance near the middle. Maybe if the sample theta was near 0.5 it would be worth it, but at this point I would go with the tried and true exact method.
Question 14:
Propose your own Dragon Jacks strategy C. Prove that your strategy beats strategy A and if possible strategy B. Note the word “prove” is left intentionally ambiguous. Make a convincing case.
So first we start off by coding up the game.
Note: I have coded up my strategies in this as well. Ignore these till later in the question.
#Function to simulate the dice. E.g. D10 = dice(10)
dice <- function(n = 10) sample(seq(1,n), 1)
#A function simulating a roll.
#If the value for a given dice is 0, the function will roll that dice.
roll <- function(d10 = 0, d12 = 0, d20 = 0){
d10 <- ifelse(d10 == 0, dice(10), d10)
d12 <- ifelse(d12 == 0, dice(12), d12)
d20 <- ifelse(d20 == 0, dice(20), d20)
return(c(d10, d12, d20))
}
playIt <- function(d10 = 0,d12 = 0,d20 = 0,count = 1, strategy = "a"){
r <- roll(d10,d12,d20) #Roll the dice, holding the specified dice.
# print(r) #uncomment to debug.
if(sum(r) != 21){ #If we didn't win, roll again, put heuristics in this step.
if(strategy == "a"){
#Do nothing.
} else if(strategy == "b"){
#Strategy b
if(r[3] <= 15 && r[3] >= 5) d20 = r[3] #if d20 in [5,15] keep it.
} else if (strategy == "c"){
#Strategy c
#Since the 10 sided is most predicatable, hold if the sum of the other two is in it's S
if((r[2] + r[3]) > 11 && (r[2] + r[3]) < 20){
d12 = r[2]
d20 = r[3]
}
} else {
#Strategy d
if((r[1] + r[3]) > 8 && (r[1] + r[3]) < 21){ #if d10 + d20 is in the range of the d12, keep em'
d10 = r[1]
d12 = 0 #reroll d12
d20 = r[3]
}
if((r[2] + r[3]) > 10 && (r[2] + r[3]) < 21){ #but if d12 + d20 is in the range of the d10 that's even better.
d10 = 0 #reroll d10
d12 = r[2]
d20 = r[3]
}
}
#Roll the dice again.
playIt(d10 = d10,d12 = d12,d20 = d20, count = count + 1, strategy = strategy)
} else { #We've won, return the count.
return(list(count = count, lastRoll = r))
}
}Okay. Now that it’s all in, let’s test it to make sure it works.
Let’s just play a normal un-strategied game.
playIt() ## $count
## [1] 4
##
## $lastRoll
## [1] 5 9 7Exploration of winning rolls:
Since we have a function written that will play dragon jacks, let’s play a bunch of games and record what are the charectaristics of the winning rolls. If obvious patterns jump out in the winning dice hand then we may potentially be able to glean some information on a good strategy for winning faster.
n <- 50000 #Number of simulations to run.
#Play n games and assemble a dataframe of the results
playLots = as.data.frame(t(replicate(n, playIt()$lastRoll)))
names(playLots) = c("d10", "d12", "d20")
summary(playLots) #What do we get from these. ## d10 d12 d20
## Min. : 1.000 Min. : 1.000 Min. : 1.000
## 1st Qu.: 3.000 1st Qu.: 3.000 1st Qu.: 6.000
## Median : 5.000 Median : 6.000 Median : 9.000
## Mean : 5.391 Mean : 6.384 Mean : 9.226
## 3rd Qu.: 8.000 3rd Qu.: 9.000 3rd Qu.:12.000
## Max. :10.000 Max. :12.000 Max. :19.000#Let's plot to really see.
ggplot(data = melt(playLots), aes(factor(variable), value)) +
geom_violin(aes(fill = factor(variable))) +
scale_fill_discrete(name = "Die") +
labs(x = "", y = "Value", title = "Dice Value at Win")## Using as id variables
From this we can see that there really isn’t much modality in the winning value for D10 and D12 dice. The probability on winning on a given value of them is uniform throughout their sample space. This indicates that the main strategy will almost assuredly depend heavily on the D20.
Is there potential information lost in looking simply at the winning values?
Let’s plot a scatter looking at each game as one point in 3-space. Maybe trends will pop out.
#Let's make a 3d scatter plot to see if that will help.
uniqueRolls = ddply(playLots,.(d10,d12,d20),nrow)
names(uniqueRolls) <- c("d10", "d12", "d20", "count")
#Scale the values to a useful range
uniqueRolls$countScalled <- sapply(uniqueRolls$count, FUN = function(X) 3*(X - min(uniqueRolls$count))/diff(range(uniqueRolls$count)))
N = max(uniqueRolls$count)
countRange <- max(uniqueRolls$count) - min(uniqueRolls$count)
colfunc <- colorRampPalette(c("blue", "red"))(countRange)
scatterplot3js(uniqueRolls$d10, uniqueRolls$d12, uniqueRolls$d2, size = uniqueRolls$countScalled,color= "steelblue",renderer="canvas")Alas, this doesn’t actually seem to help us.
Strategies:
I have proposed two new strategies.
Strategy C:
Here I simply check if the sum of the D20 and the D12 are in the range of 11 \(\to\) 20. If it is, then we put the D12 and D20 to the side just roll the D10 to finish out. Once this condition is triggered the average number of rolls to win will be E[D10].
Strategy D:
This strategy simply uses strategy C and takes it one step further. In this case first we check to see if the sum of the D10 and the D20 are in the range of 9 \(\to\) 20, if this is the case we hold both of those dice. Then however, we perform strategy C again, because if we can roll the D10 instead of the D12 we will finish E[D12] - E[D10] rolls earlier on average.
Test the strategies
n <- 100000
naive = replicate(n, playIt()$count)
summary(naive)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 6.00 14.00 20.53 28.00 298.00stratB = replicate(n, playIt(strategy = "b")$count)
summary(stratB)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 5.00 11.00 15.97 21.00 209.00stratC = replicate(n, playIt(strategy = "c")$count)
summary(stratC)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 4.00 8.00 11.22 15.00 103.00stratD = replicate(n, playIt(strategy = "d")$count)
summary(stratD)## Min. 1st Qu. Median Mean 3rd Qu. Max.
## 1.00 4.00 8.00 10.87 15.00 116.00results = data.frame(a = naive, b = stratB, c = stratC, d = stratD)
ggplot(data = melt(results), aes(factor(variable), value)) + geom_boxplot(aes(fill = factor(variable)))## Using as id variables
t.test(stratB, stratD, alternative = "greater")##
## Welch Two Sample t-test
##
## data: stratB and stratD
## t = 87.388, df = 164803.4, p-value < 2.2e-16
## alternative hypothesis: true difference in means is greater than 0
## 95 percent confidence interval:
## 5.001778 Inf
## sample estimates:
## mean of x mean of y
## 15.96724 10.86951So as we can see strategy D is the clear winner with a highly significant t distribution hypothesis test p-value. We chose a t-test because of robustness to non-normal distributions and an unknown population variance.
Bibliography:
- L. D. BROWN, T. T. CAI AND A. DASGUPTA, CONFIDENCE INTERVALS FOR A BINOMIAL PROPORTION AND ASYMPTOTIC EXPANSIONS, The Annals of Statistics, 2002, Vol. 30, No. 1, 160–201