математика с большим n (> 100)

Question

Обещаю, это не , просто еще одна проблема, связанная с выполнением домашних заданий. Я реализовал функцию для расчета вероятности получения меньше суммы s при броске n m односторонних костей. Моя функция работает для малых значений n, но я нахожу странные результаты для больших значений n. Смотрите прикрепленный сюжет. Кто-нибудь может понять, что происходит?

Моя функция вероятности

Реализована из этого Математический стек обмена

probability <- function(s, m, n) {

  i <- 0:((s-1-n) / m)
  m^(-n) * sum((-1)^i * choose(n, i) * choose(s - 1 - i * m, n))

}

Начинает ломаться с ~ n > 80

n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
p <- mapply(probability, s = s, m = m, n = n)
plot(n, p, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"))

Answer 1

Как уже упоминалось в комментариях к исходному вопросу, проблема в том, что функция вероятности запрашивает R для вычисления действительно огромных чисел (choose(80,40) = 1.075072e+23), и мы достигаем числовых пределов точности для R.

Альтернативный подход, который не включает в себя огромные числа, но вместо этого использует много чисел, состоит в запуске симуляций Монте-Карло. Это генерирует распределение сумм бросков костей и сравнивает наблюдаемую сумму с распределением. Это займет больше времени, но будет намного проще и не будет иметь проблем с числовой точностью.

mc <- Vectorize(function(s, m, n, reps = 10000) {
  x <- replicate(reps, sum(sample(m, n, replace = TRUE)))
  ecdf(x)(s-1)
})



n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces
analytic_prob <- mapply(probability, s = s, m = m, n = n)
mc_prob <- mapply(mc, s = s, m = m, n = n)


plot(n, analytic_prob, main = paste("probability of rolling less than floor(", mean(1:m),"* n) as sum of n rolls of a", m, "sided dice"),
     sub = "monte carlo in red")
points(n, mc_prob, col = "red")

Answer 2

Проблема вызвана пределами точности чисел R. Как отметили комментаторы, n выбирают значения k, которые я вычисляю выше, действительно, очень большие (choose(80,40) = 1.075072e+23).

Мы можем использовать журналы, чтобы попытаться удержать проблему в вычислительных пределах R. Это реализация подхода Рамануджана. К сожалению, ошибки в приближении усугубляются, а точность уменьшается еще быстрее. Функция вероятности требует сложения и вычитания последовательности очень больших чисел, чтобы получить окончательное значение между 0 и 1, и не допускает никаких неточностей.

0) Перепишите функцию вероятности, чтобы разбить ее на шаги

probability <- function(s, m, n) {

  # Probability of getting less than s
  i <- 0:((s-1-n) / m)

  c1 <- choose(n, i)
  c2 <- choose(s - 1 - i * m, n)

  seq <- (-1)^i * (c1 * c2)

  m^(-n) * sum(seq)

}

1) реализовать аппроксимацию log (x!)

# using the 'ramanujan' method
ramanujan <- function(n){
  n * log(n) - n + log(n * (1 + 4*n * (1 + 2*n))) / 6 + log(pi) / 2
}

# confirm Ramanujan works correctly
n <- 1:200
diff <- log(factorial(n)) - ramanujan(n)
plot(n, diff) # r returns inf for factorial(171), but up to there the numbers match

2) переписать функцию choose с использованием аппроксимации журнала.

#' This function returns log(choose(n,k)) 
log_nck <- Vectorize(function(n, k) {
  if(n <= k | n < 1 | k < 1) return(log(choose(n,k))) # logs don't like 0 or neg numbers

  return((ramanujan(n) - ramanujan(k) - ramanujan(n-k)))
})

# Check that choose function works
n <- seq(10, 100, 10)
k <- seq(5, 50, 5)
c_real <- log(choose(n, k))
c_approx <- log_nck(n, k)
# If we print them, they appear to match
print(c_real)
print(c_approx)
# and the difference shows pretty small errors. 
print(c_real - c_approx)

3) переписать функция вероятности с использованием журнала выбора.

new_probability <- function(s, m, n) {

  # Probability of getting less than s
  i <- 0:((s-1-n) / m)

  c1 <- log_nck(n, i)
  c2 <- log_nck(s - 1 - i * m, n)

  seq <- (-1)^i * exp(c1 + c2)

  return(m^(-n) * sum(seq))

}

Окончательное тестирование

n <- 1:90 # number of dice
m <- 6 # number of sides
s <- floor(mean(1:m)*n) # sum of faces

p <- mapply(probability, s = s, m = m, n = n)
newp <- mapply(new_probability, s = s, m = m, n = n)

plot(n, p, main = "Original in black, approximation in red")
points(n, newp, col = "red")