X-Git-Url: https://git.auder.net/?a=blobdiff_plain;f=pkg%2FR%2FoptimParams.R;h=f62e75a54e23864c5773f9b123483af2270b988b;hb=e92d9d9dad0780d1a7eec27220c31254c20ab60d;hp=1ef1494e74841078859d7a4ff26eadffee98e029;hpb=cf673dee64ab0ce02ccaaba0fb63a9a4589ee8aa;p=morpheus.git

diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R
index 1ef1494..f62e75a 100644
--- a/pkg/R/optimParams.R
+++ b/pkg/R/optimParams.R
@@ -10,6 +10,7 @@
 #'     \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
 #'   }
 #'
 #' @return An object 'op' of class OptimParams, initialized so that \code{op$run(x0)}
@@ -81,7 +82,9 @@ setRefClass(
 		M3 = "numeric", #M3 easier to process as a vector
 		# Dimensions
 		K = "integer",
-		d = "integer"
+		d = "integer",
+    # Weights matrix (generalized least square)
+    W = "matrix"
 	),
 
 	methods = list(
@@ -97,6 +100,7 @@ setRefClass(
 			}
 
 			d <<- length(M1)
+      W <<- diag(d+d^2+d^3) #initialize at W = Identity
 		},
 
 		expArgs = function(x)
@@ -119,7 +123,7 @@ setRefClass(
 
 		f = function(x)
 		{
-			"Sum of squares (Mi - hat_Mi)^2 where Mi is obtained from formula"
+			"Product t(Mi - hat_Mi) W (Mi - hat_Mi) with Mi(theta)"
 
 			P <- expArgs(x)
 			p <- P$p
@@ -166,7 +170,7 @@ setRefClass(
 
 			km1 = 1:(K-1)
 			grad <- #gradient on p
-				t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 +
+			  t( sweep(as.matrix(β [,km1]), 2, G1[km1], '*') - G1[K] * β [,K] ) %*% F1 +
 				t( sweep(as.matrix(β2[,km1]), 2, G2[km1], '*') - G2[K] * β2[,K] ) %*% F2 +
 				t( sweep(as.matrix(β3[,km1]), 2, G3[km1], '*') - G3[K] * β3[,K] ) %*% F3
 
@@ -213,13 +217,23 @@ setRefClass(
 		{
 			"Run optimization from x0 with solver..."
 
-			if (!is.numeric(x0) || any(is.na(x0)) || length(x0) != (d+2)*K-1
-				|| any(x0[1:(K-1)] < 0) || sum(x0[1:(K-1)]) > 1)
-			{
-				stop("x0: numeric vector, no NA, length (d+2)*K-1, sum(x0[1:(K-1) >= 0]) <= 1")
-			}
-
-			op_res = constrOptim( x0, .self$f, .self$grad_f,
+	    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) ),