Eval очень большой факториал в R

Question

Мне нужно интегрировать функцию, которая имеет факториал в своем выражении. Но если вы попытаетесь оценить факториал, когда n> 170, R вернет Inf.

Я обнаружил множество пакетов, которые позволяют вычислять очень большие числа, однако они всегда возвращают объект из класса, который Я не могу интегрироваться. Конечным результатом от интеграла всегда будет небольшое число.

Вот мой код:

integrand <- function(n, i, x) {
    (factorial(n) / (factorial(i - 1) * factorial(n - i))) *
        x^(i - 1) * (1 - x)^(n - i)
}

forder <- function(Fx, x, n, i, ...) {
    lower <- sapply(x - 1, Fx, ...)
    upper <- sapply(x, Fx, ...)
     integrate(integrand,
               lower = lower, upper = upper, n = n, i = i,
               stop.on.error = FALSE)$value
}
forder <- Vectorize(forder, "x")

##------------------------------------------------------------------------------
## Some example
y <- sort(rpois(100, 1))

## Works fine
forder(ppois, y, 170, 10, lambda = 1)

## Does not work
forder(ppois, y, 171, 10, lambda = 1)
##------------------------------------------------------------------------------

Answer 1

Как сказано в моем комментарии, вы можете заменить (factorial(n) / (factorial(i - 1) * factorial(n - i))) на i*choose(n, i). Эти две величины равны, но choose(n,i) допускает более высокие значения n.

Или вместо численного интегрирования можно использовать функцию pbeta:

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  i*choose(n, i) * (pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)) * beta(i, n-i+1)
}

Еще лучше использовать логарифмы:

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  lg <- log(i) + lchoose(n, i) + 
    log(pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)) + lbeta(i, n-i+1)
  exp(lg)
}

---

РЕДАКТИРОВАТЬ

Я не заметил этого упрощения: i*choose(n, i) * beta(i, n-i+1) = 1. Так что вы можете просто сделать:

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)
}

Answer 2

Изменение integrand на использование логарифмов, оба вызова работают. Я также публикую функции @ Идеи Стефана Лорана об использовании choose/lchoose и pbeta/beta функций.

integrand <- function(n, i, x) {
  y <- lfactorial(n) - lfactorial(i - 1) - lfactorial(n - i) +
    (i - 1)*log(x) + log(1 - x)*(n - i)
  exp(y)
}

forder <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  integrate(integrand,
            lower = lower, upper = upper, n = n, i = i,
            stop.on.error = FALSE)$value
}
forder <- Vectorize(forder, "x")

integrandSL <- function(n, i, x) {
  y <- log(i) + lchoose(n, i) + (i - 1)*log(x) + log(1 - x)*(n - i)
  exp(y)
}

forderSL <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  integrate(integrandSL,
            lower = lower, upper = upper, n = n, i = i,
            stop.on.error = FALSE)$value
}
forderSL <- Vectorize(forderSL, "x")

forderSL2 <- function(Fx, x, n, i, ...) {
  lower <- sapply(x - 1, Fx, ...)
  upper <- sapply(x, Fx, ...)
  lg <- log(i) + lchoose(n, i) + 
    log(pbeta(upper, i, n-i+1) - pbeta(lower, i, n-i+1)) + lbeta(i, n-i+1)
  exp(lg)
}

Теперь тесты. Все результаты all.equal.

##-------------------------------------------------
set.seed(1234)
y <- sort(rpois(100, 1))

res_170 <- forder(ppois, y, 170, 10, lambda = 1)
res_171 <- forder(ppois, y, 171, 10, lambda = 1)

resSL_170 <- forderSL(ppois, y, 170, 10, lambda = 1)
resSL_171 <- forderSL(ppois, y, 171, 10, lambda = 1)

resSL2_170 <- forderSL2(ppois, y, 170, 10, lambda = 1)
resSL2_171 <- forderSL2(ppois, y, 171, 10, lambda = 1)

all.equal(res_170, resSL_170)    # TRUE
all.equal(res_170, resSL2_170)   # TRUE
all.equal(res_171, resSL_171)    # TRUE
all.equal(res_171, resSL2_171)   # TRUE