-1

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?

user67275
  • 1
  • 9
  • 38
  • 64

2 Answers2

3

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
adaien
  • 1,932
  • 1
  • 12
  • 26
1

We can use, that 
exp(y)/(1+exp(y)) is the same as (1 - 1/(1+exp(y))) or also1/(1+exp(-y)) 

fun5 <- function(y,mu=mu0,lsig=lsig0) 1/(1+exp(-y)) / sqrt(2*pi)/exp(lsig) * exp(-(y-mu)^2/2/exp(lsig)^2)
integrate(fun5,-Inf,Inf,mu=2.198216,lsig=-3)$value

.

> integrate(fun5,-Inf,Inf,mu=2.198216,lsig=-3)$value
[1] 0.9
jogo
  • 12,469
  • 11
  • 37
  • 42