Brain Teaser: cómo generar 7 enteros con la misma probabilidad usando una moneda sesgada que tiene una pr (head) = p?

Question

Esta es una pregunta que me he encontrado en Glassdoor: ¿Cómo se generan 7 enteros con igual probabilidad mediante una moneda que tiene un $\mathbb{Pr}(\text{Head}) = p$?

Básicamente, usted tiene una moneda que puede o no puede ser justo, y este es el único número aleatorio de generación de proceso tiene, por lo que vienen con generador de números aleatorios que las salidas de los enteros de 1 a 7, donde la probabilidad de obtener cada uno de estos enteros es 1/7.

La eficiencia de los datos que genera el proceso de la materia.

Answer 1

Voltear la moneda dos veces. Si aterriza HH o TT, ignorarlo y darle la vuelta de nuevo dos veces.

Ahora, la moneda tiene igual probabilidad de entrada HT o TH. Si surge HT, llamar a este H1. Si surge TH, llamar a este T1.

Mantenga la obtención H1 o T1 hasta que tenga tres en una fila. Estos tres resultados dan un número, según la tabla siguiente:

H1 H1 H1 -> 1
H1 H1 T1 -> 2
H1 T1 H1 -> 3
H1 T1 T1 -> 4
T1 H1 H1 -> 5
T1 H1 T1 -> 6
T1 T1 H1 -> 7
T1 T1 T1 -> [Throw out all results so far and repeat]

Yo sostengo que esto funciona perfectamente bien, aunque tiene un gran desperdicio de lanza en el proceso!

Answer 2

Suponga que $p \in (0,1)$.

Paso 1:. Lanzamiento de la moneda de 5 veces.

Si el resultado es

$(H, H, H, T, T)$, regresar $1$ y parada.

$(H, H, T, T, H)$, regresar $2$ y parada.

$(H, T, T, H, H)$, regresar $3$ y parada.

$(T, T, H, H, H)$, regresar $4$ y parada.

$(T, H, H, H, T)$, regresar $5$ y parada.

$(H, H, T, H, T)$, regresar $6$ y parada.

$(H, T, H, T, H)$, regresar $7$ y parada.

Paso 2:. Si el resultado es ninguno de los de arriba, repita el Paso 1.

Tenga en cuenta que independientemente del valor de $p \in (0,1)$, cada uno de los siete resultados antes mencionados tiene probabilidad de $q = p^3(1-p)^2$, y el número esperado de lanzar una moneda es $\displaystyle \frac{5}{7q}$. El tosser no es necesario conocer el valor de $p$ (a excepción de la $p\neq 0$$p\neq 1$); se garantiza que los siete enteros tienen la misma probabilidad de ser devuelto por el experimento cuando se termina (y es garantizado para terminar con una probabilidad de $1$).

Answer 3

La generalización de un caso descrito por Dilip Sarwate

Algunos de los métodos descritos en las otras respuestas el uso de un esquema en el que usted lanza una secuencia de $n$ monedas en un "turno" y dependiendo del resultado que elija un número entre 1 o 7 o descartar la vuelta y tirar de nuevo.

El truco es encontrar en la expansión de las posibilidades de un múltiplo de 7 resultados con la misma probabilidad de $p^k(1-p)^{n-k}$ y coinciden con los uno contra el otro.

Debido a que el número total de resultados no es un múltiplo de 7, tenemos unos resultados que no se puede asignar a un número, y tienen una cierta probabilidad de que tenemos que descartar los resultados y empezar de nuevo.

---

El caso de uso 7 monedas de lanzamientos por turno

Intuitivamente, se podría decir que rodar los dados siete veces, sería muy interesante. Ya sólo tenemos que tirar 2 de la $2^7$ posibilidades. Es decir, el 7 veces los jefes y los tiempos 0 cabezas.

Para todos los otros $2^7-2$ posibilidades hay siempre un múltiplo de 7 casos con el mismo número de cabezas. Es decir, 7 casos con 1 cabezas, 21 casos con 2 cabezas, 35 casos con 3 cabezas, 35 casos con 4 cabezas, 21 casos con 5 cabezas, y 7 casos con 6 cabezas.

Así que si usted calcular el número (descartando 0 cabezas y 7 cabezas) $$X = \sum_{k=1}^{7} (k-1) \cdot C_k $$

with $C_k$ Bernoulli distributed variables (value 0 or 1), then X modulo 7 is a uniform variable with seven possible results.

---

Comparing different number of coin flips per turn

The question remains what the optimal number of rolls per turn would be. Rolling more dices per turn cost you more, but you reduce the probability to have to roll again.

