diff --git a/c/lib/mlrmath.c b/c/lib/mlrmath.c new file mode 100644 index 000000000..7aed6ca77 --- /dev/null +++ b/c/lib/mlrmath.c @@ -0,0 +1,146 @@ +#include +#include +#include +#include "mlrmath.h" + +#define JACOBI_TOLERANCE 1e-12 +#define JACOBI_MAXITER 20 + +static void matmul2(double C[2][2], double A[2][2], double B[2][2]); +static void matmul2t(double C[2][2], double A[2][2], double B[2][2]); + +// ---------------------------------------------------------------- +// Jacobi real-symmetric eigensolver. Loosely adapted from Numerical Recipes. +// +// Note: this is coded for n=2 (to implement PCA linear regression on 2 +// variables) but the algorithm is quite general. Changing from 2 to n is a +// matter of updating the top and bottom of the function: function signature to +// take double** matrix, double* eigenvector_1, double* eigenvector_2, and n; +// create copy-matrix and make-identity matrix functions; free temp matrices at +// the end; etc. + +void mlr_get_real_symmetric_eigensystem( + double matrix[2][2], // Input + double *peigenvalue_1, // Output: dominant eigenvalue + double *peigenvalue_2, // Output: less-dominant eigenvalue + double eigenvector_1[2], // Output: corresponding to dominant eigenvalue + double eigenvector_2[2]) // Output: corresponding to less-dominant eigenvalue +{ + double L[2][2] = { + { matrix[0][0], matrix[0][1] }, + { matrix[1][0], matrix[1][1] } + }; + double V[2][2] = { + { 1.0, 0.0 }, + { 0.0, 1.0 }, + }; + double P[2][2], PT_A[2][2]; + int n = 2; + + int found = 0; + for (int iter = 0; iter < JACOBI_MAXITER; iter++) { + double sum = 0.0; + for (int i = 1; i < n; i++) + for (int j = 0; j < i; j++) + sum += fabs(L[i][j]); + if (fabs(sum*sum) < JACOBI_TOLERANCE) { + found = 1; + break; + } + + for (int p = 0; p < n; p++) { + for (int q = p+1; q < n; q++) { + double numer = L[p][p] - L[q][q]; + double denom = L[p][q] + L[q][p]; + if (fabs(denom) < JACOBI_TOLERANCE) + continue; + double theta = numer / denom; + int sign_theta = (theta < 0) ? -1 : 1; + double t = sign_theta / (fabs(theta) + sqrt(theta*theta + 1)); + double c = 1.0 / sqrt(t*t + 1); + double s = t * c; + + for (int pi = 0; pi < n; pi++) + for (int pj = 0; pj < n; pj++) + P[pi][pj] = (pi == pj) ? 1.0 : 0.0; + P[p][p] = c; + P[p][q] = -s; + P[q][p] = s; + P[q][q] = c; + + // L = P.transpose() * L * P + // V = V * P + matmul2t(PT_A, P, L); + matmul2(L, PT_A, P); + matmul2(V, V, P); + } + } + } + + if (!found) { + fprintf(stderr, + "Jacobi eigensolver: max iterations (%d) exceeded. Non-symmetric input?\n", JACOBI_MAXITER); + exit(1); + } + + double eigenvalue_1 = L[0][0]; + double eigenvalue_2 = L[1][1]; + double abs1 = fabs(eigenvalue_1); + double abs2 = fabs(eigenvalue_2); + if (abs1 > abs2) { + *peigenvalue_1 = eigenvalue_1; + *peigenvalue_2 = eigenvalue_2; + eigenvector_1[0] = V[0][0]; // Column 0 of V + eigenvector_1[1] = V[1][0]; + eigenvector_2[0] = V[0][1]; // Column 1 of V + eigenvector_2[1] = V[1][1]; + } else { + *peigenvalue_1 = eigenvalue_2; + *peigenvalue_2 = eigenvalue_1; + eigenvector_1[0] = V[0][1]; + eigenvector_1[1] = V[1][1]; + eigenvector_2[0] = V[0][0]; + eigenvector_2[1] = V[1][0]; + } +} + +// xxx cmt mem-mgmt +static void matmul2( + double C[2][2], // Output + double A[2][2], // Input + double B[2][2]) // Input +{ + double T[2][2]; + int n = 2; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + double sum = 0.0; + for (int k = 0; k < n; k++) + sum += A[i][k] * B[k][j]; + T[i][j] = sum; + } + } + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + C[i][j] = T[i][j]; +} + +static void matmul2t( + double C[2][2], // Output + double A[2][2], // Input + double B[2][2]) // Input +{ + double T[2][2]; + int n = 2; + for (int i = 0; i < n; i++) { + for (int j = 0; j < n; j++) { + double sum = 0.0; + for (int k = 0; k < n; k++) + sum += A[k][i] * B[k][j]; + T[i][j] = sum; + } + } + for (int i = 0; i < n; i++) + for (int j = 0; j < n; j++) + C[i][j] = T[i][j]; +} diff --git a/c/lib/mlrmath.h b/c/lib/mlrmath.h new file mode 100644 index 000000000..fa6fb07df --- /dev/null +++ b/c/lib/mlrmath.h @@ -0,0 +1,11 @@ +#ifndef MLRMATH_H +#define MLRMATH_H + +void mlr_get_real_symmetric_eigensystem( + double matrix[2][2], // Input + double *peigenvalue_1, // Output: dominant eigenvalue + double *peigenvalue_2, // Output: less-dominant eigenvalue + double eigenvector_1[2], // Output: corresponding to dominant eigenvalue + double eigenvector_2[2]); // Output: corresponding to less-dominant eigenvalue + +#endif // MLRMATH_H diff --git a/c/lib/mlrstat.c b/c/lib/mlrstat.c index 114758ce9..9fed47293 100644 --- a/c/lib/mlrstat.c +++ b/c/lib/mlrstat.c @@ -22,12 +22,57 @@ double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy) void mlr_get_cov_matrix(unsigned long long n, double sumx, double sumx2, double sumy, double sumy2, double sumxy, - double* pq00, double* pq01, double* pq10, double* pq11) + double Q[2][2]) { double denominator = n - 1; - *pq00 = (sumx2 - sumx*sumx/n) / denominator; - *pq01 = (sumxy - sumx*sumy/n) / denominator; - *pq10 = *pq01; - *pq11 = (sumy2 - sumy*sumy/n) / denominator; + Q[0][0] = (sumx2 - sumx*sumx/n) / denominator; + Q[0][1] = (sumxy - sumx*sumy/n) / denominator; + Q[1][0] = Q[0][1]; + Q[1][1] = (sumy2 - sumy*sumy/n) / denominator; } +// ---------------------------------------------------------------- +// Principal component analysis can be used for linear regression: +// * Compute the covariance matrix for the x's and y's. +// * Find its eigenvalues and eigenvectors of the cov. (This is real-symmetric +// so Jacobi iteration is simple and fine.) +// * The principal eigenvector points in the direction of the fit. +// * The covariance matrix is computed on zero-mean data so the intercept +// is zero, of the form (y - nu) = m*(x - mu) where mu and nu are x and y +// means, respectively. +// * If the fit is perfect then the 2nd eigenvalue will be zero; if the fit is +// good then the 2nd eigenvalue will be smaller; if the fit is bad then +// they'll be about the same. I use 1 minus ratio of absolute values +// of 2nd to 1st eigenvalues as an indication of quality of the fit. +// +// Standard ("ordinary least-squares") linear regression is appropriate when +// the errors are thought to be all in the y's. PCA ("total least-squares") is +// appropriate when the x's and the y's are thought to both have errors. + +void mlr_get_linear_regression_pca( + // Inputs: + double eigenvalue_1, + double eigenvalue_2, + double eigenvector_1[2], + double eigenvector_2[2], + double x_mean, double y_mean, + // Outputs: + double* pm, double* pb, double* pquality) +{ + double abs_1 = fabs(eigenvalue_1); + double abs_2 = fabs(eigenvalue_2); + double quality = 1.0; + if (abs_1 == 0.0) + quality = 0.0; + else if (abs_2 > 0.0) + quality = 1.0 - abs_2 / abs_1; + + double a0 = eigenvector_1[0]; + double a1 = eigenvector_1[1]; + double m = a1 / a0; + double b = y_mean - m * x_mean; + + *pm = m; + *pb = b; + *pquality = quality; +} diff --git a/c/lib/mlrstat.h b/c/lib/mlrstat.h index 72e9cc981..13758bc5f 100644 --- a/c/lib/mlrstat.h +++ b/c/lib/mlrstat.h @@ -2,9 +2,20 @@ #define MLRSTAT_H double mlr_get_stddev(unsigned long long count, double sum, double sum2); + double mlr_get_cov(unsigned long long count, double sumx, double sumy, double sumxy); + void mlr_get_cov_matrix(unsigned long long n, - double sumx, double sumx2, double sumy, double sumy2, double sumxy, - double* pq00, double* pq01, double* pq10, double* pq11); + double sumx, double sumx2, double sumy, double sumy2, double sumxy, double Q[2][2]); + +void mlr_get_linear_regression_pca( + // Inputs: + double eigenvalue_1, + double eigenvalue_2, + double eigenvector_1[2], + double eigenvector_2[2], + double x_mean, double y_mean, + // Outputs, with quality 1 being a tight fit and quality 0 being a loose one. + double* pm, double* pb, double* pquality); #endif // MLRSTAT_H diff --git a/c/mapping/mapper_stats2.c b/c/mapping/mapper_stats2.c index c006bb786..b41877722 100644 --- a/c/mapping/mapper_stats2.c +++ b/c/mapping/mapper_stats2.c @@ -3,6 +3,7 @@ #include #include #include "lib/mlrutil.h" +#include "lib/mlrmath.h" #include "lib/mlrstat.h" #include "containers/sllv.h" #include "containers/slls.h" @@ -13,9 +14,10 @@ #include "mapping/mappers.h" #include "cli/argparse.h" -#define DO_CORR 0x11 -#define DO_COV 0x22 -#define DO_COVX 0x33 +#define DO_CORR 0x11 +#define DO_COV 0x22 +#define DO_COVX 0x33 +#define DO_LINREG_PCA 0x44 // ---------------------------------------------------------------- typedef void stats2_put_func_t(void* pvstate, double x, double y); @@ -29,6 +31,8 @@ typedef struct _stats2_t { typedef stats2_t* stats2_alloc_func_t(static_context_t* pstatx); +// xxx move to mlrstat.h/c + // ---------------------------------------------------------------- // Univariate linear regression // ---------------------------------------------------------------- @@ -216,6 +220,7 @@ typedef struct _stats2_corr_cov_state_t { double sumxy; double sumy2; int do_which; + // xxx do_verbose; static_context_t* pstatx; } stats2_corr_cov_state_t; void stats2_corr_cov_put(void* pvstate, double x, double y) { @@ -241,19 +246,62 @@ void stats2_corr_cov_get(void* pvstate, char* name1, char* name2, lrec_t* poutre lrec_put(poutrec, key10, "", LREC_FREE_ENTRY_KEY); lrec_put(poutrec, key11, "", LREC_FREE_ENTRY_KEY); } else { - double q00, q01, q10, q11; + double Q[2][2]; mlr_get_cov_matrix(pstate->count, - pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, - &q00, &q01, &q10, &q11); - char* val00 = mlr_alloc_string_from_double(q00, pstate->pstatx->ofmt); - char* val01 = mlr_alloc_string_from_double(q01, pstate->pstatx->ofmt); - char* val10 = mlr_alloc_string_from_double(q10, pstate->pstatx->ofmt); - char* val11 = mlr_alloc_string_from_double(q11, pstate->pstatx->ofmt); + pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, Q); + char* val00 = mlr_alloc_string_from_double(Q[0][0], pstate->pstatx->ofmt); + char* val01 = mlr_alloc_string_from_double(Q[0][1], pstate->pstatx->ofmt); + char* val10 = mlr_alloc_string_from_double(Q[1][0], pstate->pstatx->ofmt); + char* val11 = mlr_alloc_string_from_double(Q[1][1], pstate->pstatx->ofmt); lrec_put(poutrec, key00, val00, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); lrec_put(poutrec, key01, val01, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); lrec_put(poutrec, key10, val10, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); lrec_put(poutrec, key11, val11, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); } + } else if (pstate->do_which == DO_LINREG_PCA) { + char* keym = mlr_paste_4_strings(name1, "_", name1, "_pca_m"); + char* keyb = mlr_paste_4_strings(name1, "_", name2, "_pca_b"); + char* keyq = mlr_paste_4_strings(name2, "_", name1, "_pca_quality"); + char* keyl1 = mlr_paste_4_strings(name2, "_", name1, "_pca_eival1"); + char* keyl2 = mlr_paste_4_strings(name2, "_", name1, "_pca_eival2"); + char* keyv11 = mlr_paste_4_strings(name2, "_", name1, "_pca_eivec11"); + char* keyv12 = mlr_paste_4_strings(name2, "_", name1, "_pca_eivec12"); + char* keyv21 = mlr_paste_4_strings(name2, "_", name1, "_pca_eivec21"); + char* keyv22 = mlr_paste_4_strings(name2, "_", name1, "_pca_eivec22"); + if (pstate->count < 2LL) { + lrec_put(poutrec, keym, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyb, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyq, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyl1, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyl2, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyv11, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyv12, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyv21, "", LREC_FREE_ENTRY_KEY); + lrec_put(poutrec, keyv22, "", LREC_FREE_ENTRY_KEY); + } else { + double Q[2][2]; + mlr_get_cov_matrix(pstate->count, + pstate->sumx, pstate->sumx2, pstate->sumy, pstate->sumy2, pstate->sumxy, Q); + + double l1, l2; // Eigenvalues + double v1[2], v2[2]; // Eigenvectors + mlr_get_real_symmetric_eigensystem(Q, &l1, &l2, v1, v2); + + double x_mean = pstate->sumx / pstate->count; + double y_mean = pstate->sumy / pstate->count; + double m, b, q; + mlr_get_linear_regression_pca(l1, l2, v1, v2, x_mean, y_mean, &m, &b, &q); + + lrec_put(poutrec, keym, mlr_alloc_string_from_double(m, pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyb, mlr_alloc_string_from_double(b, pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyq, mlr_alloc_string_from_double(q, pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyl1, mlr_alloc_string_from_double(l1, pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyl2, mlr_alloc_string_from_double(l2, pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyv11, mlr_alloc_string_from_double(v1[0], pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyv12, mlr_alloc_string_from_double(v1[1], pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyv21, mlr_alloc_string_from_double(v2[0], pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + lrec_put(poutrec, keyv22, mlr_alloc_string_from_double(v2[1], pstate->pstatx->ofmt), LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE); + } } else { char* suffix = (pstate->do_which == DO_CORR) ? "corr" : "cov"; char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix); @@ -296,6 +344,9 @@ stats2_t* stats2_cov_alloc(static_context_t* pstatx) { stats2_t* stats2_covx_alloc(static_context_t* pstatx) { return stats2_corr_cov_alloc(DO_COVX, pstatx); } +stats2_t* stats2_linreg_pca_alloc(static_context_t* pstatx) { + return stats2_corr_cov_alloc(DO_LINREG_PCA, pstatx); +} // ---------------------------------------------------------------- typedef struct _stats2_lookup_t { @@ -304,11 +355,12 @@ typedef struct _stats2_lookup_t { static_context_t* pstatx; } stats2_lookup_t; static stats2_lookup_t stats2_lookup_table[] = { - {"linreg", stats2_linreg_alloc}, - {"r2", stats2_r2_alloc}, - {"corr", stats2_corr_alloc}, - {"cov", stats2_cov_alloc}, - {"covx", stats2_covx_alloc}, + {"linreg", stats2_linreg_alloc}, + {"r2", stats2_r2_alloc}, + {"corr", stats2_corr_alloc}, + {"cov", stats2_cov_alloc}, + {"covx", stats2_covx_alloc}, + {"linregpca", stats2_linreg_pca_alloc}, }; static int stats2_lookup_table_length = sizeof(stats2_lookup_table) / sizeof(stats2_lookup_table[0]);