X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FoptimParams.R;h=a45f71ac4f64769f7da037324efef533ffa58ba0;hb=b357ac5457105ee87ad8650d61ad3c3845e0df95;hp=c061fcf0b7d8690e901b02d8ae4943c7c940701a;hpb=d08fef424150599b8095727c0f9870ca9535fb65;p=morpheus.git

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index c061fcf..a45f71a 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -31,12 +31,12 @@
 #' 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"))
 {
 	# 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)
   if (!is.numeric(K) || K!=floor(K) || K < 2)
@@ -86,11 +86,11 @@ setRefClass(
 				stop("Missing arguments")
 
       # Precompute empirical moments
-      M <- computeMoments(optargs$X,optargs$Y)
+      M <- computeMoments(X, Y)
       M1 <- as.double(M[[1]])
       M2 <- as.double(M[[2]])
       M3 <- as.double(M[[3]])
-      Mhat <<- matrix(c(M1,M2,M3), ncol=1)
+      Mhat <<- c(M1, M2, M3)
 
 			n <<- nrow(X)
 			d <<- length(M1)
@@ -117,15 +117,15 @@ setRefClass(
 
     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(θ)),
-        pn=as.integer(n), pd=as.integer(d),
-        W=as.double(W), PACKAGE="morpheus")$W, nrow=dim, ncol=dim) )
+      #require(MASS)
+      dd <- d + d^2 + d^3
+      W <<- MASS::ginv( matrix( .C("Compute_Omega",
+        X=as.double(X), Y=Y, M=Moments(θ), pn=as.integer(n), pd=as.integer(d),
+        W=as.double(W), PACKAGE="morpheus")$W, nrow=dd, ncol=dd ) )
       NULL #avoid returning W
     },
 
-    M <- function(θ)
+    Moments = function(θ)
     {
       "Vector of moments, of size d+d^2+d^3"
 
@@ -138,17 +138,18 @@ setRefClass(
 			β2 <- apply(β, 2, function(col) col %o% col)
 			β3 <- apply(β, 2, function(col) col %o% col %o% col)
 
-			matrix(c(
+			c(
 				β  %*% (p * .G(li,1,λ,b)),
 				β2 %*% (p * .G(li,2,λ,b)),
-				β3 %*% (p * .G(li,3,λ,b))), ncol=1)
+				β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
     },
 
@@ -156,19 +157,19 @@ setRefClass(
 		{
 			"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$β
+			p <- θ$p
+			β <- θ$β
 			λ <- sqrt(colSums(β^2))
 			μ <- sweep(β, 2, λ, '/')
-			b <- L$b
+			b <- θ$b
 
       res <- matrix(nrow=nrow(W), ncol=0)
 
@@ -186,11 +187,10 @@ setRefClass(
       # 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] )))
+        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] ))
 
-      # TODO: understand derivatives order and match the one in optim init param
 			for (i in 1:d)
 			{
 				# i determines the derivated matrix dβ[2,3]
@@ -217,14 +217,14 @@ setRefClass(
 				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)))
+				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, '*') )))
+				sweep(β,  2, p * G2, '*'),
+				sweep(β2, 2, p * G3, '*'),
+				sweep(β3, 2, p * G4, '*') ))
 
 			res
 		},
@@ -249,17 +249,21 @@ setRefClass(
       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 <= ε ?
+      # TODO: stopping condition? N iterations? Delta <= epsilon ?
+      for (loop in 1:10)
+      {
+        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)) )
+
+        computeW(expArgs(op_res$par))
+        # debug:
+        #print(W)
+        print(op_res$value)
+        print(expArgs(op_res$par))
+      }
 
 			expArgs(op_res$par)
 		}