The image below shows a manual computations for the first few numbers of coin flips per turn. (possibly there might be an analytical solution, but I believe it is safe to say that a system with 7 coin flips provides the best method regarding the expectation value for the necessary number of coin flips)

# plot an empty canvas
plot(-100,-100,
     xlab="flips per turn",
     ylab="E(total flips)",
     ylim=c(7,400),xlim=c(0,20),log="y")
title("expectation value for total number of coin flips
(number of turns times flips per turn)")

# loop 1
# different values p from fair to very unfair 
# since this is symmetric only from 0 to 0.5 is necessary 

# loop 2
# different values for number of flips per turn
# we can only use a multiple of 7 to assign 
#   so the modulus will have to be discarded
#   from this we can calculate the probability that the turn succeeds
#   the expected number of flips is 
#       the flips per turn 
#             divided by 
#       the probability for the turn to succeed 

for (p in c(0.5,0.2,0.1,0.05)) {
  Ecoins <- rep(0,16)
  for (dr in (5:20)){
    Pdiscards = 0
    for (i in c(0:dr)) { 
      Pdiscards = Pdiscards + p^(i)*(1-p)^(dr-i) * (choose(dr,i) %% 7)
    }
    Ecoins[dr-4] = dr/(1-Pdiscards)
  }
  lines(5:20, Ecoins)
  points(5:20, Ecoins, pch=21, col="black", bg="white", cex=0.5)
  text(5, Ecoins[1], paste0("p = ",p), pos=2)
}

---

Using an early stopping rule

note: the calculations below, for the expectation value of number of flips, are for a fair coin $p=0.5$, it would become a mess to do this for different $p$, but the principle remains the same (although different book-keeping of the cases is needed)

We should be able to choose the cases (instead of the formula for $$ X) de tal modo que no podría ser capaz de detener el anterior.

• Con 5 coin flips tenemos para las seis posibles diferentes desordenada conjuntos de cabezas y colas:

  1+5+10+10+5+1 conjuntos ordenados

  Y podemos usar los grupos con diez casos (que es el grupo con 2 cabezas o el grupo con 2 colas) para elegir (con igual probabilidad) de un número. Esto ocurre en 14 de los 2^5=32 casos. Esto nos deja con:

  1+5+3+3+5+1 conjuntos ordenados

• Con un extra (6-t) tirón de la moneda tenemos para los siete posibles diferentes desordenada conjuntos de cabezas y colas:

  1+6+8+6+8+6+1 conjuntos ordenados

  Y podemos usar los grupos con ocho casos (que es el grupo con 3 cabezas o el grupo con 3 colas) para elegir (con igual probabilidad) de un número. Esto ocurre en 14 de 2*(2^5-14)=36 casos. Esto nos deja con:

  1+6+1+6+1+6+1 conjuntos ordenados

• Con otro (7-th) extra tirón de la moneda tenemos a los ocho posibles diferentes desordenada conjuntos de cabezas y colas:

  1+7+7+7+7+7+7+1 conjuntos ordenados

  Y podemos usar los grupos con siete casos (todos excepto el de todas las colas y todos los jefes de los casos) a elegir (con igual probabilidad) de un número. Esto ocurre en 42 de los 44 casos. Esto nos deja con:

  1+0+0+0+0+0+0+1 conjuntos ordenados

  (podemos continuar con esto, pero solo en la 49 ª paso hace esto nos da una ventaja)

Así que la probabilidad de seleccionar un número

• a las 5 de la volteretas en el es $\frac{14}{32} = \frac{7}{16}$
• a las 6 de la volteretas en el es $\frac{9}{16}\frac{14}{36} = \frac{7}{32}$
• a las 7 de la volteretas en el es $\frac{11}{32}\frac{42}{44} = \frac{231}{704}$
• no en 7 volteretas en el es $1-\frac{7}{16}-\frac{7}{32}-\frac{231}{704} = \frac{2}{2^7}$

Esto hace que la expectativa de valor para el número de lanzamientos en un turno, con la condición de que no hay éxito y p=0.5:

$$5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704} = 5.796875 $$

The expectation value for the total number of flips (until there is a success), conditional that p=0.5, becomes:

$$\left(5 \cdot \frac{7}{16}+ 6 \cdot \frac{7}{32} + 7 \cdot \frac{231}{704}\right) \frac{2^7}{2^7-2} = \frac{53}{9} = 5.88889 $$

---

The answer by NcAdams uses a variation of this stopping-rule strategy (each time come up with two new coin flips) but is not optimally selecting out all the flips.

