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

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index 934a757..039070c 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -31,20 +31,32 @@
 #' o$f( o$linArgs(par0) )
 #' o$f( o$linArgs(par1) )
 #' @export
-optimParams = function(X, Y, K, link=c("logit","probit"))
+optimParams <- function(X, Y, K, link=c("logit","probit"), M=NULL)
 {
-	# Check arguments
+  # 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))
+  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)
+  link <- match.arg(link)
   if (!is.numeric(K) || K!=floor(K) || K < 2)
     stop("K: integer >= 2")
 
-	# Build and return optimization algorithm object
-	methods::new("OptimParams", "li"=link, "X"=X,
-    "Y"=as.integer(Y), "K"=as.integer(K))
+  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)
@@ -60,209 +72,221 @@ optimParams = function(X, Y, K, link=c("logit","probit"))
 #' @field W Weights matrix (iteratively refined)
 #'
 setRefClass(
-	Class = "OptimParams",
+  Class = "OptimParams",
 
-	fields = list(
-		# Inputs
-		li = "character", #link function
-		X = "matrix",
-		Y = "numeric",
+  fields = list(
+    # Inputs
+    li = "character", #link function
+    X = "matrix",
+    Y = "numeric",
     Mhat = "numeric", #vector of empirical moments
-		# Dimensions
-		K = "integer",
+    # Dimensions
+    K = "integer",
     n = "integer",
-		d = "integer",
+    d = "integer",
     # Weights matrix (generalized least square)
     W = "matrix"
-	),
+  ),
 
-	methods = list(
-		initialize = function(...)
-		{
-			"Check args and initialize K, d, W"
+  methods = list(
+    initialize = function(...)
+    {
+      "Check args and initialize K, d, W"
 
       callSuper(...)
-			if (!hasArg("X") || !hasArg("Y") || !hasArg("K") || !hasArg("li"))
-				stop("Missing arguments")
-
-      # Precompute empirical moments
-      M <- computeMoments(optargs$X,optargs$Y)
-      M1 <- as.double(M[[1]])
-      M2 <- as.double(M[[2]])
-      M3 <- as.double(M[[3]])
-      Mhat <<- matrix(c(M1,M2,M3), ncol=1)
-
-			n <<- nrow(X)
-			d <<- length(M1)
-      W <<- diag(d+d^2+d^3) #initialize at W = Identity
-		},
-
-		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)])),
-				"β" = matrix(v[K:(K+d*K-1)], ncol=K),
-				"b" = v[(K+d*K):(K+(d+1)*K-1)])
-		},
-
-		linArgs = function(L)
-		{
-			"Linearize vectors+matrices from list L into a vector"
-
-			c(L$p[1:(K-1)], as.double(L$β), L$b)
-		},
+      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(θ)
     {
-      dim <- d + d^2 + d^3
-      W <<- solve( matrix( .C("Compute_Omega",
-        X=as.double(X), Y=as.double(Y), M=as.double(M(θ)),
+      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=dim, ncol=dim) )
-      NULL #avoid returning W
+        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd )
+      MASS::ginv(Omega)
     },
 
-    M <- function(θ)
+    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)
-
-			matrix(c(
-				β  %*% (p * .G(li,1,λ,b)),
-				β2 %*% (p * .G(li,2,λ,b)),
-				β3 %*% (p * .G(li,3,λ,b))), ncol=1)
+      β <- θ$β
+      λ <- 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(Mi - hat_Mi) W (Mi - hat_Mi) with Mi(theta)"
+      "Product t(hat_Mi - Mi) W (hat_Mi - Mi) with Mi(theta)"
 
-			A <- M(θ) - Mhat
+      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"
+    grad_f = function(θ)
+    {
+      "Gradient of f, dimension (K-1) + d*K + K = (d+2)*K - 1"
 
-      -2 * t(grad_M(θ)) %*% getW(θ) %*% (Mhat - M(θ))
-    }
+      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"
 
-      L <- expArgs(θ)
-			p <- L$p
-			β <- L$β
-			λ <- sqrt(colSums(β^2))
-			μ <- sweep(β, 2, λ, '/')
-			b <- L$b
+      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)
+      # 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)
+      # 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(
-        t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ),
-        t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ),
-        t( 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(t(dβ), t(dβ2), t(dβ3)))
-			}
+      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(
-				t( sweep(β,  2, p * G2, '*') ),
-				t( sweep(β2, 2, p * G3, '*') ),
-				t( sweep(β3, 2, p * G4, '*') )))
+      res <- cbind(res, rbind(
+        sweep(β,  2, p * G2, '*'),
+        sweep(β2, 2, p * G3, '*'),
+        sweep(β3, 2, p * G4, '*') ))
 
-			res
-		},
+      res
+    },
 
-		run = function(θ0)
-		{
-			"Run optimization from θ0 with solver..."
+    run = function(θ0)
+    {
+      "Run optimization from θ0 with solver..."
 
-	    if (!is.list(θ0))
-		    stop("θ0: list")
+      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$β)) || ncol(θ0$β) != K)
-				stop("θ0$β: matrix, no NA, ncol == K")
+      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 (length(θ0$p) != K-1 || sum(θ0$p) > 1)
-        stop("θ0$p should contain positive integers and sum to < 1")
-      # Next test = heuristic to detect missing b (when matrix is called "beta")
-      if (is.null(θ0$b) || all(θ0$b == θ0$β))
+      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 (any(is.na(θ0$b)))
-        stop("θ0$b cannot have missing values")
-
-			op_res = constrOptim( linArgs(θ0), .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)) )
-
-      # debug:
-      print(computeW(expArgs(op_res$par)))
-      # We get a first non-trivial estimation of W
-      # TODO: loop, this redefine f, so that we can call constrOptim again...
-      # Stopping condition? N iterations? Delta <= ε ?
-
-			expArgs(op_res$par)
-		}
-	)
+      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)
@@ -276,9 +300,9 @@ 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...)
+  # 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({
@@ -301,24 +325,24 @@ setRefClass(
 # Derivatives list: g^(k)(x) for links 'logit' and 'probit'
 #
 .deriv <- list(
-	"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)
-	)
+  "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)
@@ -327,6 +351,6 @@ setRefClass(
 #
 .zin <- function(x)
 {
-	x[is.nan(x)] <- 0.
-	x
+  x[is.nan(x)] <- 0.
+  x
 }