X-Git-Url: https://git.auder.net/?p=morpheus.git;a=blobdiff_plain;f=pkg%2FR%2FoptimParams.R;h=039070c5280b847af077056b2cfcc89b6c73b589;hp=4f886ac3fbd60139c8350e2222cb382c791b49d0;hb=2591fa8343c69ddb94dec5e55871d34c55eff9a3;hpb=98b8a5ddffdce7e0b63746d4b58bb923049dca7d

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index 4f886ac..039070c 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -1,17 +1,9 @@
-#' Optimize parameters
-#'
-#' Optimize the parameters of a mixture of logistic regressions model, possibly using
-#' \code{mu <- computeMu(...)} as a partial starting point.
+#' Wrapper function for OptimParams class
 #'
 #' @param K Number of populations.
 #' @param link The link type, 'logit' or 'probit'.
-#' @param optargs a list with optional arguments:
-#'   \itemize{
-#'     \item 'M' : list of moments of order 1,2,3: will be computed if not provided.
-#'     \item 'X,Y' : input/output, mandatory if moments not given
-#'     \item 'exact': use exact formulas when available?
-#'     \item weights Weights on moments when minimizing sum of squares
-#'   }
+#' @param X Data matrix of covariables
+#' @param Y Output as a binary vector
 #'
 #' @return An object 'op' of class OptimParams, initialized so that \code{op$run(x0)}
 #'   outputs the list of optimized parameters
@@ -20,8 +12,8 @@
 #'     \item β: regression matrix, size dxK
 #'     \item b: intercepts, size K
 #'   }
-#'   x0 is a vector containing respectively the K-1 first elements of p, then β by
-#'   columns, and finally b: \code{x0 = c(p[1:(K-1)],as.double(β),b)}.
+#'   θ0 is a vector containing respectively the K-1 first elements of p, then β by
+#'   columns, and finally b: \code{θ0 = c(p[1:(K-1)],as.double(β),b)}.
 #'
 #' @seealso \code{multiRun} to estimate statistics based on β, and
 #'   \code{generateSampleIO} for I/O random generation.
@@ -30,228 +22,276 @@
 #' # Optimize parameters from estimated μ
 #' io = generateSampleIO(10000, 1/2, matrix(c(1,-2,3,1),ncol=2), c(0,0), "logit")
 #' μ = computeMu(io$X, io$Y, list(K=2))
-#' M <- computeMoments(io$X, io$Y)
-#' o <- optimParams(2, "logit", list(M=M))
-#' x0 <- c(1/2, as.double(μ), c(0,0))
-#' par0 <- o$run(x0)
+#' o <- optimParams(io$X, io$Y, 2, "logit")
+#' θ0 <- list(p=1/2, β=μ, b=c(0,0))
+#' par0 <- o$run(θ0)
 #' # Compare with another starting point
-#' x1 <- c(1/2, 2*as.double(μ), c(0,0))
-#' par1 <- o$run(x1)
+#' θ1 <- list(p=1/2, β=2*μ, b=c(0,0))
+#' par1 <- o$run(θ1)
 #' o$f( o$linArgs(par0) )
 #' o$f( o$linArgs(par1) )
 #' @export