The answer by Clid might be similar as well although there might be an uneven selection rule that each two coin flips a number might be chosen but not necessarily with equal probability (a discrepancy which is being repaired during later coin flips)

---

Comparison with other methods

Other methods using a similar principle are the one by NcAdams and AdamO.

The principle is: A decision for a number between 1 and 7 is made after a certain number of heads and tails. After an $x$ number of flips, for each decision that leads to a number $i$ there is a similar, equally probable, decision that leads to a number $j$ (the same number of heads and tails but just in a different order). Some series of heads and tails can lead to a decision to start over.

For such type of methods the one placed here is the most efficient because it makes the decisions as early as possible (as soon as there is a possibility for 7 equal probability sequences of heads and tails, after the $x$-th flip, we can use them to make a decision on a number and we do not need to flip further if we encounter one of those cases).

This is demonstrated by the image and simulation below:

#### mathematical part #####
set.seed(1)


#plotting this method
p <- seq(0.001,0.999,0.001)
tot <- (5*7*(p^2*(1-p)^3+p^3*(1-p)^2)+
       6*7*(p^2*(1-p)^4+p^4*(1-p)^2)+
       7*7*(p^1*(1-p)^6+p^2*(1-p)^5+p^3*(1-p)^4+p^4*(1-p)^3+p^5*(1-p)^2+p^6*(1-p)^1)+
        7*1*(0+p^7+(1-p)^7) )/
             (1-p^7-(1-p)^7)
plot(p,tot,type="l",log="y",
     xlab="p",
     ylab="expactation value number of flips"
     )

#plotting method by AdamO
tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
          (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
          (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
          35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
          1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
           *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
           5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
           15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot,col=2,lty=2)

#plotting method by NcAdam
lines(p,3*8/7/(p*(1-p)),col=3,lty=2)

legend(0.2,500,
       c("this method calculation","AdamO","NcAdams","this method simulation"),
       lty=c(1,2,2,0),pch=c(NA,NA,NA,1),col=c(1,2,3,1))


##### simulation part ######

#creating decision table
mat<-matrix(as.numeric(intToBits(c(0:(2^5-1)))),2^5,byrow=1)[,c(1:12)]
colnames(mat) <- c("b1","b2","b3","b4","b5","b6","b7","sum5","sum6","sum7","decision","exit")

# first 5 rolls
mat[,8] <- sapply(c(1:2^5), FUN = function(x) {sum(mat[x,1:5])})

mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,8]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],12] = rep(5,7) # we can stop for 7 cases with 3 heads
mat[which((mat[,8]==3)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 6th roll
mat <- rbind(mat,mat)
mat[c(33:64),6] <- rep(1,32)
mat[,9] <- sapply(c(1:2^6), FUN = function(x) {sum(mat[x,1:6])})

mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 2 heads
mat[which((mat[,9]==2)&(mat[,11]==0))[1:7],11] = c(1:7)   
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],12] = rep(6,7) # we can stop for 7 cases with 4 heads
mat[which((mat[,9]==4)&(mat[,11]==0))[1:7],11] = c(1:7)    

# extra 7th roll
mat <- rbind(mat,mat)
mat[c(65:128),7] <- rep(1,64)
mat[,10] <- sapply(c(1:2^7), FUN = function(x) {sum(mat[x,1:7])})

for (i in 1:6) {
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],12] = rep(7,7) # we can stop for 7 cases with i heads
  mat[which((mat[,10]==i)&(mat[,11]==0))[1:7],11] = c(1:7)   
}


mat[1,12] = 7           # when we did not have succes we still need to count the 7 coin tosses
mat[2^7,12] = 7


draws = rep(0,100)
num = rep(0,100)
# plotting simulation
for (p in seq(0.05,0.95,0.05)) {
  n <- rep(0,1000)
  for (i in 1:1000) {
    coinflips <- rbinom(7,1,p)  # draw seven numbers
    I <- mat[,1:7]-matrix(rep(coinflips,2^7),2^7,byrow=1) == rep(0,7)                      # compare with the table
    Imatch = I[,1]*I[,2]*I[,3]*I[,4]*I[,5]*I[,6]*I[,7]        # compare with the table 
      draws[i] <- mat[which(Imatch==1),11]                 # result which number
      num[i]   <- mat[which(Imatch==1),12]                 # result how long it took
  }
  Nturn <- mean(num)                   #how many flips we made
  Sturn <- (1000-sum(draws==0))/1000   #how many numbers we got (relatively)
  points(p,Nturn/Sturn)
}

another image which is scaled by $p*(1-p)$ para una mejor comparación:

Answer 4

