sim <- function(n=10){ # n : number of players p=vector('list',length=n) # probability list p[[1]]=1 # p[[m]][[i]] probability of m players to i winners for(m in 2:n){ k=1:(m-1) p[[m]]=c(3^m-sum(3*choose(m,k)),3*choose(m,k))/3^m }
# simulation of number of winners among n players NW <- function(n,i=1){ # n:players, i:total plays j=1 nw= sample(0:(n-1),1,prob=p[[n]]) # number of winners while(nw==0){ # while no winner,repeats j=j+1 nw=sample(0:(n-1),1,prob=p[[n]]) } c(nw,i-1+j) # (number of winners, total plays) }
wj=NW(n,1) if(wj[1]==1) return(1) # single winner at initial play while(wj[1]!=1){ # repeats till single winner determined wj=NW(wj[1],wj[2]) } return(wj[2]) } j10=mean(replicate(1e6,sim(10))) j10
