George Lucas Movie Experiment

#
#see http://stats.stackexchange.com/a/35707/9013 for more details
#

#total viewers of both "Episode A" and "Episode B" movies
totalViewers <- 20000

#movie visit simulations, each item representing the length of all "Episode A" viewer chains (each chain gets interrupted by an "Episode B" viewer at some point, we sum up the total length of all "Episode A" chains)
simulationCount <- 3000
simulations <- vector(length=simulationCount)

#run the simulator
for(i in 1:simulationCount)
{
    #simulate a movie visit: generate a vector of {-1,1},  "1" means the viewer watched "Episode A" and "-1" means "Episode B" (50% probability for each)
    movieVisit <- 2*rbinom(totalViewers, 1, 0.5) - 1

    #calculate the "Episode A" viewer chains (each positive number in the vector is the length of the "Episode A" viewers chain until it's interrupted by an "Episode B" viewer)
    episodeAViewers <- cumsum(movieVisit)

    #remember the total length of all "Episode A" viewer chains for this simulation
    simulations[i] <- sum(episodeAViewers>=0)/totalViewers
}

#plot the histogram
hist(simulations, freq=FALSE, xlab="Proportion of time one movie is the leader", main="George Lucas experiment")
curve(dbeta(x, 1/2, 1/2), 0, 1, col=2, add=TRUE)