Divida un cuadro en siete regiones de área igual, cada una etiquetada con un número entero. Tirar la moneda en la caja de tal manera que tenga la misma probabilidad de aterrizar en cada región.

Esto es solo una broma. Es esencialmente el mismo procedimiento que estimar$\pi$ usando un procedimiento físico de Monte Carlo, como dejar caer granos de arroz sobre papel con un círculo dibujado sobre él.

Esta es una de las únicas respuestas que funciona para el caso de$p = 1$ o$p=0$.

Answer 5

EDIT: se basa en otros comentarios.

He aquí una idea interesante:

establecer la lista de {1,2,3,4,5,6,7}. Tirar la moneda para cada elemento de la lista de forma secuencial. Si cae de cabeza de lado para un elemento en particular, eliminar el número de la lista. Si todos los números de una determinada iteración de la lista se quitan, repetir el muestreo. Así hasta que sólo un número restos.

drop.one <- function(x, p) {
  drop <- runif(length(x)) < p
  if (all(drop))
    return(x)
  return(x[!drop])
}

sample.recur <- function(x, p) {
  if (length(x) > 1)
    return(sample.recur(drop.one(x, p), p))
  return(x)
}

# x <- c(1:7,7:1)
x <- 1:7
p <- 0.01

out <- replicate(1e5, sample.recur(x, p))

round(prop.table(table(out)), 2)

me da una distribución aproximadamente uniforme

> round(prop.table(table(out)), 2)
out
   1    2    3    4    5    6    7 
0.14 0.14 0.15 0.14 0.14 0.14 0.14 

Es interesante tener en cuenta (si no he hecho un grave error) que esto produce un resultado diferente de la generación de $N$ binomial resultados como la suma de 13 lanzamientos de la moneda (recuento 0 cabezas como un resultado) y el mapeo de la {0,1,2,...,12,13} índice en la lista anterior de {1,2,3,...,3,2,1}. No sé bien cómo probar que mi método funciona.

---

La evaluación de la expectativa de valor para el número de monedas de lanza

La expectativa de valor para el número de monedas de los lanzamientos pueden ser calculados utilizando la matriz de transición de abajo (respondiendo a la pregunta cuando empezamos con $x$ no eliminado los números, a continuación, ¿cuál es la probabilidad de obtener a $y$ no eliminado los números)

$$M = \begin{bmatrix} q^7 & 0 & 0 & 0 &0 & 0 & 1 & 1\\ 7p^1q^6 & q^6 & 0 & 0 & 0 & 0 & 0 & 0 \\ 21p^2q^5 & 6p^1q^5 & q^5 & 0 & 0 & 0 & 0 & 0\\ 35 p^3q^4 & 15 p^2q^4 & 5q^4 & q^4 & 0 & 0 & 0 & 0\\ 35 p^4q^3 & 20 p^3q^3 & 10 p^2q^3 & 4 p^1q^3 & q^3 & 0 & 0 & 0\\ 21p^5q^2 & 15 p^4q^2 & 10 p^3q^2 & 6 p^2q^2 & 3 p^1q^2 & q^2 & 0 & 0\\ 7p^6q^1 & 6 p^5q^1 & 5 p^4q^1 & 4 p^3q^1 & 3p^2q^1 & 2p^1q^1 & 0 & 0\\ p^7 & p^6 & p^5 & p^4 & p^3 & p^2 & 0 & 0 \end{bmatrix}$$

The eigenvector associated with the eigenvalue 1 (which can be found by solving $(M-I)v=0$) depicts how much time is relatively spend in what state. Then 7th state is how often you will be able to draw a number from 1 to 7. The other states tell how many coin flips it costs.

Below is the image which compares with the answer from NcAdams which has expectation value for coin throws being $E(n) = \frac{24}{7}p(1-p)$

Remarkable is that the method performs better roughly for $p>2/3$. Pero también el rendimiento no es simétrica. Una simétrica y un mejor rendimiento general se podría hacer cuando un probabilística de la regla de conmutación sería que los cambios de la regla de decisión de las colas a los jefes cuando los jefes pasa a ser improbable.

La solución se encontró con wxMaxima

