Reputation: 112
I need to integrate a function which has a factorial in its expression. But, if you try to evaluate factorial, when n > 170, R returns Inf.
I found a lot of packages that allow you to calculate very large numbers, however, they always returns an object from a class that I can't integrate. The final result from the integral always will be a small number.
Here's my code:
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)
##------------------------------------------------------------------------------
Upvotes: 3
Views: 294
Reputation: 76683
Changing the integrand
to use logarithms both calls work. I also post functions @StéphaneLaurent's ideas of using choose/lchoose
and the pbeta/beta
functions.
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)
}
Now the tests. All results are 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
Upvotes: 1
Reputation: 84719
As said in my comment, you can replace (factorial(n) / (factorial(i - 1) * factorial(n - i)))
with i*choose(n, i)
. These two quantities are equal but choose(n,i)
allows higher values of n
.
Or you can use the pbeta
function instead of doing a numerical integration:
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)
}
Even better, use logarithms:
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 didn't notice this simplification: i*choose(n, i) * beta(i, n-i+1) = 1
. So you can simply do:
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)
}
Upvotes: 6