| 1 | #include <stdlib.h> |
| 2 | |
| 3 | //index matrix (by columns) |
| 4 | int mi(int i, int j, int d1, int d2) |
| 5 | { |
| 6 | return j*d1+i; |
| 7 | } |
| 8 | |
| 9 | //index 3-tensor (by columns, matrices ordered by last dim) |
| 10 | int ti(int i, int j, int k, int d1, int d2, int d3) |
| 11 | { |
| 12 | return k*d1*d2 + j*d1 + i; |
| 13 | } |
| 14 | |
| 15 | // Emprical cross-moment of order 2 between X size nxd and Y size n |
| 16 | void Moments_M2(double* X, double* Y, int* pn, int* pd, double* M2) |
| 17 | { |
| 18 | int n=*pn, d=*pd; |
| 19 | //double* M2 = (double*)calloc(d*d,sizeof(double)); |
| 20 | |
| 21 | // M2 = E[Y*X^*2] - E[Y*e^*2] = E[Y (X^*2 - I)] |
| 22 | for (int j=0; j<d; j++) |
| 23 | { |
| 24 | for (int i=0; i<n; i++) |
| 25 | { |
| 26 | M2[mi(j,j,d,d)] -= Y[i] / n; |
| 27 | for (int k=0; k<d; k++) |
| 28 | M2[mi(j,k,d,d)] += Y[i] * X[mi(i,j,n,d)]*X[mi(i,k,n,d)] / n; |
| 29 | } |
| 30 | } |
| 31 | } |
| 32 | |
| 33 | // Emprical cross-moment of order 3 between X size nxd and Y size n |
| 34 | void Moments_M3(double* X, double* Y, int* pn, int* pd, double* M3) |
| 35 | { |
| 36 | int n=*pn, d=*pd; |
| 37 | //double* M3 = (double*)calloc(d*d*d,sizeof(double)); |
| 38 | |
| 39 | // M3 = E[Y*X^*3] - E[Y*e*X*e] - E[Y*e*e*X] - E[Y*X*e*e] |
| 40 | for (int j=0; j<d; j++) |
| 41 | { |
| 42 | for (int k=0; k<d; k++) |
| 43 | { |
| 44 | for (int i=0; i<n; i++) |
| 45 | { |
| 46 | double tensor_elt = Y[i]*X[mi(i,k,n,d)] / n; |
| 47 | M3[ti(j,k,j,d,d,d)] -= tensor_elt; |
| 48 | M3[ti(j,j,k,d,d,d)] -= tensor_elt; |
| 49 | M3[ti(k,j,j,d,d,d)] -= tensor_elt; |
| 50 | for (int o=0; o<d; o++) |
| 51 | 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; |
| 52 | } |
| 53 | } |
| 54 | } |
| 55 | } |
| 56 | |
| 57 | void Compute_Omega(double* X, double* Y, int* pn, int* pd, double* W) |
| 58 | { |
| 59 | // TODO: formula 1/N sum( t(g(Zi,theta)) g(Zi,theta) ) |
| 60 | // = 1/N sum( t( (XiYi-...) - theta[i] ) ( ... ) ) |
| 61 | // --> similar to Moments_M2 and M3 above |
| 62 | } |