The following code is to calculate the expectation of a random variable of logit-normal distribution with parameters mu and sigma (mu is mu, and lsig is the logarithm of sigma).
fun5 = 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)
}
el = 17
integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
We should integrate this function from negative infinity to positive infinity, but I don't know how to do so and know only to integrate for finite interval. Thus, I try to integrate from 'wide enough' interval (from -el to +el). When the 'el' is greater than 0.5, it seems to work reasonably (true value of this integration is 0.585). But when el is 14 and 15, this works weirdly. Does anybody know why this happens?
> el = 10
> integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
[1] 0.585
> el = 13
> integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
[1] 0.585
> el = 14
> integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
[1] 2.975338e-05
> el = 15
> integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
[1] 1.134474e-05
> el = 16
> integrate(fun5,-el,el,mu=0.3434108,lsig=-3.5)$value
[1] 0.585