-
- exactComp <- FALSE #TODO: global, or argument...
-
- if (exactComp && link == "probit")
- {
- # Use exact computations
- sapply( seq_along(λ), function(k) {
- .exactProbitIntegral(order, λ[k], b[k])
- })
- }
-
- else
- {
- # Numerical integration
- sapply( seq_along(λ), function(k) {
- res <- NULL
- tryCatch({
- # Fast code, may fail:
- res <- stats::integrate(
- function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
- lower=-Inf, upper=Inf )$value
- }, error = function(e) {
- # Robust slow code, no fails observed:
- sink("/dev/null") #pracma package has some useless printed outputs...
- res <- pracma::integral(
- function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
- xmin=-Inf, xmax=Inf, method="Kronrod")
- sink()
- })
- res
- })
- }
-}
-
-# TODO: check these computations (wrong atm)
-.exactProbitIntegral <- function(order, λ, b)
-{
- c1 = (1/sqrt(2*pi)) * exp( -.5 * b/((λ^2+1)^2) )
- if (order == 1)
- return (c1)
- c2 = b - λ^2 / (λ^2+1)
- if (order == 2)
- return (c1 * c2)
- if (order == 3)
- return (c1 * (λ^2 - 1 + c2^2))
- if (order == 4)
- return ( (c1*c2/((λ^2+1)^2)) * (-λ^4*((b+1)^2+1) -
- 2*λ^3 + λ^2*(2-2*b*(b-1)) + 6*λ + 3 - b^2) )
- if (order == 5) #only remaining case...
- return ( c1 * (3*λ^4+c2^4+6*c1^2*(λ^2-1) - 6*λ^2 + 6) )
+ sapply( seq_along(λ), function(k) {
+ res <- NULL
+ tryCatch({
+ # Fast code, may fail:
+ res <- stats::integrate(
+ function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
+ lower=-Inf, upper=Inf )$value
+ }, error = function(e) {
+ # Robust slow code, no fails observed:
+ sink("/dev/null") #pracma package has some useless printed outputs...
+ res <- pracma::integral(
+ function(z) .deriv[[link]][[order]](λ[k]*z+b[k]) * exp(-z^2/2) / sqrt(2*pi),
+ xmin=-Inf, xmax=Inf, method="Kronrod")
+ sink()
+ })
+ res
+ })