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
I'm wondering if you understand that R designers allow -Inf and Inf as bounds and in fact encourages user to employ them, especially when one or both of those bounds is far from what might be called the "predominant support":
> integrate(fun5,-Inf,Inf, mu=0.3434108, lsig=-3.5)$value
[1] 0.585