mirror of
https://github.com/johnkerl/miller.git
synced 2026-08-01 20:12:19 +00:00
3rd/4th moments
This commit is contained in:
John Kerl
parent
b21e7cb979
commit
aad4ef13d2
4 changed files with 171 additions and 18 deletions
|
|
@ -81,13 +81,68 @@ void mlr_get_linear_regression_ols(unsigned long long n, double sumx, double sum
|
|||
// ----------------------------------------------------------------
|
||||
double mlr_get_var(unsigned long long n, double sum, double sum2) {
|
||||
double mean = sum / n;
|
||||
double numerator = sum2 - 2.0*mean*sum + n*mean*mean;
|
||||
double numerator = sum2 - mean*(2.0*sum - n*mean);
|
||||
if (numerator < 0.0) // round-off error
|
||||
numerator = 0.0;
|
||||
double denominator = n - 1LL;
|
||||
return numerator / denominator;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
// Unbiased estimator:
|
||||
// (1/n) sum{(xi-mean)**3}
|
||||
// -----------------------------
|
||||
// [(1/(n-1)) sum{(xi-mean)**2}]**1.5
|
||||
|
||||
// mean = sumx / n; n mean = sumx
|
||||
|
||||
// sum{(xi-mean)^3}
|
||||
// = sum{xi^3 - 3 mean xi^2 + 3 mean^2 xi - mean^3}
|
||||
// = sum{xi^3} - 3 mean sum{xi^2} + 3 mean^2 sum{xi} - n mean^3
|
||||
// = sumx3 - 3 mean sumx2 + 3 mean^2 sumx - n mean^3
|
||||
// = sumx3 - 3 mean sumx2 + 3n mean^3 - n mean^3
|
||||
// = sumx3 - 3 mean sumx2 + 2n mean^3
|
||||
// = sumx3 - mean*(3 sumx2 + 2n mean^2)
|
||||
|
||||
// sum{(xi-mean)^2}
|
||||
// = sum{xi^2 - 2 mean xi + mean^2}
|
||||
// = sum{xi^2} - 2 mean sum{xi} + n mean^2
|
||||
// = sumx2 - 2 mean sumx + n mean^2
|
||||
// = sumx2 - 2 n mean^2 + n mean^2
|
||||
// = sumx2 - n mean^2
|
||||
|
||||
double mlr_get_skewness(unsigned long long n, double sumx, double sumx2, double sumx3) {
|
||||
double mean = sumx / n;
|
||||
double numerator = sumx3 - mean*(3*sumx2 - 2*n*mean*mean);
|
||||
numerator = numerator / n;
|
||||
double denominator = (sumx2 - n*mean*mean) / (n-1);
|
||||
denominator = pow(denominator, 1.5);
|
||||
return numerator / denominator;
|
||||
}
|
||||
|
||||
// Unbiased:
|
||||
// (1/n) sum{(x-mean)**4}
|
||||
// ----------------------- - 3
|
||||
// [(1/n) sum{(x-mean)**2}]**2
|
||||
|
||||
// sum{(xi-mean)^4}
|
||||
// = sum{xi^4 - 4 mean xi^3 + 6 mean^2 xi^2 - 4 mean^3 xi + mean^4}
|
||||
// = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 4 mean^3 sum{xi} + n mean^4
|
||||
// = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 4 n mean^4 + n mean^4
|
||||
// = sum{xi^4} - 4 mean sum{xi^3} + 6 mean^2 sum{xi^2} - 3 n mean^4
|
||||
// = sum{xi^4} - mean*(4 sum{xi^3} - 6 mean sum{xi^2} + 3 n mean^3)
|
||||
// = sumx4 - mean*(4 sumx3 - 6 mean sumx2 + 3 n mean^3)
|
||||
// = sumx4 - mean*(4 sumx3 - mean*(6 sumx2 - 3 n mean^2))
|
||||
|
||||
double mlr_get_kurtosis(unsigned long long n, double sumx, double sumx2, double sumx3, double sumx4) {
|
||||
double mean = sumx / n;
|
||||
double numerator = sumx4 - mean*(4*sumx3 - mean*(6*sumx2 - 3*n*mean*mean));
|
||||
numerator = numerator / n;
|
||||
double denominator = (sumx2 - n*mean*mean) / n;
|
||||
denominator = denominator * denominator;
|
||||
return numerator / denominator - 3.0;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
// Non-streaming implementation:
|
||||
//
|
||||
|
|
|
|||
|
|
@ -8,6 +8,10 @@ double mlr_get_var(unsigned long long n, double sum, double sum2);
|
|||
|
||||
double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy);
|
||||
|
||||
double mlr_get_skewness(unsigned long long n, double sumx, double sumx2, double sumx3);
|
||||
|
||||
double mlr_get_kurtosis(unsigned long long n, double sumx, double sumx2, double sumx3, double sumx4);
|
||||
|
||||
void mlr_get_cov_matrix(unsigned long long n,
|
||||
double sumx, double sumx2, double sumy, double sumy2, double sumxy, double Q[2][2]);
|
||||
|
||||
|
|
|
|||
|
|
@ -63,6 +63,8 @@ static stats1_t* stats1_stddev_var_meaneb_alloc(char* value_field_name, char* st
|
|||
static stats1_t* stats1_stddev_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_var_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_meaneb_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_skewness_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_kurtosis_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_min_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_max_alloc(char* value_field_name, char* stats1_name);
|
||||
static stats1_t* stats1_percentile_alloc(char* value_field_name, char* stats1_name);
|
||||
|
|
@ -78,15 +80,17 @@ typedef struct _acc_lookup_t {
|
|||
char* desc;
|
||||
} stats1_lookup_t;
|
||||
static stats1_lookup_t stats1_lookup_table[] = {
|
||||
{"count", stats1_count_alloc, "Count instances of fields"},
|
||||
{"mode", stats1_mode_alloc, "Find most-frequently-occurring values for fields; first-found wins tie"},
|
||||
{"sum", stats1_sum_alloc, "Compute sums of specified fields"},
|
||||
{"mean", stats1_mean_alloc, "Compute averages (sample means) of specified fields"},
|
||||
{"stddev", stats1_stddev_alloc, "Compute sample standard deviation of specified fields"},
|
||||
{"var", stats1_var_alloc, "Compute sample variance of specified fields"},
|
||||
{"meaneb", stats1_meaneb_alloc, "Estimate error bars for averages (assuming no sample autocorrelation)"},
|
||||
{"min", stats1_min_alloc, "Compute minimum values of specified fields"},
|
||||
{"max", stats1_max_alloc, "Compute maximum values of specified fields"},
|
||||
{"count", stats1_count_alloc, "Count instances of fields"},
|
||||
{"mode", stats1_mode_alloc, "Find most-frequently-occurring values for fields; first-found wins tie"},
|
||||
{"sum", stats1_sum_alloc, "Compute sums of specified fields"},
|
||||
{"mean", stats1_mean_alloc, "Compute averages (sample means) of specified fields"},
|
||||
{"stddev", stats1_stddev_alloc, "Compute sample standard deviation of specified fields"},
|
||||
{"var", stats1_var_alloc, "Compute sample variance of specified fields"},
|
||||
{"meaneb", stats1_meaneb_alloc, "Estimate error bars for averages (assuming no sample autocorrelation)"},
|
||||
{"skewness", stats1_skewness_alloc, "Compute sample skewness of specified fields"},
|
||||
{"kurtosis", stats1_kurtosis_alloc, "Compute sample kurtosis of specified fields"},
|
||||
{"min", stats1_min_alloc, "Compute minimum values of specified fields"},
|
||||
{"max", stats1_max_alloc, "Compute maximum values of specified fields"},
|
||||
};
|
||||
static int stats1_lookup_table_length = sizeof(stats1_lookup_table) / sizeof(stats1_lookup_table[0]);
|
||||
|
||||
|
|
@ -432,7 +436,7 @@ static lrec_t* mapper_stats1_emit(mapper_stats1_state_t* pstate, lrec_t* poutrec
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_count_state_t {
|
||||
typedef struct _stats1_count_state_t {
|
||||
unsigned long long count;
|
||||
char* output_field_name;
|
||||
} stats1_count_state_t;
|
||||
|
|
@ -460,7 +464,7 @@ static stats1_t* stats1_count_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_mode_state_t {
|
||||
typedef struct _stats1_mode_state_t {
|
||||
lhmsi_t* pcounts_for_value;
|
||||
char* output_field_name;
|
||||
} stats1_mode_state_t;
|
||||
|
|
@ -502,7 +506,7 @@ static stats1_t* stats1_mode_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_sum_state_t {
|
||||
typedef struct _stats1_sum_state_t {
|
||||
double sum;
|
||||
char* output_field_name;
|
||||
} stats1_sum_state_t;
|
||||
|
|
@ -529,7 +533,7 @@ static stats1_t* stats1_sum_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_mean_state_t {
|
||||
typedef struct _stats1_mean_state_t {
|
||||
double sum;
|
||||
unsigned long long count;
|
||||
char* output_field_name;
|
||||
|
|
@ -566,7 +570,7 @@ static stats1_t* stats1_mean_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_stddev_var_meaneb_state_t {
|
||||
typedef struct _stats1_stddev_var_meaneb_state_t {
|
||||
unsigned long long count;
|
||||
double sumx;
|
||||
double sumx2;
|
||||
|
|
@ -621,7 +625,95 @@ static stats1_t* stats1_meaneb_alloc(char* value_field_name, char* stats1_name)
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_min_state_t {
|
||||
typedef struct _stats1_skewness_state_t {
|
||||
unsigned long long count;
|
||||
double sumx;
|
||||
double sumx2;
|
||||
double sumx3;
|
||||
char* output_field_name;
|
||||
} stats1_skewness_state_t;
|
||||
static void stats1_skewness_dingest(void* pvstate, double val) {
|
||||
stats1_skewness_state_t* pstate = pvstate;
|
||||
pstate->count++;
|
||||
pstate->sumx += val;
|
||||
pstate->sumx2 += val*val;
|
||||
pstate->sumx3 += val*val*val;
|
||||
}
|
||||
|
||||
static void stats1_skewness_emit(void* pvstate, char* value_field_name, char* stats1_name, lrec_t* poutrec) {
|
||||
stats1_skewness_state_t* pstate = pvstate;
|
||||
if (pstate->count < 2LL) {
|
||||
lrec_put(poutrec, pstate->output_field_name, "", LREC_FREE_ENTRY_KEY);
|
||||
} else {
|
||||
double output = mlr_get_skewness(pstate->count, pstate->sumx, pstate->sumx2, pstate->sumx3);
|
||||
char* val = mlr_alloc_string_from_double(output, MLR_GLOBALS.ofmt);
|
||||
lrec_put(poutrec, pstate->output_field_name, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
|
||||
}
|
||||
}
|
||||
|
||||
static stats1_t* stats1_skewness_alloc(char* value_field_name, char* stats1_name) {
|
||||
stats1_t* pstats1 = mlr_malloc_or_die(sizeof(stats1_t));
|
||||
stats1_skewness_state_t* pstate = mlr_malloc_or_die(sizeof(stats1_skewness_state_t));
|
||||
pstate->count = 0LL;
|
||||
pstate->sumx = 0.0;
|
||||
pstate->sumx2 = 0.0;
|
||||
pstate->sumx3 = 0.0;
|
||||
pstate->output_field_name = mlr_paste_3_strings(value_field_name, "_", stats1_name);
|
||||
|
||||
pstats1->pvstate = (void*)pstate;
|
||||
pstats1->psingest_func = NULL;
|
||||
pstats1->pdingest_func = stats1_skewness_dingest;
|
||||
pstats1->pemit_func = stats1_skewness_emit;
|
||||
return pstats1;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _stats1_kurtosis_state_t {
|
||||
unsigned long long count;
|
||||
double sumx;
|
||||
double sumx2;
|
||||
double sumx3;
|
||||
double sumx4;
|
||||
char* output_field_name;
|
||||
} stats1_kurtosis_state_t;
|
||||
static void stats1_kurtosis_dingest(void* pvstate, double val) {
|
||||
stats1_kurtosis_state_t* pstate = pvstate;
|
||||
pstate->count++;
|
||||
pstate->sumx += val;
|
||||
pstate->sumx2 += val*val;
|
||||
pstate->sumx3 += val*val*val;
|
||||
pstate->sumx4 += val*val*val*val;
|
||||
}
|
||||
|
||||
static void stats1_kurtosis_emit(void* pvstate, char* value_field_name, char* stats1_name, lrec_t* poutrec) {
|
||||
stats1_kurtosis_state_t* pstate = pvstate;
|
||||
if (pstate->count < 2LL) {
|
||||
lrec_put(poutrec, pstate->output_field_name, "", LREC_FREE_ENTRY_KEY);
|
||||
} else {
|
||||
double output = mlr_get_kurtosis(pstate->count, pstate->sumx, pstate->sumx2, pstate->sumx3, pstate->sumx4);
|
||||
char* val = mlr_alloc_string_from_double(output, MLR_GLOBALS.ofmt);
|
||||
lrec_put(poutrec, pstate->output_field_name, val, LREC_FREE_ENTRY_KEY|LREC_FREE_ENTRY_VALUE);
|
||||
}
|
||||
}
|
||||
|
||||
static stats1_t* stats1_kurtosis_alloc(char* value_field_name, char* stats1_name) {
|
||||
stats1_t* pstats1 = mlr_malloc_or_die(sizeof(stats1_t));
|
||||
stats1_kurtosis_state_t* pstate = mlr_malloc_or_die(sizeof(stats1_kurtosis_state_t));
|
||||
pstate->count = 0LL;
|
||||
pstate->sumx = 0.0;
|
||||
pstate->sumx2 = 0.0;
|
||||
pstate->sumx3 = 0.0;
|
||||
pstate->output_field_name = mlr_paste_3_strings(value_field_name, "_", stats1_name);
|
||||
|
||||
pstats1->pvstate = (void*)pstate;
|
||||
pstats1->psingest_func = NULL;
|
||||
pstats1->pdingest_func = stats1_kurtosis_dingest;
|
||||
pstats1->pemit_func = stats1_kurtosis_emit;
|
||||
return pstats1;
|
||||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _stats1_min_state_t {
|
||||
int have_min;
|
||||
double min;
|
||||
char* output_field_name;
|
||||
|
|
@ -660,7 +752,7 @@ static stats1_t* stats1_min_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_max_state_t {
|
||||
typedef struct _stats1_max_state_t {
|
||||
int have_max;
|
||||
double max;
|
||||
char* output_field_name;
|
||||
|
|
@ -699,7 +791,7 @@ static stats1_t* stats1_max_alloc(char* value_field_name, char* stats1_name) {
|
|||
}
|
||||
|
||||
// ----------------------------------------------------------------
|
||||
typedef struct _acc_percentile_state_t {
|
||||
typedef struct _stats1_percentile_state_t {
|
||||
percentile_keeper_t* ppercentile_keeper;
|
||||
} stats1_percentile_state_t;
|
||||
static void stats1_percentile_dingest(void* pvstate, double val) {
|
||||
|
|
|
|||
|
|
@ -22,6 +22,8 @@ MINOR: brew version bump?!?
|
|||
|
||||
MINOR: skewness & kurtosis -> stats1
|
||||
|
||||
MINOR: filter/put option for no auto-convert: all inputs should stay MT_STRING until converted.
|
||||
|
||||
MINOR: faqent re R/mysql/etc inouts
|
||||
MINOR: hold-and-fit regressor doc: 'then put' for residuals; note avoids two-pass & the saving of fit parameters
|
||||
MINOR: faqent re histo w/ min/max is effectively 2-pass (unless you have prior knowledge about the data).
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue