0

I want to calculate numerical integral for the following function

library(cubature)
library(dplyr)
fun4=function(y1,y2,b0,b1,x)
{
  y=c(y1,y2)
  num=   lapply(1:length(b0), function(u) exp(b0[[u]] + t(b1[[u]])%*%y ) )
  num_l=lapply(1:length(x), function(u)  nth(num[[u]], x[u])  )
  den =  lapply(1:length(b0), function(u) 1+ sum(num[[u]]) )
  pp =  lapply(1:length(b0), function(u)  num_l[[u]] / den[[u]]  )
  gpp=unlist(pp)
  res=prod(gpp)
}

I used hcubature but I got NaN

  x2=c(1,3)
  b01=c(1,3,-3)
  b02=c(2,-2,0)
  b11=matrix(c(1,2,3,-2,1,0),ncol=3,nrow=2)
  b12=matrix(c(2,-2,3,1,0,4),ncol=3,nrow=2)
  beta0=list(b01,b02)
  beta1=list(b11,b12)
  fintg = function(X){
    y1=X[1]
    y2= X[2]
    z= X[3]
   log(z)*fun4(y1,y2,beta0,beta1,x2) 
  }
  loLimt = c(-Inf, -Inf,0)
  hiLimt = c( Inf, Inf,Inf)
  hcubature(fintg, loLimt, hiLimt)$integral

I'm not sure why I got NaN, but I looked at this question which has similar problem NaN and one of the answers mention that the function exp(y)/1+exp(y) is rounded to NaN if y is too big. A part of my function has something similar to that function, exp(b0 + t(b1)%*%y )/ 1+sum(exp(b0 + t(b1)%*%y )) . I do not know how to add ifelse condition in my function or if there is another way to solve NaN problem.

Adam
  • 131
  • 6
  • 1
    If the only possibility for `NaN` values occurring is that they result from `Inf/(1+Inf)`, you could just fill in the limit, i.e., `gpp <- ifelse(is.nan(gpp), 1, gpp)`. – Roland Feb 26 '21 at 06:54

1 Answers1

0

Your num is throwing an Inf, which causes the NaN. Would replacing the Inf with 1 work for you?

fun4=function(y1,y2,b0,b1,x)
{
  y=c(y1,y2)
  num=  lapply(1:length(b0), function(u) exp(b0[[u]] + t(b1[[u]])%*%y ) )

  num2=unlist(num)
  if(Inf %in% abs(num2)) num2[num2==Inf | num2==-Inf]=1
  num=list(as.matrix(num2[1:3]),as.matrix(num2[4:6]))
  
  num_l=lapply(1:length(x), function(u)  nth(num[[u]], x[u])  )
  den =  lapply(1:length(b0), function(u) 1+ sum(num[[u]]) )
  pp =  lapply(1:length(b0), function(u)  num_l[[u]] / den[[u]]  )
  gpp=unlist(pp)
  res=prod(gpp)
  res
}
rdodhia
  • 350
  • 2
  • 9