sim <- function(n){ # 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 } counter=0 # play counter
# simulation of number of winners among n players NW <- function(m){ # m:players -> (winners,junkens) till any winner j=1 nw = sample(0:(m-1),1,prob=p[[m]]) # number of winners while(nw==0){ # while no winner,repeats j=j+1 nw=sample(0:(m-1),1,prob=p[[m]]) } c(winner=nw,junkens=j) # (number of winners, total plays) } wj=NW(n) if(wj[1]==1) return(wj[2]) # single winner at initial series counter=wj[2] while(wj[1]!=1){ # repeats till single winner determined wj=NW(wj[1]) counter=counter+wj[2] } return(counter) }
