1

My goal is to find solution to the following integration in R. Integral to solve

I have tried the following code, but it does not work.

my_func <- function(x){
  n=45
  alpha_const=1
  beta_const=1
  term_11 <- c()
  
  for (i in 0:n) {
    term_11[i+1] <- (n-i) / (x+n-i)
  }
  
  return( (x^(alpha_const-1)) * (exp(-beta_const*x)) * prod(term_11) )
}

integrate(my_func,0,Inf)

The solution I get is

0 with absolute error < 0
There were 46 warnings (use warnings() to see them)

When I looked in to the warnings,

Warning messages:
1: In term_11[i + 1] <- (n - i)/(x + n - i) :
  number of items to replace is not a multiple of replacement length
2: In term_11[i + 1] <- (n - i)/(x + n - i) :
  number of items to replace is not a multiple of replacement length
3: In term_11[i + 1] <- (n - i)/(x + n - i) :
  number of items to replace is not a multiple of replacement length
4: In term_11[i + 1] <- (n - i)/(x + n - i) :
  number of items to replace is not a multiple of replacement length
5: In term_11[i + 1] <- (n - i)/(x + n - i) :
  number of items to replace is not a multiple of replacement length
6: In term_11[i + 1] <- (n - i)/(x + n - i) :

which went on till "46: In term...".
I think I am not handling the product term in the integrand properly. Is it possible to solve this integral numerically? How do I work with the product term here?

1 Answers1

1

First, check the definition of your function to integrate as when i=n the factor in the product is zero, so the function is identically zero.

Leaving out that factor and rewriting your code in better R style, I get:

my_func <- function(x){
  n=45
  range <- 0:(n-1)
  alpha_const=1
  beta_const=1
  
  prod_term <- exp(sum(log(n-range)-log(x+n-range)))
  
  return( (x^(alpha_const-1)) * (exp(-beta_const*x)) * prod_term )
}

 integrate(my_func, 0, Inf)
1.688208e-07 with absolute error < 2.1e-08
kjetil b halvorsen
  • 1,206
  • 2
  • 18
  • 28
  • 2
    The main problem was that the function was not vectorized. You fixed that. – Stéphane Laurent Dec 11 '22 at 20:21
  • I understood the issue of "i=n" case. When I compare the result above with the result I get after modifying the code as per @akrun I see that it is different, 1.154001e-114 with absolute error < 1.4e-115. Any idea why that would be the case? – Arnoneel Sinha Dec 11 '22 at 20:36