y=sort(c(3,1,2,5))
theta=c(.62,.72,by=.01)
par(mfrow=c(2,1))
for(p in theta){
  Nreps=1000000
  StrkLengths=matrix(nrow=Nreps, ncol=4)
  for(i in 1:Nreps){
    StrkID=0
    samp.prev=sample(0:1,1,prob=c(1-p,p))
    if(samp.prev==1){StrkID=StrkID+1}
    rec=cbind(StrkID,samp.prev)
    colnames(rec)=c("StreakID","Result")

    while(StrkID<5){
      samp.new=sample(0:1,1,prob=c(1-p,p))
      if(samp.prev==0 & samp.new==1){StrkID=StrkID+1}
      if(StrkID<5){
        rec=rbind(rec,cbind(StrkID,samp.new))
        samp.prev=samp.new
      }
    }
    StrkLengths[i,]=table(rec[-which(rec[,"Result"]==0),1])
  }
  StrkLengths=t(apply(StrkLengths,1,sort))
  matches=which(apply(StrkLengths,1,toString)==toString(y))
  hist(StrkLengths, breaks=seq(0,max(StrkLengths),by=1), col="Grey",
       right=F, main=paste0("theta=",p, " Nmatch=", length(matches)))
  t=1:max(StrkLengths)
  lines(t+0.5,4*Nreps*(p^(t-1)-p^(t)), type="b", lwd=3)
}