# File "Rcoins".

seqlength = 100
numreps = 1000

listheads = listchanges = listlongest = listthree = NULL

for (i in 1:numreps) {

    # Generate a random sample of Bernoulli(1/2) random variables.
    seq = rep(0,seqlength)
    for (j in 1:seqlength) {
        if (runif(1) < 0.5) {
	    seq[j] = 1
	}
    }
    # OR, INSTEAD COULD SIMPLY SAY:  seq = rbinom(seqlength,1,0.5)

    # Compute numheads.
    numheads =  sum(seq)
    listheads = c(listheads, numheads)

    # Compute numchanges
    numchanges = 0
    for (j in 2:seqlength) {
	if (seq[j] != seq[j-1]) {
	    numchanges = numchanges + 1
	}
    }
    listchanges = c(listchanges, numchanges)

    # Compute numlongest.
    numlongest = 0
    currun = 0  # count the length of the current run
    curvalue = 2  # special value, neither heads nor tails
    for (j in 1:seqlength) {
        if (seq[j] == curvalue) {
	    # current run continues!
	    currun = currun + 1
	    # check if we've set a new record
	    if (currun > numlongest) {
		numlongest = currun
	    }
	} else {
	    # current run ends; start new run!
	    curvalue = seq[j]
	    currun = 1
	}
    }
    listlongest = c(listlongest, numlongest)

    # Compute numthree
    numthree = 0
    for (j in 1:(seqlength-2)) {
	if ( (seq[j] == seq[j+1]) && (seq[j+1] == seq[j+2]) ) {
	    numthree = numthree + 1
	}
    }
    listthree = c(listthree, numthree)

}

cat("\n")
cat("heads:    mean = ", mean(listheads), "   sd = ", sd(listheads), "\n\n")
cat("changes:    mean = ", mean(listchanges), "   sd = ", sd(listchanges), "\n\n")
cat("longest:    mean = ", mean(listlongest), "   sd = ", sd(listlongest), "\n\n")
cat("three:    mean = ", mean(listthree), "   sd = ", sd(listthree), "\n\n")

hist(listheads, breaks=0:100)
hist(listchanges, breaks=0:100)
hist(listlongest, breaks=0:max(listlongest))
hist(listthree, breaks=0:max(listthree))