b634c98ddab1c8713171b1437e8597bbf98e5a11
3 // Index matrix (by columns)
4 #define mi(i, j, d1, d2) (j*d1 + i)
6 // Index 3-tensor (by columns, matrices ordered by last dim)
7 #define ti(i, j, k, d1, d2, d3) (k*d1*d2 + j*d1 + i)
9 // Empirical cross-moment of order 2 between X size nxd and Y size n
10 void Moments_M2(double* X
, double* Y
, int* pn
, int* pd
, double* M2
)
13 //double* M2 = (double*)calloc(d*d,sizeof(double));
15 // M2 = E[Y*X^*2] - E[Y*e^*2] = E[Y (X^*2 - I)]
16 for (int j
=0; j
<d
; j
++)
18 for (int i
=0; i
<n
; i
++)
20 M2
[mi(j
,j
,d
,d
)] -= Y
[i
] / n
;
21 for (int k
=0; k
<d
; k
++)
22 M2
[mi(j
,k
,d
,d
)] += Y
[i
] * X
[mi(i
,j
,n
,d
)]*X
[mi(i
,k
,n
,d
)] / n
;
27 // Empirical cross-moment of order 3 between X size nxd and Y size n
28 void Moments_M3(double* X
, double* Y
, int* pn
, int* pd
, double* M3
)
31 //double* M3 = (double*)calloc(d*d*d,sizeof(double));
33 // M3 = E[Y*X^*3] - E[Y*e*X*e] - E[Y*e*e*X] - E[Y*X*e*e]
34 for (int j
=0; j
<d
; j
++)
36 for (int k
=0; k
<d
; k
++)
38 for (int i
=0; i
<n
; i
++)
40 double tensor_elt
= Y
[i
]*X
[mi(i
,k
,n
,d
)] / n
;
41 M3
[ti(j
,k
,j
,d
,d
,d
)] -= tensor_elt
;
42 M3
[ti(j
,j
,k
,d
,d
,d
)] -= tensor_elt
;
43 M3
[ti(k
,j
,j
,d
,d
,d
)] -= tensor_elt
;
44 for (int o
=0; o
<d
; o
++)
45 M3
[ti(j
,k
,o
,d
,d
,d
)] += Y
[i
] * X
[mi(i
,j
,n
,d
)]*X
[mi(i
,k
,n
,d
)]*X
[mi(i
,o
,n
,d
)] / n
;
51 // W = 1/N sum( t(g(Zi,theta)) g(Zi,theta) )
52 // with g(Zi, theta) = i-th contribution to all moments (size dim) - real moments
53 void Compute_Omega(double* X
, int* Y
, double* M
, int* pn
, int* pd
, double* W
)
56 int dim
= d
+ d
*d
+ d
*d
*d
;
57 //double* W = (double*)malloc(dim*dim*sizeof(double));
60 for (int j
=0; j
<dim
; j
++)
62 for (int k
=0; k
<dim
; k
++)
65 double* g
= (double*)malloc(dim
*sizeof(double));
66 for (int i
=0; i
<n
; i
++)
69 for (int j
=0; j
<d
; j
++)
70 g
[j
] = Y
[i
] * X
[mi(i
,j
,n
,d
)] - M
[j
];
71 for (int j
=d
; j
<d
+(d
*d
); j
++)
73 int idx1
= (j
-d
) % d
; //num row
74 int idx2
= ((j
-d
) - idx1
) / d
; //num col
78 g
[j
] += Y
[i
] * X
[mi(i
,idx1
,n
,d
)]*X
[mi(i
,idx2
,n
,d
)] - M
[j
];
80 for (int j
=d
+d
*d
; j
<dim
; j
++)
82 int idx1
= (j
-d
-d
*d
) % d
; //num row
83 int idx2
= ((j
-d
-d
*d
- idx1
) / d
) %d
; //num col
84 int idx3
= (((j
-d
-d
*d
- idx1
) / d
) - idx2
) / d
; //num "depth"
87 g
[j
] -= Y
[i
] * X
[mi(i
,idx3
,n
,d
)];
89 g
[j
] -= Y
[i
] * X
[mi(i
,idx2
,n
,d
)];
91 g
[j
] -= Y
[i
] * X
[mi(i
,idx1
,n
,d
)];
92 g
[j
] += Y
[i
] * X
[mi(i
,idx1
,n
,d
)]*X
[mi(i
,idx2
,n
,d
)]*X
[mi(i
,idx3
,n
,d
)] - M
[j
];
94 // Add 1/n t(gi) %*% gi to W
95 for (int j
=0; j
<dim
; j
++)
97 for (int k
=0; k
<dim
; k
++)
98 W
[j
*dim
+k
] += g
[j
] * g
[k
] / n
;