-optimParams = function(K, link=c("logit","probit"), optargs=list())
+optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL)
 {
-	# Check arguments
-	link <- match.arg(link)
-	if (!is.list(optargs))
-		stop("optargs: list")
-	if (!is.numeric(K) || K < 2)
-		stop("K: integer >= 2")
-
-	M <- optargs$M
-	if (is.null(M))
-	{
-		if (is.null(optargs$X) || is.null(optargs$Y))
-			stop("If moments are not provided, X and Y are required")
-		M <- computeMoments(optargs$X,optargs$Y)
-	}
-
-  weights <- optargs$weights
-  if (is.null(weights))
-    weights <- rep(1, K)
-
-	# Build and return optimization algorithm object
-	methods::new("OptimParams", "li"=link, "M1"=as.double(M[[1]]),
-		"M2"=as.double(M[[2]]), "M3"=as.double(M[[3]]),
-    "weights"=weights, "K"=as.integer(K))
+  # Check arguments
+  if (!is.matrix(X) || any(is.na(X)))
+    stop("X: numeric matrix, no NAs")
+  if (!is.numeric(Y) || any(is.na(Y)) || any(Y!=0 & Y!=1))
+    stop("Y: binary vector with 0 and 1 only")
+  link <- match.arg(link)
+  if (!is.numeric(K) || K!=floor(K) || K < 2)
+    stop("K: integer >= 2")
+
+  if (is.null(M))
+  {
+    # Precompute empirical moments
+    Mtmp <- computeMoments(X, Y)
+    M1 <- as.double(Mtmp[[1]])
+    M2 <- as.double(Mtmp[[2]])
+    M3 <- as.double(Mtmp[[3]])
+    M <- c(M1, M2, M3)
+  }
+  else
+    M <- c(M[[1]], M[[2]], M[[3]])
+
+  # Build and return optimization algorithm object
+  methods::new("OptimParams", "li"=link, "X"=X,
+    "Y"=as.integer(Y), "K"=as.integer(K), "Mhat"=as.double(M))
 }
 
