If there is a God, he’s a great mathematician (Paul Dirac)
Imagine you toss a coin 12 times and you count how many heads and tails you are obtaining after each throwing (the coin is equilibrated so the probability of head or tail is the same). At some point, it can happen that number of heads and number of tails are the same. For example, if you obtain the sequence T-H-T-T-H-T-H-H-T-T-H-H, after the second throwing, number of heads is equal to number of tails (and both equal to one). It happens again after the 8th throwing and after last one. In this example, the last throwing where equallity occurs is the number 12. Obviously, equallity can only be observed in even throwings.
If you repeat the experiment 10.000 times you will find something like this if you draw the relative frequency of the last throwing where cumulated number of heads is equal to the one of tails:
From my point of view there are three amazing things in this plot:
- It is symmetrical, so
prob(n)=prob(12-n) - The least likely throwing to obtain the last equality is the central one.
- As a corollary, the most likely is not obtaining any equality (number of heads never are the same than number of tails) or obtaining last equality in the last throwing: two extremely different scenarios with the same chances to be observed.
Behind the simplicity of tossing coins there is a beautiful universe of mathematical surprises.
library(dplyr)
library(ggplot2)
library(scales)
tosses=12
iter=10000
results=data.frame(nmax=numeric(0), count=numeric(0), iter=numeric(0))
tmp=data.frame(nmax=numeric(0))
for (j in 1:iter)
{
data.frame(x=sample(c(-1,1), size=tosses, replace=TRUE)) %>%
add_rownames(var = "n") %>%
mutate(cumsum = cumsum(x)) %>% filter(cumsum==0) %>%
summarize(nmax=max(as.numeric(n))) %>% rbind(tmp)->tmp
}
tmp %>%
group_by(nmax) %>%
summarize(count=n()) %>%
mutate(nmax=ifelse(is.finite(nmax), nmax, 0), iter=iter) %>%
rbind(results)->results
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=20),
legend.key = element_blank(),
plot.title = element_text(size = 30)
)
ggplot(results, aes(x=nmax, y=count/iter)) +
geom_line(size=2, color="green4")+
geom_point(size=8, fill="green4", colour="darkolivegreen1",pch=21)+
scale_x_continuous(breaks = seq(0, tosses, by=2))+
scale_y_continuous(labels=percent, limits=c(0, .25))+
labs(title="What happens when you toss a coin 12 times?",
x="Last throwing where cumulated #tails = #heads",
y="Probability (estimated)")+opts
I think,that’s because the formulae of permutation and combination -> C(n,m)=C(n,n-m)
You didn’t mention that this is also an expensive way to draw a smiley 🙂
Jajajaja … It’s true!
For those who are unfamiliar with it, check out the arcsine law for random walks. For example a good discussion is in Feller’s Introduction to Probability Theory and its Applications, Vol.1
I know it’s petty, but should catch the attempt to process -Inf values inside “tmp” somewhere in that code.
This can be done with straight R code. In the code sample that follows, the for {…} structure is eliminated in favor of apply(…), which I think runs faster but I haven’t benchmarked it.
tosses=12
iter=10000
sampl=sample(c(-1,1),size=tosses*iter,replace=TRUE)
sampl.m=matrix(sampl,nrow=tosses,ncol=iter)
results=apply(sampl.m,2,function(x) max((cumsum(x)==0)*(1:12)))
plot(table(results)/iter,xlab=”Last Throw…”,ylab=”Probability…”,main=”What Happens…?”)
In the last statement I’ve included a basic plot statement but the author’s (very sophisticated!) ggplot commands could be used as well.
Also, the third line could be combined with the fourth line and still maintain readability but i figured the increment in clarity was worth its inclusion.
I know some people love pipes–a holdover from the old Unix days, I suppose. But for my money simple programming with a minimum of piping, recursion or jumping is greatly to be desired.
Finally, the exact probability distribution is given by,
p(k)=n(0,k)*n(0,12-k) / 2^12
Where n(0,k) is the equal to the number of random walk paths that make a return to the origin at epoch k. And n(0,12-k) is equal to the number of random walks that do not return to the origin by epoch 12-k. Again, I refer the reader to Feller for the mathematical details.
The numerical values of the pdf are,
0 2 4 6 8 10 12
0.225586, 0.123047, 0.102539, 0.0976563, 0.102539, 0.123047, 0.225586
Note that these can be approximated by the arcsine distribution.
Thank you for the code. Very clear!
Thanks, I really hope it is! I think R’s basic package set does a great job. There are lots of specialized packages and I’ve tried them. Mostly (because of advanced age) I forget they are there and write my own code.
Generally I find the apply(…) functions to be faster and more clear than f{…} structures. I read somewhere the apply(…) functions are optimized for array operations or some such. I started programming in the days of 80-column Hollerith cards. In those days you knew the cost of every machine function and resource.
Also I was a little skimpy on the math details because I don’t know how to write mathematical expressions for this type of blog. Can you do it?
Also, I noticed you use ggplot and not ggplot2. Not to imply it is wrong but why don’t you use ggplot2?
Is there a random-walk func in some common R-package?
I came across a randomwalk package by Bruno Le Floch. Don’t know anything about it but looks like you can do various 2D walks with and without constraints. But I recommend to all William Feller’s Introduction to Probability Theory and its Applications. Volume 1 deals mainly with discrete processes and Volume 2 with continuous. Note that Vol. 2 employs almost no measure theory yet manages to cover lot of important topics. In both volumes there are lots of interesting and workable problems.
For those who don’t have Feller, that n(0,x) is equal to binomial( n, n/2) *2^(-n) ; the probability of finding the last such x involves a sum over all “not” cases.
If I were on a desert island with just a laptop, internet connection and my choice of five books of probability and statistics the books would be,
1) William Feller, Introduction to Probability and it’s Applications Vol. I&II (counts as one!)
2) Harald Cramer, Mathematical Methods of Statistics
3) S. Wilks, Mathematical Statistics
4) Herbert A. David, Order Statistics
5) Maurice Kendall, Rank Correlation Methods
That’s it! All the rest can be done.
But Douglas, if you had an internet connection, you could download most if not all of those 🙂
I know you can download Cramer’s book. How about the rest? In any event, I find it difficult to read on my computer. I need a book in hand.
I haven’t looked for all of them, but both volumes of Feller were easy to find (albeit as scanned images)