From: Benjamin Auder Date: Mon, 9 Dec 2019 13:48:50 +0000 (+0100) Subject: Completed draft for version with matrix W X-Git-Url: https://git.auder.net/doc/html/css/img/current/pieces/config.php?a=commitdiff_plain;h=d08fef424150599b8095727c0f9870ca9535fb65;p=morpheus.git Completed draft for version with matrix W --- diff --git a/pkg/R/optimParams.R b/pkg/R/optimParams.R index 934a757..c061fcf 100644 --- a/pkg/R/optimParams.R +++ b/pkg/R/optimParams.R @@ -190,6 +190,7 @@ setRefClass( 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] ))) + # 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] diff --git a/pkg/R/utils.R b/pkg/R/utils.R index 6ac9bec..6d1c361 100644 --- a/pkg/R/utils.R +++ b/pkg/R/utils.R @@ -74,23 +74,6 @@ normalize = function(X) computeMoments = function(X, Y) list( colMeans(Y * X), .Moments_M2(X,Y), .Moments_M3(X,Y) ) -# Computes the Omega matrix for generalized least square method -# -# @param X matrix of covariates (of size n*d) -# @param Y vector of responses (of size n) -# @param theta list with p, beta, b -# -# @return Matrix of size dimxdim where dim=d+d^2+d^3 -# -.Moments_M3 = function(X, Y) -{ - n = nrow(X) - d = ncol(X) - M3 = array(0,dim=c(d,d,d)) - array( .C("Moments_M3", X=as.double(X), Y=as.double(Y), pn=as.integer(n), - pd=as.integer(d), M3=as.double(M3), PACKAGE="morpheus")$M3, dim=c(d,d,d) ) -} - # Find the optimal assignment (permutation) between two sets (minimize cost) # # @param distances The distances matrix, in columns (distances[i,j] is distance between i diff --git a/pkg/src/functions.c b/pkg/src/functions.c index 41065bd..8731b9f 100644 --- a/pkg/src/functions.c +++ b/pkg/src/functions.c @@ -1,18 +1,18 @@ #include -//index matrix (by columns) +// Index matrix (by columns) int mi(int i, int j, int d1, int d2) { - return j*d1+i; + return j*d1 + i; } -//index 3-tensor (by columns, matrices ordered by last dim) +// Index 3-tensor (by columns, matrices ordered by last dim) int ti(int i, int j, int k, int d1, int d2, int d3) { return k*d1*d2 + j*d1 + i; } -// Emprical cross-moment of order 2 between X size nxd and Y size n +// Empirical cross-moment of order 2 between X size nxd and Y size n void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2) { int n=*pn, d=*pd; @@ -30,7 +30,7 @@ void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2) } } -// Emprical cross-moment of order 3 between X size nxd and Y size n +// Empirical cross-moment of order 3 between X size nxd and Y size n void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3) { int n=*pn, d=*pd; @@ -54,17 +54,49 @@ void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3) } } +// W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) ) +// with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments void Compute_Omega(double* X, double* Y, double* M, int* pn, int* pd, double* W) { int n=*pn, d=*pd; - //double* W = (double*)calloc(d+d*d+d*d*d,sizeof(double)); - - // TODO: formula 1/N sum( t(g(Zi,theta)) g(Zi,theta) ) - // = 1/N sum( t( (XiYi-...) - M[i] ) ( ... ) ) - // --> similar to Moments_M2 and M3 above - for (int j=0; j< + //int dim = d+d*d+d*d*d + //double* W = (double*)calloc(dim*dim,sizeof(double)); + double* g = (double*)malloc(dim * sizeof(double)); for (int i=0; i