Как смоделировать доход с помощью испытаний Бернулли / вероятностей успеха в R?

Я хотел бы смоделировать сценарии доходов для: price и est_p (оценочная вероятность) из следующих df:

df <- data.frame(price        = c(200, 100, 600, 20, 100),
                 est_p        = c(0.9, 0.2, 0.8, 0.5, 0.6),
                 actual_sale  = c(FALSE, TRUE, TRUE, TRUE, TRUE))

Выручка - сумма price, где actual_sale равно TRUE:

print(actual1 <- sum(df$price[df$actual_sale])) # Actual Revenue

[1] 820

Я создал функцию для имитации испытаний Бернулли со значениями est_p и price:

bernoulli <- function(df) {
        sapply(seq(nrow(df)), function(x) {
                prc <- df$price[x]
                p   <- df$est_p[x]
                sample(c(prc, 0), size = 1000, replace = T, prob = c(p, 1 - p))
                })
}

И применил его к образцу df:

set.seed(100)
distr1 <- rowSums(bernoulli(df))
quantile(distr1)

  0%  25%  50%  75% 100% 
   0  700  820  920 1020 

Выглядит нормально, фактическое значение = медиана! Но когда я применяю ту же функцию к увеличенному (реплицированному в 1000 раз) образцу - df1000, фактический доход выходит за пределы смоделированных значений:

df1000 <- do.call("rbind", replicate(1000, df, simplify = FALSE))

print(actual2 <- sum(df1000$price[df1000$actual_sale])) 

[1] 820000

distr2 <- rowSums(bernoulli(df1000))
quantile(distr2)

    0%    25%    50%    75%   100% 
726780 744300 750050 754920 775800

Почему фактическая выручка выходит за рамки смоделированных значений? Где я сделал ошибку и как правильно решить эту проблему?