Say you have M rounds. We start by calculating how big is the probability of being at zero at round M. As you said, you have three results for each round: -1, 0 and 1 with 0.25, 0.5 and 0.25 probabilities respectively. In order to arrive at zero, one must obtain an equal number (call it k) of +1 and -1, while the remaining can be zeroes. Of course, k can range from 0 to M/2 (/ being the integer division). So you have M-2*k zeroes, k +1s and k -1s. You can arrange your zeroes whenever you want in your M rounds and next you arrange k +1s in the 2*k remaining spots. You finally sum for k ranging from zero to M/2. In formula:
P_0(M) = Sum_{K=0}^{M/2} C(M,M-2K) * (1/2)^(M-2K) * C(2K,K) * (1/4)^(2K)
With a little algebra, and by noting that the ratio P_0(M+1)/P_0(M) is equal to 1-1/(2M), we arrive to:
P_0(M) = Prod_{i=1}^{M} (1-1/(2M)) = Gamma(M+1/2)/sqrt(pi)/Gamma(M+1)
The last step has been obtained with the help of Wolfram Alpha. Gamma is the
gamma function.
We now consider the expected number of total times we are at zero in M rounds. Call this quantity E_T[M]. We can exploit linearity of the expected value without having to deal with mutual dependencies:
E_T[M] = Sum_{i=1}^M E_0[i] = Sum_{i=0}^M P_0(i)
where E_0[i] is the expected number of times the i-th round is zero, which is obviously equal to P_0(i). Luckily for us, even the sum above can be simplified using the properties of the gamma function (and Wolfram Alpha of course), leading to:
E_T[M] = 2 * (M+1) * Gamma(M+3/2) / sqrt(pi) / Gamma(M+2) - 1
Now to lead changes. It must be noted that:
1 - the first non zero value always represents a lead change;
2 - for every zero outside an eventual initial sequence of zeroes and not in the last spot, the probability of having a lead change is 1/4.
In formula:
E_LC[M] = E_L0[M] + 1/4*(E_T(M-1)-E_L0[M])
where E_L0[M] is the expected number of leading zeroes and E_LC[M] is the expected number of lead changes in M rounds. Here is an implementation in the R language:
Code:
expectedzeroes<-function(M) {
exp(log(2)+log(M+1)+lgamma(M+3/2)-0.5*log(pi)-lgamma(M+2))-1
}
expectedLeadingZeros<-function(M) {
2*(1-(1/2)^(M+1))-1
}
expectedChanges<-function(M) {
3/4*expectedLeadingZeros(M) + 1/4*expectedzeroes(M-1)
}
Also, some routine for a simulation:
Code:
simulationLeadChanges<-function(n,nsim=100000) {
simula<-replicate(nsim,sample(-1:1,n,prob=c(1,2,1),replace=TRUE))
simula2<-sign(apply(simula,2,cumsum))
res<-apply(simula2,2,function(x) length(rle(x[x!=0])$values))
c(mean(res),sd(res))
}
And a test (it will took awhile for every simulation to run):
Code:
n<-seq(10,1000,by=10)
calc<-expectedChanges(n)
fromsim<-matrix(0,nrow=length(n),ncol=2)
for (i in seq_len(nrow(fromsim))) {
fromsim[i,]<-simulationLeadChanges(n[i],nsim=10000)
cat("Done n:",n[i],"\n")
}
plot(n,calc,ty="l",lwd=2)
points(n,fromsim[,1],col="blue")
