Remove weights for now; planning generalization with matrix W
[morpheus.git] / pkg / R / optimParams.R
index 2eada8f..f62e75a 100644 (file)
@@ -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