= R code computing the replication probabilities for the example of 92 stimuli randomly drawn with replacement into sequences of length 6 =

{{{
#
# From Daniel Molinari July 2015
#

# S.x.y = y+1 instances of x replicates
# 0 replicates   abcdef

  S.0 <- choose(92,6)*factorial(6)


# 1-0 replicates aabcde

  S.1.0 <- choose(6,2)*92*91*90*89*88 +
           choose(6,2)*choose(4,2)*92*91*90*89 +
           choose(6,2)*choose(4,2)*choose(2,2)*92*91*90  

# 1.1 PW addition using Laurence S formula aabbcd 

 S.1.1 <- choose(92,2)*choose(90, 2)*15*6*2 


# 1-2 replicates aabbbc
  
  S.1.2 <- choose(6,2)*choose(4,3)*92*91*90


# 1-3 replicates aabbbb
  
  S.1.3 <- choose(6,2)*choose(4,4)*92*91 


# 2-0 replicates aaabcd
  
  S.2.0 <- choose(6,3)*92*91*90 

# PW ADDING IN LAURENCE S FORMULA 2-1 replicates aaabbb

  S.2.1 <- factorial(2)*choose(92,2)*choose(6,3)/92^6 

# 2-2 replicates aaabbb
  
  S.2.2 <- choose(6,3)*choose(3,3)*92*91 


# 3-0 replicates aaaabc
  
  S.3.0 <- choose(6,4)*92*91*90

# PW addition 
# Yes, the order counts.  The correct computation is 
# Sequences( 3 instances 1 replicate) =  3 ! (92 Choose 3)*(6 choose 2)*(4 choose 2)   

  S.3.1 <-  factorial(3)*choose(92,3)*choose(6,2)*choose(4,2)   


# 4-0 replicates aaaaab
  
  S.4.0 <- choose(6,5)*92*91 


# 5-0 replicates aaaaaa
  
  S.5.0 <- choose(6,6)*92  


# Compute cumulative probability

  S <- S.0 + 
       S.1.0 + S.1.2 + S.1.3 +
       S.2.0 + S.2.2 +
       S.3.0 +
       S.4.0 +
       S.5.0

   P  <- S / (92^6)

  S <- S.0 + 
       S.1.0 + S.1.1 + S.1.2 + S.1.3 +
       S.2.0 + S.2.1 + S.2.2 +
       S.3.0 + S.3.1
       S.4.0 +
       S.5.0

   P1 <- S / (92^6)
}}}

yields

{{{
  P    # 1.002883
  P1   # 1.007972  
}}}