mirror of
https://github.com/johnkerl/miller.git
synced 2026-07-18 00:45:47 +00:00
add covariance-matrix option to stats2
This commit is contained in:
John Kerl
parent
dd925cb914
commit
1f7c2a5b71
5 changed files with 136 additions and 84 deletions
|
|
@ -3,19 +3,31 @@
|
|||
|
||||
// xxx cmt intended for streaming applications. otherwise the formulas are
|
||||
// different (and more intuitive).
|
||||
double mlr_get_stddev(unsigned long long count, double sum, double sum2) {
|
||||
double mean = sum / count;
|
||||
double numerator = sum2 - 2.0*mean*sum + count*mean*mean;
|
||||
double mlr_get_stddev(unsigned long long n, double sum, double sum2) {
|
||||
double mean = sum / n;
|
||||
double numerator = sum2 - 2.0*mean*sum + n*mean*mean;
|
||||
if (numerator < 0.0) // round-off error
|
||||
numerator = 0.0;
|
||||
double denominator = count - 1LL;
|
||||
double denominator = n - 1LL;
|
||||
return sqrt(numerator / denominator);
|
||||
}
|
||||
|
||||
double mlr_get_cov(unsigned long long count, double sumx, double sumy, double sumxy) {
|
||||
double meanx = sumx / count;
|
||||
double meany = sumy / count;
|
||||
double numerator = sumxy - meanx * sumy - meany * sumx + count * meanx * meany;
|
||||
double denominator = count - 1;
|
||||
double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy) {
|
||||
double meanx = sumx / n;
|
||||
double meany = sumy / n;
|
||||
double numerator = sumxy - meanx*sumy - meany*sumx + n*meanx*meany;
|
||||
double denominator = n - 1;
|
||||
return numerator / denominator;
|
||||
}
|
||||
|
||||
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 denominator = n - 1;
|
||||
*pq00 = (sumx2 - sumx*sumx/n) / denominator;
|
||||
*pq01 = (sumxy - sumx*sumy/n) / denominator;
|
||||
*pq10 = *pq01;
|
||||
*pq11 = (sumy2 - sumy*sumy/n) / denominator;
|
||||
}
|
||||
|
||||
|
|
|
|||
|
|
@ -3,5 +3,8 @@
|
|||
|
||||
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);
|
||||
|
||||
#endif // MLRSTAT_H
|
||||
|
|
|
|||
|
|
@ -71,7 +71,6 @@ sllv_t* mapper_sort_func(lrec_t* pinrec, context_t* pctx, void* pvstate) {
|
|||
}
|
||||
}
|
||||
|
||||
// xxx move to mlrutil?
|
||||
static int string_list_compare(const void* pva, const void* pvb) {
|
||||
const lhmslve_t* pea = pva;
|
||||
const lhmslve_t* peb = pvb;
|
||||
|
|
|
|||
|
|
@ -13,6 +13,10 @@
|
|||
#include "mapping/mappers.h"
|
||||
#include "cli/argparse.h"
|
||||
|
||||
#define DO_CORR 0x11
|
||||
#define DO_COV 0x22
|
||||
#define DO_COVX 0x33
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef void stats2_put_func_t(void* pvstate, double x, double y);
|
||||
typedef void stats2_get_func_t(void* pvstate, char* name1, char* name2, lrec_t* poutrec);
|
||||
|
|
@ -132,79 +136,6 @@ stats2_t* stats2_linreg_alloc(static_context_t* pstatx) {
|
|||
return pstats2;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
// def find_sample_covariance(xs, ys):
|
||||
// N = len(xs)
|
||||
// mean_x = find_mean(xs)
|
||||
// mean_y = find_mean(ys)
|
||||
//
|
||||
// sum = 0.0
|
||||
// for k in range(0, N):
|
||||
// sum += (xs[k] - mean_x) * (ys[k] - mean_y)
|
||||
//
|
||||
// return sum / (N-1.0)
|
||||
|
||||
// Corr(X,Y) = Cov(X,Y) / sigma_X sigma_Y.
|
||||
typedef struct _stats2_corr_cov_state_t {
|
||||
unsigned long long count;
|
||||
double sumx;
|
||||
double sumy;
|
||||
double sumx2;
|
||||
double sumxy;
|
||||
double sumy2;
|
||||
int do_corr;
|
||||
static_context_t* pstatx;
|
||||
} stats2_corr_cov_state_t;
|
||||
void stats2_corr_cov_put(void* pvstate, double x, double y) {
|
||||
stats2_corr_cov_state_t* pstate = pvstate;
|
||||
pstate->count++;
|
||||
pstate->sumx += x;
|
||||
pstate->sumy += y;
|
||||
pstate->sumx2 += x*x;
|
||||
pstate->sumxy += x*y;
|
||||
pstate->sumy2 += y*y;
|
||||
}
|
||||
|
||||
void stats2_corr_cov_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
|
||||
stats2_corr_cov_state_t* pstate = pvstate;
|
||||
char* suffix = (pstate->do_corr) ? "corr" : "cov";
|
||||
char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
|
||||
if (pstate->count < 2LL) {
|
||||
lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);
|
||||
} else {
|
||||
double output = mlr_get_cov(pstate->count, pstate->sumx, pstate->sumy, pstate->sumxy);
|
||||
if (pstate->do_corr) {
|
||||
double sigmax = mlr_get_stddev(pstate->count, pstate->sumx, pstate->sumx2);
|
||||
double sigmay = mlr_get_stddev(pstate->count, pstate->sumy, pstate->sumy2);
|
||||
output = output / sigmax / sigmay;
|
||||
}
|
||||
char* val = mlr_alloc_string_from_double(output, pstate->pstatx->ofmt);
|
||||
lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
|
||||
}
|
||||
}
|
||||
stats2_t* stats2_corr_cov_alloc(int do_corr, static_context_t* pstatx) {
|
||||
stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_t));
|
||||
stats2_corr_cov_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_corr_cov_state_t));
|
||||
pstate->count = 0LL;
|
||||
pstate->sumx = 0.0;
|
||||
pstate->sumy = 0.0;
|
||||
pstate->sumx2 = 0.0;
|
||||
pstate->sumxy = 0.0;
|
||||
pstate->sumy2 = 0.0;
|
||||
pstate->do_corr = do_corr;
|
||||
pstate->pstatx = pstatx;
|
||||
pstats2->pvstate = (void*)pstate;
|
||||
pstats2->pput_func = &stats2_corr_cov_put;
|
||||
pstats2->pget_func = &stats2_corr_cov_get;
|
||||
return pstats2;
|
||||
}
|
||||
stats2_t* stats2_corr_alloc(static_context_t* pstatx) {
|
||||
return stats2_corr_cov_alloc(TRUE, pstatx);
|
||||
}
|
||||
stats2_t* stats2_cov_alloc(static_context_t* pstatx) {
|
||||
return stats2_corr_cov_alloc(FALSE, pstatx);
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
// http://en.wikipedia.org/wiki/Pearson_product-moment_correlation_coefficient
|
||||
// Alternatively, just use sqrt(corr) as defined above.
|
||||
|
|
@ -264,6 +195,108 @@ stats2_t* stats2_r2_alloc(static_context_t* pstatx) {
|
|||
return pstats2;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
// def find_sample_covariance(xs, ys):
|
||||
// N = len(xs)
|
||||
// mean_x = find_mean(xs)
|
||||
// mean_y = find_mean(ys)
|
||||
//
|
||||
// sum = 0.0
|
||||
// for k in range(0, N):
|
||||
// sum += (xs[k] - mean_x) * (ys[k] - mean_y)
|
||||
//
|
||||
// return sum / (N-1.0)
|
||||
|
||||
// Corr(X,Y) = Cov(X,Y) / sigma_X sigma_Y.
|
||||
typedef struct _stats2_corr_cov_state_t {
|
||||
unsigned long long count;
|
||||
double sumx;
|
||||
double sumy;
|
||||
double sumx2;
|
||||
double sumxy;
|
||||
double sumy2;
|
||||
int do_which;
|
||||
static_context_t* pstatx;
|
||||
} stats2_corr_cov_state_t;
|
||||
void stats2_corr_cov_put(void* pvstate, double x, double y) {
|
||||
stats2_corr_cov_state_t* pstate = pvstate;
|
||||
pstate->count++;
|
||||
pstate->sumx += x;
|
||||
pstate->sumy += y;
|
||||
pstate->sumx2 += x*x;
|
||||
pstate->sumxy += x*y;
|
||||
pstate->sumy2 += y*y;
|
||||
}
|
||||
|
||||
void stats2_corr_cov_get(void* pvstate, char* name1, char* name2, lrec_t* poutrec) {
|
||||
stats2_corr_cov_state_t* pstate = pvstate;
|
||||
if (pstate->do_which == DO_COVX) {
|
||||
char* key00 = mlr_paste_4_strings(name1, "_", name2, "_q00");
|
||||
char* key01 = mlr_paste_4_strings(name1, "_", name2, "_q01");
|
||||
char* key10 = mlr_paste_4_strings(name1, "_", name2, "_q10");
|
||||
char* key11 = mlr_paste_4_strings(name1, "_", name2, "_q11");
|
||||
if (pstate->count < 2LL) {
|
||||
lrec_put(poutrec, key00, "", LREC_FREE_ENTRY_KEY);
|
||||
lrec_put(poutrec, key01, "", LREC_FREE_ENTRY_KEY);
|
||||
lrec_put(poutrec, key10, "", LREC_FREE_ENTRY_KEY);
|
||||
lrec_put(poutrec, key11, "", LREC_FREE_ENTRY_KEY);
|
||||
} else {
|
||||
double q00, q01, q10, q11;
|
||||
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);
|
||||
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 {
|
||||
char* suffix = (pstate->do_which == DO_CORR) ? "corr" : "cov";
|
||||
char* key = mlr_paste_5_strings(name1, "_", name2, "_", suffix);
|
||||
if (pstate->count < 2LL) {
|
||||
lrec_put(poutrec, key, "", LREC_FREE_ENTRY_KEY);
|
||||
} else {
|
||||
double output = mlr_get_cov(pstate->count, pstate->sumx, pstate->sumy, pstate->sumxy);
|
||||
if (pstate->do_which == DO_CORR) {
|
||||
double sigmax = mlr_get_stddev(pstate->count, pstate->sumx, pstate->sumx2);
|
||||
double sigmay = mlr_get_stddev(pstate->count, pstate->sumy, pstate->sumy2);
|
||||
output = output / sigmax / sigmay;
|
||||
}
|
||||
char* val = mlr_alloc_string_from_double(output, pstate->pstatx->ofmt);
|
||||
lrec_put(poutrec, key, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
|
||||
}
|
||||
}
|
||||
}
|
||||
stats2_t* stats2_corr_cov_alloc(int do_which, static_context_t* pstatx) {
|
||||
stats2_t* pstats2 = mlr_malloc_or_die(sizeof(stats2_t));
|
||||
stats2_corr_cov_state_t* pstate = mlr_malloc_or_die(sizeof(stats2_corr_cov_state_t));
|
||||
pstate->count = 0LL;
|
||||
pstate->sumx = 0.0;
|
||||
pstate->sumy = 0.0;
|
||||
pstate->sumx2 = 0.0;
|
||||
pstate->sumxy = 0.0;
|
||||
pstate->sumy2 = 0.0;
|
||||
pstate->do_which = do_which;
|
||||
pstate->pstatx = pstatx;
|
||||
pstats2->pvstate = (void*)pstate;
|
||||
pstats2->pput_func = &stats2_corr_cov_put;
|
||||
pstats2->pget_func = &stats2_corr_cov_get;
|
||||
return pstats2;
|
||||
}
|
||||
stats2_t* stats2_corr_alloc(static_context_t* pstatx) {
|
||||
return stats2_corr_cov_alloc(DO_CORR, pstatx);
|
||||
}
|
||||
stats2_t* stats2_cov_alloc(static_context_t* pstatx) {
|
||||
return stats2_corr_cov_alloc(DO_COV, pstatx);
|
||||
}
|
||||
stats2_t* stats2_covx_alloc(static_context_t* pstatx) {
|
||||
return stats2_corr_cov_alloc(DO_COVX, pstatx);
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _stats2_lookup_t {
|
||||
char* name;
|
||||
|
|
@ -272,9 +305,10 @@ typedef struct _stats2_lookup_t {
|
|||
} 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},
|
||||
{"r2", stats2_r2_alloc},
|
||||
{"covx", stats2_covx_alloc},
|
||||
};
|
||||
static int stats2_lookup_table_length = sizeof(stats2_lookup_table) / sizeof(stats2_lookup_table[0]);
|
||||
|
||||
|
|
@ -486,3 +520,6 @@ mapper_setup_t mapper_stats2_setup = {
|
|||
.pusage_func = mapper_stats2_usage,
|
||||
.pparse_func = mapper_stats2_parse_cli
|
||||
};
|
||||
|
||||
// 1/(n-1) sumx2 - sumx**2 / n
|
||||
// 1/(n-1) sumxy - sumx*sumy / n
|
||||
|
|
|
|||
|
|
@ -8,6 +8,7 @@ NEATEN
|
|||
* mk perfcomp dir
|
||||
* "index-numbered" -> "implicitly index-numbered" in mlrwik
|
||||
* source hygiene: xxx's, top-of-header comments, readme re memory management, etc.
|
||||
* "level_1", "level_2", etc. are just confusing. come up with more descriptive names.
|
||||
|
||||
================================================================
|
||||
FEATURES
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue