Integral function; hcubature runs forever

General
1

I used hcubature for double integrals over two parameters inside a function that I defined but it runs forever . How I can solve this problem to get the result of the integral over the function "myFunction" that I created. The following is the code.

library(TailRank)
library(numDeriv)
library(cubature)
set.seed(8765)
obs=15
try=10  
u=0.175
v=0.075
rbet=cbind(rbb(obs,try,u,v))  
rbet
combinations=combn(rbet, 2)
combinationst=t(combinations)
pair=combinationst[sample(nrow(combinationst),size=25,replace=FALSE),]
pair
myFunction=function(parameter)
{
  sumlog=sum(lchoose(try,rbet)+lbeta((parameter[1]*(1-parameter[2]))/parameter[2]+(rbet),(1-parameter[1])*(1-parameter[2])/parameter[2]+(try-rbet))-lbeta((parameter[1]*(1-parameter[2]))/parameter[2],(1-parameter[1])*(1-parameter[2])/parameter[2]))  
  z=0
  for (i in 1:nrow(pair))
  {
    a1=pair[i,1]
    logf1=lchoose(try,a1)+lbeta((parameter[1]*(1-parameter[2]))/parameter[2]+(a1),(1-parameter[1])*(1-parameter[2])/parameter[2]+(try-a1))-lbeta((parameter[1]*(1-parameter[2]))/parameter[2],(1-parameter[1])*(1-parameter[2])/parameter[2])   
    i1=seq(0,a1,1)
    a2=pair[i,2]
    logf2=lchoose(try,a2)+lbeta((parameter[1]*(1-parameter[2]))/parameter[2]+(a2),(1-parameter[1])*(1-parameter[2])/parameter[2]+(try-a2))-lbeta((parameter[1]*(1-parameter[2]))/parameter[2],(1-parameter[1])*(1-parameter[2])/parameter[2])  
    i2=seq(0,a2,1)
    myFun=function(x)
    { 
      f1=0
      f1=sum(choose(try,i1)*beta((x[1]*(1-x[2]))/x[2]+i1,(1-x[1])*(1-x[2])/x[2]+(try-i1))/beta((x[1]*(1-x[2]))/x[2],(1-x[1])*(1-x[2])/x[2]))
      f2=0
      f2=sum(choose(try,i2)*beta((x[1]*(1-x[2]))/x[2]+i2,(1-x[1])*(1-x[2])/x[2]+(try-i2))/beta((x[1]*(1-x[2]))/x[2],(1-x[1])*(1-x[2])/x[2])) 
      return(c(f1,f2))
    }
    Rjac=jacobian(func=myFun,x=c(parameter[1],parameter[2]),method="Richardson")
    labsdetjac=log(abs(det(Rjac)))
    subtotal=sumlog+labsdetjac-(logf1+logf2)
    total=z+subtotal
  }
  exptotal=exp(total)
  return(exptotal)
}
out=hcubature(myFunction, c(0.001,0.001), c(0.7,0.8),tol=1e-4)
result=out$integral
print(result)
2

A brief glance at the code for the integrand suggests to me that it might be hard to vectorise, which means it will be quite slow in R. It's probably a good case for using Rcpp, and calling cubature on the compiled code that is called many times.

But first, I would make sure that there isn't an issue with the code itself – cubature has a maxEval argument that you could set to check that it's just taking a long time, and nothing else is preventing it from arriving at a result.

1 Like
3

Thank you for your suggestion. I used tol=1e-4 and for some values of lower and upper limits of the integrals, I could get the results of the integrals directly but for another limits such as the one in the program (c(0.001,0.001), c(0.7,0.8)), it runs forever. I worked with your suggestion and used maxEval argument in this case and it works.If maxEval=10, the result of the integral= 2.096755e-21 but when maxEval=100, the result= 3.347127e-21.

So, how can I generally choose the appropriate number for maxEval that works for different upper and lower limits of the multidimensional integrals to avoid the problem that integral function runs forever?
Is it appropriate to choose tol to be maximum (i.e. tol=1e-5) or it will not affect the integral calculations?
Thanks