Search code examples
rintegralnumerical-integration

How to double integrate with +inf?


So i have this code in R that tries to double integrate a function f(x,y) using integral2, where one domain is (0.5 , +inf), but integral2 doesn't support it (because of that infinite). My question is: do you know a way to make it work?

 fprob8 <- function(a , b)
{
  f <- function(x , y)
  {
    densGamma <- function(x)
    {
      numitor <- b ^ a * fgam(a)
      numarator <- x ^ (a - 1) * exp(-x / b)
      return (numarator / numitor)
    }

    densBeta <- function(x)
    {
      numitor <- fbet(a , b)
      numarator <- (1 - x) ^ (b - 1) * x ^ (a - 1)
      return (numarator/numitor)
    }

    return (densGamma(x) * densBeta(y)) 
  }

  ymin <- function(x) # limita superioara pentru cea de-a doua integrala
  {
    return (min(x - 0.5),1)
  }

  I <- integral2(f, 0.5 , Inf , ymin , 1)
  return(I)
}

Solution

  • Use the cubature package. It evaluates multiple integrals, allowing infinite bounds.

    Here is an example:

    library(cubature)
    f <- function(x) exp(-x[1]-x[2])
    pcubature(f, c(0,0), c(Inf,Inf))$integral
    # 1
    
    library(pracma)
    f2 <- function(x,y) f(c(x,y))
    integral2(f2, 0, 1000, 0, 1000, vectorized = FALSE)$Q
    # 1