-# Encapsulated optimization for p (proportions), β and b (regression parameters)
-#
-# @field li Link, 'logit' or 'probit'
-# @field M1 Estimated first-order moment
-# @field M2 Estimated second-order moment (flattened)
-# @field M3 Estimated third-order moment (flattened)
-# @field weights Vector of moments' weights
-# @field K Number of populations
-# @field d Number of dimensions
-#
+#' Encapsulated optimization for p (proportions), β and b (regression parameters)
+#'
+#' Optimize the parameters of a mixture of logistic regressions model, possibly using
+#' \code{mu <- computeMu(...)} as a partial starting point.
+#'
+#' @field li Link function, 'logit' or 'probit'
+#' @field X Data matrix of covariables
+#' @field Y Output as a binary vector
+#' @field K Number of populations
+#' @field d Number of dimensions
+#' @field W Weights matrix (iteratively refined)
+#'
 setRefClass(
-	Class = "OptimParams",
-
-	fields = list(
-		# Inputs
-		li = "character", #link 'logit' or 'probit'
-		M1 = "numeric", #order-1 moment (vector size d)
-		M2 = "numeric", #M2 easier to process as a vector
-		M3 = "numeric", #M3 easier to process as a vector
-		# Dimensions
-		K = "integer",
-		d = "integer"
-	),
-
-	methods = list(
-		initialize = function(...)
-		{
-			"Check args and initialize K, d"
-
-			callSuper(...)
-			if (!hasArg("li") || !hasArg("M1") || !hasArg("M2") || !hasArg("M3")
-				|| !hasArg("K"))
-			{
-				stop("Missing arguments")
-			}
-
-			d <<- length(M1)
-		},
-
-		expArgs = function(x)
-		{
-			"Expand individual arguments from vector x"
-
-			list(
-				# p: dimension K-1, need to be completed
-				"p" = c(x[1:(K-1)], 1-sum(x[1:(K-1)])),
-				"β" = matrix(x[K:(K+d*K-1)], ncol=K),
-				"b" = x[(K+d*K):(K+(d+1)*K-1)])
-		},
-
-		linArgs = function(o)
-		{
-			" Linearize vectors+matrices into a vector x"
-
-			c(o$p[1:(K-1)], as.double(o$β), o$b)
-		},
-
-		f = function(x)
-		{
-			"Sum of squares (Mi - hat_Mi)^2 where Mi is obtained from formula"
-
-			P <- expArgs(x)
-			p <- P$p
-			β <- P$β
-			λ <- sqrt(colSums(β^2))
-			b <- P$b
-
-			# Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1
-			β2 <- apply(β, 2, function(col) col %o% col)
-			β3 <- apply(β, 2, function(col) col %o% col %o% col)
-
-			return(
-				weights[1] * sum( ( β  %*% (p * .G(li,1,λ,b)) - M1 )^2 ) +
-				weights[2] * sum( ( β2 %*% (p * .G(li,2,λ,b)) - M2 )^2 ) +
-				weights[3] * sum( ( β3 %*% (p * .G(li,3,λ,b)) - M3 )^2 ) )
-		},
-
-		grad_f = function(x)
-		{
-			"Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
-
-			P <- expArgs(x)
-			p <- P$p
-			β <- P$β
-			λ <- sqrt(colSums(β^2))
-			μ <- sweep(β, 2, λ, '/')
-			b <- P$b
-
-			# Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1
-			β2 <- apply(β, 2, function(col) col %o% col)
-			β3 <- apply(β, 2, function(col) col %o% col %o% col)
-
-			# Some precomputations
-			G1 = .G(li,1,λ,b)
-			G2 = .G(li,2,λ,b)
-			G3 = .G(li,3,λ,b)
-			G4 = .G(li,4,λ,b)
-			G5 = .G(li,5,λ,b)
-
-			# (Mi - hat_Mi)^2 ' == (Mi - hat_Mi)' 2(Mi - hat_Mi) = Mi' Fi
-			F1 = as.double( 2 * ( β  %*% (p * G1) - M1 ) )
-			F2 = as.double( 2 * ( β2 %*% (p * G2) - M2 ) )
-			F3 = as.double( 2 * ( β3 %*% (p * G3) - M3 ) )
-
-			km1 = 1:(K-1)
-			grad <- #gradient on p
-				weights[1] * t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 +
-				weights[2] * t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ) %*% F2 +
-				weights[3] * t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ) %*% F3
-
-			grad_β <- matrix(nrow=d, ncol=K)
-			for (i in 1:d)
-			{
-				# i determines the derivated matrix dβ[2,3]
-
-				dβ_left <- sweep(β, 2, p * G3 * β[i,], '*')
-				dβ_right <- matrix(0, nrow=d, ncol=K)
-				block <- i
-				dβ_right[block,] <- dβ_right[block,] + 1
-				dβ <- dβ_left + sweep(dβ_right, 2,  p * G1, '*')
-
-				dβ2_left <- sweep(β2, 2, p * G4 * β[i,], '*')
-				dβ2_right <- do.call( rbind, lapply(1:d, function(j) {
-					sweep(dβ_right, 2, β[j,], '*')
-				}) )
-				block <- ((i-1)*d+1):(i*d)
-				dβ2_right[block,] <- dβ2_right[block,] + β
-				dβ2 <- dβ2_left + sweep(dβ2_right, 2, p * G2, '*')
-
-				dβ3_left <- sweep(β3, 2, p * G5 * β[i,], '*')
-				dβ3_right <- do.call( rbind, lapply(1:d, function(j) {
-					sweep(dβ2_right, 2, β[j,], '*')
-				}) )
-				block <- ((i-1)*d*d+1):(i*d*d)
-				dβ3_right[block,] <- dβ3_right[block,] + β2
-				dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*')
-
-				grad_β[i,] <-
-          weights[1] * t(dβ) %*% F1 +
-          weights[2] * t(dβ2) %*% F2 +
-          weights[3] * t(dβ3) %*% F3
-			}
-			grad <- c(grad, as.double(grad_β))
-
-			grad = c(grad, #gradient on b
-				weights[1] * t( sweep(β,  2, p * G2, '*') ) %*% F1 +
-				weights[2] * t( sweep(β2, 2, p * G3, '*') ) %*% F2 +
-				weights[3] * t( sweep(β3, 2, p * G4, '*') ) %*% F3 )
-
-			grad
-		},
-
-		run = function(x0)
-		{
-			"Run optimization from x0 with solver..."
-
-	    if (!is.list(x0))
-		    stop("x0: list")
-      if (is.null(x0$β))
-        stop("At least x0$β must be provided")
-			if (!is.matrix(x0$β) || any(is.na(x0$β)) || ncol(x0$β) != K)
-				stop("x0$β: matrix, no NA, ncol == K")
-      if (is.null(x0$p))
-        x0$p = rep(1/K, K-1)
-      else if (length(x0$p) != K-1 || sum(x0$p) > 1)
-        stop("x0$p should contain positive integers and sum to < 1")
-      # Next test = heuristic to detect missing b (when matrix is called "beta")
-      if (is.null(x0$b) || all(x0$b == x0$β))
-        x0$b = rep(0, K)
-      else if (any(is.na(x0$b)))
-        stop("x0$b cannot have missing values")
-
-			op_res = constrOptim( linArgs(x0), .self$f, .self$grad_f,
-				ui=cbind(
-					rbind( rep(-1,K-1), diag(K-1) ),
-					matrix(0, nrow=K, ncol=(d+1)*K) ),
-				ci=c(-1,rep(0,K-1)) )
-
-			expArgs(op_res$par)
-		}
-	)
+  Class = "OptimParams",
+
+  fields = list(
+    # Inputs
+    li = "character", #link function
+    X = "matrix",
+    Y = "numeric",
+    Mhat = "numeric", #vector of empirical moments
+    # Dimensions
+    K = "integer",
+    n = "integer",
+    d = "integer",
+    # Weights matrix (generalized least square)
+    W = "matrix"
+  ),
+
+  methods = list(
+    initialize = function(...)
+    {
+      "Check args and initialize K, d, W"
+
+      callSuper(...)
+      if (!hasArg("X") || !hasArg("Y") || !hasArg("K")
+        || !hasArg("li") || !hasArg("Mhat"))
+      {
+        stop("Missing arguments")
+      }
+
+      n <<- nrow(X)
+      d <<- ncol(X)
+      # W will be initialized when calling run()
+    },
+
+    expArgs = function(v)
+    {
+      "Expand individual arguments from vector v into a list"
+
+      list(
+        # p: dimension K-1, need to be completed
+        "p" = c(v[1:(K-1)], 1-sum(v[1:(K-1)])),
+        "β" = t(matrix(v[K:(K+d*K-1)], ncol=d)),
+        "b" = v[(K+d*K):(K+(d+1)*K-1)])
+    },
+
+    linArgs = function(L)
+    {
+      "Linearize vectors+matrices from list L into a vector"
+
+      # β linearized row by row, to match derivatives order
+      c(L$p[1:(K-1)], as.double(t(L$β)), L$b)
+    },
+
+    computeW = function(θ)
+    {
+      require(MASS)
+      dd <- d + d^2 + d^3
+      M <- Moments(θ)
+      Omega <- matrix( .C("Compute_Omega",
+        X=as.double(X), Y=as.integer(Y), M=as.double(M),
+        pn=as.integer(n), pd=as.integer(d),
+        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
+      MASS::ginv(Omega)
+    },
+
+    Moments = function(θ)
+    {
+      "Vector of moments, of size d+d^2+d^3"
+
+      p <- θ$p
+      β <- θ$β
+      λ <- sqrt(colSums(β^2))
+      b <- θ$b
+
+      # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1
+      β2 <- apply(β, 2, function(col) col %o% col)
+      β3 <- apply(β, 2, function(col) col %o% col %o% col)
+
+      c(
+        β  %*% (p * .G(li,1,λ,b)),
+        β2 %*% (p * .G(li,2,λ,b)),
+        β3 %*% (p * .G(li,3,λ,b)))
+    },
+
+    f = function(θ)
+    {
+      "Product t(hat_Mi - Mi) W (hat_Mi - Mi) with Mi(theta)"
+
+      L <- expArgs(θ)
+      A <- as.matrix(Mhat - Moments(L))
+      t(A) %*% W %*% A
+    },
+
+    grad_f = function(θ)
+    {
+      "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
+
+      L <- expArgs(θ)
+      -2 * t(grad_M(L)) %*% W %*% as.matrix(Mhat - Moments(L))
+    },
+
+    grad_M = function(θ)
+    {
+      "Gradient of the vector of moments, size (dim=)d+d^2+d^3 x K-1+K+d*K"
+
+      p <- θ$p
+      β <- θ$β
+      λ <- sqrt(colSums(β^2))
+      μ <- sweep(β, 2, λ, '/')
+      b <- θ$b
+
+      res <- matrix(nrow=nrow(W), ncol=0)
+
+      # Tensorial products β^2 = β2 and β^3 = β3 must be computed from current β1
+      β2 <- apply(β, 2, function(col) col %o% col)
+      β3 <- apply(β, 2, function(col) col %o% col %o% col)
+
+      # Some precomputations
+      G1 = .G(li,1,λ,b)
+      G2 = .G(li,2,λ,b)
+      G3 = .G(li,3,λ,b)
+      G4 = .G(li,4,λ,b)
+      G5 = .G(li,5,λ,b)
+
+      # Gradient on p: K-1 columns, dim rows
+      km1 = 1:(K-1)
+      res <- cbind(res, rbind(
+        sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K],
+        sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K],
+        sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ))
+
+      for (i in 1:d)
+      {
+        # i determines the derivated matrix dβ[2,3]
+
+        dβ_left <- sweep(β, 2, p * G3 * β[i,], '*')
+        dβ_right <- matrix(0, nrow=d, ncol=K)
+        block <- i
+        dβ_right[block,] <- dβ_right[block,] + 1
+        dβ <- dβ_left + sweep(dβ_right, 2,  p * G1, '*')
+
+        dβ2_left <- sweep(β2, 2, p * G4 * β[i,], '*')
+        dβ2_right <- do.call( rbind, lapply(1:d, function(j) {
+          sweep(dβ_right, 2, β[j,], '*')
+        }) )
+        block <- ((i-1)*d+1):(i*d)
+        dβ2_right[block,] <- dβ2_right[block,] + β
+        dβ2 <- dβ2_left + sweep(dβ2_right, 2, p * G2, '*')
+
+        dβ3_left <- sweep(β3, 2, p * G5 * β[i,], '*')
+        dβ3_right <- do.call( rbind, lapply(1:d, function(j) {
+          sweep(dβ2_right, 2, β[j,], '*')
+        }) )
+        block <- ((i-1)*d*d+1):(i*d*d)
+        dβ3_right[block,] <- dβ3_right[block,] + β2
+        dβ3 <- dβ3_left + sweep(dβ3_right, 2, p * G3, '*')
+
+        res <- cbind(res, rbind(dβ, dβ2, dβ3))
+      }
+
+      # Gradient on b
+      res <- cbind(res, rbind(
+        sweep(β,  2, p * G2, '*'),
+        sweep(β2, 2, p * G3, '*'),
+        sweep(β3, 2, p * G4, '*') ))
+
+      res
+    },
+
+    run = function(θ0)
+    {
+      "Run optimization from θ0 with solver..."
+
+      if (!is.list(θ0))
+        stop("θ0: list")
+      if (is.null(θ0$β))
+        stop("At least θ0$β must be provided")
+      if (!is.matrix(θ0$β) || any(is.na(θ0$β))
+        || nrow(θ0$β) != d || ncol(θ0$β) != K)
+      {
+        stop("θ0$β: matrix, no NA, nrow = d, ncol = K")
+      }
+      if (is.null(θ0$p))
+        θ0$p = rep(1/K, K-1)
+      else if (!is.numeric(θ0$p) || length(θ0$p) != K-1
+        || any(is.na(θ0$p)) || sum(θ0$p) > 1)
+      {
+        stop("θ0$p: length K-1, no NA, positive integers, sum to <= 1")
+      }
+      if (is.null(θ0$b))
+        θ0$b = rep(0, K)
+      else if (!is.numeric(θ0$b) || length(θ0$b) != K || any(is.na(θ0$b)))
+        stop("θ0$b: length K, no NA")
+
+      # (Re)Set W to identity, to allow several run from the same object
+      W <<- diag(d+d^2+d^3)
+
+      loopMax <- 2 #TODO: loopMax = 3 ? Seems not improving...
+      x_init <- linArgs(θ0)
+      for (loop in 1:loopMax)
+      {
+        op_res = constrOptim( x_init, .self$f, .self$grad_f,
+          ui=cbind(
+            rbind( rep(-1,K-1), diag(K-1) ),
+            matrix(0, nrow=K, ncol=(d+1)*K) ),
+          ci=c(-1,rep(0,K-1)) )
+        if (loop < loopMax) #avoid computing an extra W
+          W <<- computeW(expArgs(op_res$par))
+        x_init <- op_res$par
+        #print(op_res$value) #debug
+        #print(expArgs(op_res$par)) #debug
+      }
+
+      expArgs(op_res$par)
+    }
+  )
 )
 
 # Compute vectorial E[g^{(order)}(<β,x> + b)] with x~N(0,Id) (integral in R^d)
 #                 = E[g^{(order)}(z)] with z~N(b,diag(λ))
