# Cartas numeradas de 1 a 10.
# Retira-se uma carta por vez ate sair uma par (sucesso).
# Conta-se o numero de cartas retiradas ate o sucesso.

get_n <- function() {
  cartas <- 1:10
  n <- 0
  carta <- 11
  while(carta %% 2 != 0) {
    carta <- sample(cartas, 1)
    n = n + 1
    if(carta %% 2 == 0) {
      p <- 1/length(cartas)
    } else {
      cartas <- cartas[cartas != carta]
    }
  }
  return(cbind(n, p))
}

N <- 100000
nn <- NULL
for (i in 1:N) nn <- rbind(nn, get_n())

(x <- prop.table(table(nn[,1])))
## 
##       1       2       3       4       5       6 
## 0.50085 0.27759 0.13813 0.05966 0.01952 0.00425
barplot(x, space = 0, xlab = 'Nº de retiradas até sair carta par.')
y <- fitdistrplus::fitdist(nn[,1] - 1, distr = 'geom')

summary(y)
## Fitting of the distribution ' geom ' by maximum likelihood 
## Parameters : 
##       estimate  Std. Error
## prob 0.5458039 0.001163206
## Loglikelihood:  -126225.8   AIC:  252453.6   BIC:  252463.1
plot(y)

x <- prop.table(table(nn[,2]))
names(x) <- round(as.numeric(names(x)), 4)
x
##     0.1  0.1111   0.125  0.1429  0.1667     0.2 
## 0.50085 0.27759 0.13813 0.05966 0.01952 0.00425
barplot(x, space = 0, xlab = 'Histograma de p')