M: matrix(
 [(1-p)^7,        0,          0,0,0,0,1,1], 
 [7* p*(1-p)^6,   (1-p)^6,        0,0,0,0,0,0], 
 [21*p^2*(1-p)^5, 6*p*(1-p)^5,    (1-p)^5,0,0,0,0,0], 
 [35*p^3*(1-p)^4, 15*p^2*(1-p)^4, 5*p*(1-p)^4,(1-p)^4,0,0,0,0], 
 [35*p^4*(1-p)^3, 20*p^3*(1-p)^3, 10*p^2*(1-p)^3,4*p*(1-p)^3,(1-p)^3,0,0,0], 
 [21*p^5*(1-p)^2, 15*p^4*(1-p)^2, 10*p^3*(1-p)^2,6*p^2*(1-p)^2,3*p*(1-p)^2,(1-p)^2,0,0], 
 [7* p^6*(1-p)^1, 6*p^5*(1-p),    5*p^4*(1-p),4*p^3*(1-p),3*p^2*(1-p),2*(1-p)*p,0,0], 
 [p^7,        p^6,        p^5,p^4,p^3,p^2,0,0]
);
z: nullspace(M-diagmatrix(8,1));
x : apply (addcol, args (z));
t : [7,6,5,4,3,2,0,0];
plot2d(t.x/x[7],[p,0,1],logy);

Los cálculos en R

# plotting empty canvas
plot(-100,-100,
     xlab="p",
     ylab="E(total flips)",
     ylim=c(10,1000),xlim=c(0,1),log="y")

# plotting simulation
for (p in seq(0.1,0.9,0.05)) {

  n <- rep(0,10000)
  for (i in 1:10000) {
    success  = 0
    tests = c(1,1,1,1,1,1,1)     # start with seven numbers in the set
    count = 0
    while(success==0) {
      for (j in 1:7)  {
        if (tests[j]==1) {
          count = count + 1
          if  (rbinom(1,1,p) == 1) {
            tests[j] <- 0        # elliminate number when we draw heads
          }
        }
      }
      if (sum(tests)==1) {
        n[i] = count
        success = 1              # end     when 1 is left over
      }
      if (sum(tests)==0) {
        tests = c(1,1,1,1,1,1,1) # restart when 0 are left over
      }
    }
  }
  points(p,mean(n))
}

# plotting formula
p <- seq(0.001,0.999,0.001)

tot <- (7*(p^20-20*p^19+189*p^18-1121*p^17+4674*p^16-14536*p^15+34900*p^14-66014*p^13+99426*p^12-119573*p^11+114257*p^10-85514*p^9+48750*p^8-20100*p^7+5400*p^6-720*p^5)+6*
    (-7*p^21+140*p^20-1323*p^19+7847*p^18-32718*p^17+101752*p^16-244307*p^15+462196*p^14-696612*p^13+839468*p^12-806260*p^11+610617*p^10-357343*p^9+156100*p^8-47950*p^7+9240*p^6-840*p^5)+5*
    (21*p^22-420*p^21+3969*p^20-23541*p^19+98154*p^18-305277*p^17+733257*p^16-1389066*p^15+2100987*p^14-2552529*p^13+2493624*p^12-1952475*p^11+1215900*p^10-594216*p^9+222600*p^8-61068*p^7+11088*p^6-1008*p^5)+4*(-
    35*p^23+700*p^22-6615*p^21+39235*p^20-163625*p^19+509425*p^18-1227345*p^17+2341955*p^16-3595725*p^15+4493195*p^14-4609675*p^13+3907820*p^12-2745610*p^11+1592640*p^10-750855*p^9+278250*p^8-76335*p^7+13860*p^6-
    1260*p^5)+3*(35*p^24-700*p^23+6615*p^22-39270*p^21+164325*p^20-515935*p^19+1264725*p^18-2490320*p^17+4027555*p^16-5447470*p^15+6245645*p^14-6113275*p^13+5102720*p^12-3597370*p^11+2105880*p^10-999180*p^9+371000
   *p^8-101780*p^7+18480*p^6-1680*p^5)+2*(-21*p^25+420*p^24-3990*p^23+24024*p^22-103362*p^21+340221*p^20-896679*p^19+1954827*p^18-3604755*p^17+5695179*p^16-7742301*p^15+9038379*p^14-9009357*p^13+7608720*p^12-
 5390385*p^11+3158820*p^10-1498770*p^9+556500*p^8-152670*p^7+27720*p^6-2520*p^5))/(7*p^27-147*p^26+1505*p^25-10073*p^24+49777*p^23-193781*p^22+616532*p^21-1636082*p^20+3660762*p^19-6946380*p^18+11213888*p^17-
  15426950*p^16+18087244*p^15-18037012*p^14+15224160*p^13-10781610*p^12+6317640*p^11-2997540*p^10+1113000*p^9-305340*p^8+55440*p^7-5040*p^6)
lines(p,tot)

#plotting comparison with alternative method
lines(p,3*8/7/(p*(1-p)),lty=2)

legend(0.2,500,
       c("simulation","calculation","comparison"),
       lty=c(0,1,2),pch=c(1,NA,NA))