[R] Integrate: compound distribution

[tsakk]

Dear All,

thanks in advance for your attention. I'm struggling with this for a few
weeks now, and figured it's time to get some help. I want to

integrate(f(x), lower=-Inf, upper=Inf)

with f(x) =
((gamma(K+1)/(gamma(r+1)*gamma(K-r+1)))*(q(x)^r)*(((1-q(x))^(K-r))*phi(x),

where phi(x) is the standard normal pdf, and q(x) the logistic CDF (with
inverse scale parameter, so a little reparameterized). In short: it's a
compound binomial probability. It can be shown that this integral can be
expressed as a sum of moments of the System-Bounded Johnson distribution,
for which it is proven that no analytical expressions exist. Moreover, it
can be written as an infinite Taylor series in terms of moments of the
lognormal distribution, but this series doesn't converge in general.
Numerical integration is the only way to go.

The problem is that integrate() is not accurate enough. I cannot compensate
using rel.tol, since there is always a trade-off between "integral is
probably divergent" and "roundoff error detected". I've read about the Rmpfr
package, but this doesn't seem like an elegant way to proceed.

I'm looking forward to your input. Thanks in advance! Cheers,

Thomas





--
View this message in context: http://r.789695.n4.nabble.com/Integrate-compound-distribution-tp4637735.html
Sent from the R help mailing list archive at Nabble.com.