+# by numerically evaluating the integral.
 #
 # @param link Link, 'logit' or 'probit'
 # @param order Order of derivative
@@ -260,82 +300,49 @@ setRefClass(
 #
 .G <- function(link, order, λ, b)
 {
-	# NOTE: weird "integral divergent" error on inputs:
-	# link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969)
-	# Switch to pracma package for that (but it seems slow...)
-
-	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) )
+  # NOTE: weird "integral divergent" error on inputs:
+  # link="probit"; order=2; λ=c(531.8099,586.8893,523.5816); b=c(-118.512674,-3.488020,2.109969)
+  # Switch to pracma package for that (but it seems slow...)
+  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
+  })
 }
 
 # Derivatives list: g^(k)(x) for links 'logit' and 'probit'
 #
 .deriv <- list(
-	"probit"=list(
-		# 'probit' derivatives list;
-		# TODO: exact values for the integral E[g^(k)(λz+b)]
-		function(x) exp(-x^2/2)/(sqrt(2*pi)),                     #g'
-		function(x) exp(-x^2/2)/(sqrt(2*pi)) *  -x,               #g''
-		function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^2 - 1),        #g^(3)
-		function(x) exp(-x^2/2)/(sqrt(2*pi)) * (-x^3 + 3*x),      #g^(4)
-		function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^4 - 6*x^2 + 3) #g^(5)
-	),
-	"logit"=list(
-		# Sigmoid derivatives list, obtained with http://www.derivative-calculator.net/
-		# @seealso http://www.ece.uc.edu/~aminai/papers/minai_sigmoids_NN93.pdf
-		function(x) {e=exp(x); .zin(e                                    /(e+1)^2)}, #g'
-		function(x) {e=exp(x); .zin(e*(-e   + 1)                         /(e+1)^3)}, #g''
-		function(x) {e=exp(x); .zin(e*( e^2 - 4*e    + 1)                /(e+1)^4)}, #g^(3)
-		function(x) {e=exp(x); .zin(e*(-e^3 + 11*e^2 - 11*e   + 1)       /(e+1)^5)}, #g^(4)
-		function(x) {e=exp(x); .zin(e*( e^4 - 26*e^3 + 66*e^2 - 26*e + 1)/(e+1)^6)}  #g^(5)
-	)
+  "probit"=list(
+    # 'probit' derivatives list;
+    # NOTE: exact values for the integral E[g^(k)(λz+b)] could be computed
+    function(x) exp(-x^2/2)/(sqrt(2*pi)),                     #g'
+    function(x) exp(-x^2/2)/(sqrt(2*pi)) *  -x,               #g''
+    function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^2 - 1),        #g^(3)
+    function(x) exp(-x^2/2)/(sqrt(2*pi)) * (-x^3 + 3*x),      #g^(4)
+    function(x) exp(-x^2/2)/(sqrt(2*pi)) * ( x^4 - 6*x^2 + 3) #g^(5)
+  ),
+  "logit"=list(
+    # Sigmoid derivatives list, obtained with http://www.derivative-calculator.net/
+    # @seealso http://www.ece.uc.edu/~aminai/papers/minai_sigmoids_NN93.pdf
+    function(x) {e=exp(x); .zin(e                                    /(e+1)^2)}, #g'
+    function(x) {e=exp(x); .zin(e*(-e   + 1)                         /(e+1)^3)}, #g''
+    function(x) {e=exp(x); .zin(e*( e^2 - 4*e    + 1)                /(e+1)^4)}, #g^(3)
+    function(x) {e=exp(x); .zin(e*(-e^3 + 11*e^2 - 11*e   + 1)       /(e+1)^5)}, #g^(4)
+    function(x) {e=exp(x); .zin(e*( e^4 - 26*e^3 + 66*e^2 - 26*e + 1)/(e+1)^6)}  #g^(5)
+  )
 )
 
 # Utility for integration: "[return] zero if [argument is] NaN" (Inf / Inf divs)
@@ -344,6 +351,6 @@ setRefClass(
 #
 .zin <- function(x)
 {
-	x[is.nan(x)] <- 0.
-	x
+  x[is.nan(x)] <- 0.
+  x
 }