rnannumerical-integration

'non-finite' function value error in 'integrate' function in R

I defined a function called 'fun5' as follows:

function(y,mu=mu0,lsig=lsig0) {
  res = exp(y)/(1+exp(y)) * 1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2)
  return(res)

, then integrated the function from negative infinity to positive infinity with two parameters. 

integrate(fun5,-Inf,Inf,mu=2.198216,lsig=-3)$value

This integral gives the expectation of a random variable which has logit-normal distribution with mu = 2.198216 and sigma = exp(-3).

This error occurred.

Error in integrate(fun5, -Inf, Inf, mu = 2.198216, lsig = -3) : 
  non-finite function value

Since the function 'fun5' is a random variable between 0 and 1 multiplied by probability density, it should be positive everywhere, though it might be very close to zero. I don't understand why it has non-finite value somewhere.

Could anybody give an advice?

Solution

  • The problem is that the function 

    exp(y)/(1+exp(y))

    is rounded to NaN when y is too big. You can avoid this replacing it with 1 when y is too big. This function will play the trick:

    fun5<-function(y,mu=mu0,lsig=-lsig0) {
res = ifelse(y<100, exp(y)/(1+exp(y)) * 1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2),
             1/sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2))
return(res)}

    and now this will work

    integrate(fun5,-Inf,Inf,mu=2.198216,lsig=-3)$value
[1] 0.9