sim<-function(n,nrep){ BestSoFar<-BestOverall<-numeric(n) for(i in 1:nrep){ Order<-sample(n) BestOverall<-BestOverall+(Order==n) BestSoFar<-BestSoFar+(Order==cummax(Order)) } if(nrep==1) rbind(Order=Order,BestOverall=BestOverall,BestSoFar=BestSoFar) else rbind(BestOverall=BestOverall,BestSoFar=BestSoFar,Ratio=BestOverall/BestSoFar, Theoretical=(1:n)/n) } sim(10,1) sim(10,1) sim(10,1) sim(10,1e6)