6

I would like to add a parabola line denoting 95% confidence limits to this coin toss plot using R:

x  <- sample(c(-1,1), 60000, replace = TRUE)
plot.ts(cumsum(x), ylim=c(-250,250))

Here is an example of what I'm looking for:graph

UPDATE: @bill_080's answer is excellent. However I have already calculated 100,000 coin tosses:

str(100ktoss)
num [1:100000] -1 1 1 1 -1 -1 1 -1 -1 -1 ...

and I really want to just add the 95% limit to that plot:toss

plot.ts(cumsum(100ktoss))

It took several hours to calculate my 100K coin tosses and when I try and replicate with @bill_080's code I run out of memory (for 100,000).

FINAL UPDATE: Okay. Last problem. I have a plot of several rounds of cummulative hits, on a single graph with the start of each round clamped at zero (actually 1 or -1 depending on if it was a win or lose).

>str(1.ts)  
Time-Series [1:35] from 1 to 35: 1 2 1 2 3 4 5 4 5 6 ...  
>str(2.ts)  
Time-Series [1:150] from 36 to 185: -1 0 1 0 -1 -2 -1 0 1 2 ...  

I would like to add the same 95% limit to each segment, like thus. Now solved:

@bill_080 Many thanks. This is the the end product:

cum

Bobby
  • 11,419
  • 5
  • 44
  • 69
Frank Zafka
  • 829
  • 9
  • 30

1 Answers1

9

Try this. All loops are for loops, so you can easily add more calculations.

#Set the number of bets and number of trials and % lines
numbet <- 6000 #6000 bets
numtri <- 1000 #Run 1000 trials of the 6000 bets
perlin <- 0.05 #Show the +/- 5% lines on the graph
rantri <- 60 #The 60th trial (just a random trial to be drawn)

#Fill a matrix where the rows are the cumulative bets and the columns are the trials
xcum <- matrix(NA, nrow=numbet, ncol=numtri)
for (i in 1:numtri) {
  x <- sample(c(-1,1), numbet, replace = TRUE)
  xcum[,i] <- cumsum(x)
}

#Plot the trials as transparent lines so you can see the build up
matplot(xcum, type="l", xlab="Number of Bets", ylab="Cumulative Sum", main="Cumulative Results", col=rgb(0.01, 0.01, 0.01, 0.02))
grid()

#Sort the trials of each bet so you can pick out the desired %
xcumsor <- xcum
for (i in 1:numbet) {
  xcumsor[i,] <- xcum[i,order(xcum[i,])]
}

#Draw the upper/lower limit lines and the 50% probability line     
lines(xcumsor[, perlin*numtri], type="l", lwd=2, col=rgb(1, 0.0, 0.0)) #Lower limit
lines(xcumsor[, 0.5*numtri], type="l", lwd=3, col=rgb(0, 1, 0.0)) #50% Line
lines(xcumsor[, (1-perlin)*numtri], type="l", lwd=2, col=rgb(1, 0.0, 0.0)) #Upper limit

#Show one of the trials
lines(xcum[, rantri], type="l", lwd=1, col=rgb(1, 0.8, 0)) #Random trial

#Draw the legend
legend("bottomleft", legend=c("Various Trials", "Single Trial", "50% Probability", "Upper/Lower % Limts"), bg="white", lwd=c(1, 1, 3, 2), col=c("darkgray", "orange", "green", "red"))

enter image description here

Edit 1 ==========================================================

If you're just trying to draw the +/- 5% lines, it's just a square root function. Here's the code:

#Set the bet sequence and the % lines
betseq <- 1:100000 #1 to 100,000 bets
perlin <- 0.05 #Show the +/- 5% lines on the graph

#Calculate the Upper and Lower limits using perlin
#qnorm() gives the multiplier for the square root
upplim <- qnorm(1-perlin)*sqrt(betseq)
lowlim <- qnorm(perlin)*sqrt(betseq)

#Get the range for y
yran <- range(upplim, lowlim)

#Plot the upper and lower limit lines
plot(betseq, upplim, ylim=yran, type="l", xlab="", ylab="")
lines(betseq, lowlim)

enter image description here

Edit 2 ==================================================

To add the parabolas at the right locations, it is probably easier if you define a function. Keep in mind that because the new function (dralim) uses lines, the plot has to exist before you call dralim. Using some of the same variables as the code in Edit 1:

#Set the bet sequence and the % lines
betseq <- 0:700 #0 to 700 bets
perlin <- 0.05 #Show the +/- 5% lines on the graph

#Define a function that plots the upper and lower % limit lines
dralim <- function(stax, endx, perlin) {
  lines(stax:endx, qnorm(1-perlin)*sqrt((stax:endx)-stax))
  lines(stax:endx, qnorm(perlin)*sqrt((stax:endx)-stax))
}

#Build the plot area and draw the vertical dashed lines
plot(betseq, rep(0, length(betseq)), type="l", ylim=c(-50, 50), main="", xlab="Trial Number", ylab="Cumulative Hits")
abline(h=0)
abline(v=35, lty="dashed") #Seg 1
abline(v=185, lty="dashed") #Seg 2
abline(v=385, lty="dashed") #Seg 3
abline(v=485, lty="dashed") #Seg 4
abline(v=585, lty="dashed") #Seg 5

#Draw the % limit lines that correspond to the vertical dashed lines by calling the
#new function dralim.
dralim(0, 35, perlin) #Seg 1
dralim(36, 185, perlin) #Seg 2
dralim(186, 385, perlin) #Seg 3
dralim(386, 485, perlin) #Seg 4
dralim(486, 585, perlin) #Seg 5
dralim(586, 701, perlin) #Seg 6

enter image description here

bill_080
  • 4,692
  • 1
  • 23
  • 30
  • @bill_080 Really, really nice. If you check out my profile you will see I have a thing for coin toss plots and yours is a beauty. The only thing it is maybe a bit too complicated for me. :) – Frank Zafka May 27 '11 at 20:06
  • 1
    @RSoul, I added another plot with the associated code. It shows the upper/lower limits for the 5% lines. It's just the square root of the bet number times a multiplier for the % probability that you want. In the code, the `qnorm()` function gives you the multiplier. For 5%, `qnorm(0.05)` gives -1.644 and `qnorm(0.95)` gives 1.644. – bill_080 May 27 '11 at 20:56
  • @bill_080 That is excellent. Exactly what I wanted. This is why I have acknowledged the users of StackOverflow.com in my thesis. Many thanks. I may use your more interesting plot. I've yet to decide. – Frank Zafka May 27 '11 at 21:04
  • @bill_080 I can calculate the separate parabola lines for each round, but not put them in the right place on the x-axis. – Frank Zafka May 28 '11 at 16:07
  • 1
    @RSoul, I updated my answer with Edit 2. In the future, it's probably better if you ask a new question for your Final Update. You can still reference the original question. – bill_080 May 28 '11 at 20:58
  • That is absolutely brilliant. It has all worked perfectly and you've given me lots to work with. I wish I could donate more points or something. Many thanks @bill_080. – Frank Zafka May 28 '11 at 21:43