misc neatens

This commit is contained in:
John Kerl 2015-05-10 12:17:15 -04:00
parent 80093e6efc
commit ae64bf0163
6 changed files with 90 additions and 61 deletions

View file

@ -109,7 +109,6 @@ lrec_t* lrec_parse_dkvp(char* line, char ifs, char ips, int allow_repeat_ifs) {
lrec_put_no_free(prec, key, value);
}
p++;
if (allow_repeat_ifs) {
while (*p == ifs)

View file

@ -20,7 +20,7 @@ static lrec_t* reader_dkvp_func(FILE* input_stream, void* pvstate, context_t* pc
if (line == NULL)
return NULL;
else
return lrec_parse_dkvp(line, pstate->ifs, pstate->ips, FALSE);
return lrec_parse_dkvp(line, pstate->ifs, pstate->ips, pstate->allow_repeat_ifs);
}
// No-op for stateless readers such as this one.

View file

@ -1,6 +1,57 @@
#include <math.h>
#include "lib/mlrstat.h"
// ================================================================
// These are intended for streaming (i.e. single-pass) applications. Otherwise
// the formulas look different (and are more intuitive).
// ================================================================
// ----------------------------------------------------------------
// Univariate linear regression
// ----------------------------------------------------------------
// There are N (xi, yi) pairs.
//
// minimize E = sum (yi - m xi - b)^2
//
// Set the two partial derivatives to zero and solve for m and b:
//
// DE/Dm = sum 2 (yi - m xi - b) (-xi) = 0
// DE/Db = sum 2 (yi - m xi - b) (-1) = 0
//
// sum (yi - m xi - b) (xi) = 0
// sum (yi - m xi - b) = 0
//
// sum (xi yi - m xi^2 - b xi) = 0
// sum (yi - m xi - b) = 0
//
// m sum(xi^2) + b sum(xi) = sum(xi yi)
// m sum(xi) + b N = sum(yi)
//
// [ sum(xi^2) sum(xi) ] [ m ] = [ sum(xi yi) ]
// [ sum(xi) N ] [ b ] = [ sum(yi) ]
//
// [ m ] = [ sum(xi^2) sum(xi) ]^-1 [ sum(xi yi) ]
// [ b ] [ sum(xi) N ] [ sum(yi) ]
//
// = [ N -sum(xi) ] [ sum(xi yi) ] * 1/D
// [ -sum(xi) sum(xi^2)] [ sum(yi) ]
//
// where
//
// D = N sum(xi^2) - sum(xi)^2.
//
// So
//
// N sum(xi yi) - sum(xi) sum(yi)
// m = --------------------------------
// D
//
// -sum(xi)sum(xi yi) + sum(xi^2) sum(yi)
// b = ----------------------------------------
// D
//
// ----------------------------------------------------------------
void mlr_get_linear_regression_ols(unsigned long long n, double sumx, double sumx2, double sumxy, double sumy,
double* pm, double* pb)
{
@ -27,8 +78,7 @@ void mlr_get_linear_regression_ols(unsigned long long n, double sumx, double sum
//
// output = [m, b, math.sqrt(var_m), math.sqrt(var_b)]
// This is intended for streaming (i.e. single-pass) applications. Otherwise
// the formulas look different (and are more intuitive).
// ----------------------------------------------------------------
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;
@ -38,6 +88,7 @@ double mlr_get_stddev(unsigned long long n, double sum, double sum2) {
return sqrt(numerator / denominator);
}
// ----------------------------------------------------------------
double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy) {
double meanx = sumx / n;
double meany = sumy / n;
@ -46,6 +97,7 @@ double mlr_get_cov(unsigned long long n, double sumx, double sumy, double sumxy)
return numerator / denominator;
}
// ----------------------------------------------------------------
void mlr_get_cov_matrix(unsigned long long n,
double sumx, double sumx2, double sumy, double sumy2, double sumxy,
double Q[2][2])
@ -59,17 +111,22 @@ void mlr_get_cov_matrix(unsigned long long n,
// ----------------------------------------------------------------
// 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.
// is zero. The fit equation is 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.
// they'll be about the same. I use 1 - |lambda2|/|lambda1| 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

View file

@ -287,6 +287,9 @@ char* acc_mode_get(void* pvstate, char* pfree_flags) {
// use it on subsequent rows. assumptions:
// * the address doesn't change
// * the content we use (namely, ofmt) isn't row-dependent
// Option 1:
// * modify make_acc to special-case p{n}. needs multi-level hashmap keys
// * do it outside make_acc; requires separate hash maps for percentiles/deciles/quartiles/etc.
acc_t* acc_mode_alloc(static_context_t* pstatx) {
acc_t* pacc = mlr_malloc_or_die(sizeof(acc_t));
acc_mode_state_t* pstate = mlr_malloc_or_die(sizeof(acc_mode_state_t));
@ -317,6 +320,10 @@ static acc_lookup_t acc_lookup_table[] = {
static int acc_lookup_table_length = sizeof(acc_lookup_table) / sizeof(acc_lookup_table[0]);
// xxx make this a hashmap?
// xxx what if acc_name is p50? need:
// * here and here alone is cross-dependence between accumulators
// * if there are min,p10,p50,avg,p90,max then the values array should be
// shared between p10,p50,p90
static acc_t* make_acc(char* acc_name, static_context_t* pstatx) {
for (int i = 0; i < acc_lookup_table_length; i++)
if (streq(acc_name, acc_lookup_table[i].name))
@ -345,11 +352,6 @@ typedef struct _mapper_stats1_state_t {
// ["s","t"] |--> "x" |--> "sum" |--> acc_t* (as void*)
// level_1 level_2 level_3
// lhmslv_t lhmsv_t lhmsv_t
// acc_t implements interface:
// void init();
// void dacc(double dval);
// void sacc(char* sval);
// char* get();
// ----------------------------------------------------------------
sllv_t* mapper_stats1_func(lrec_t* pinrec, context_t* pctx, void* pvstate) {

View file

@ -31,50 +31,6 @@ typedef struct _stats2_t {
typedef stats2_t* stats2_alloc_func_t(static_context_t* pstatx, int do_verbose);
// xxx move to mlrstat.h/c
// ----------------------------------------------------------------
// Univariate linear regression
// ----------------------------------------------------------------
// There are N (xi, yi) pairs.
//
// E = sum (yi - m xi - b)^2
//
// DE/Dm = sum 2 (yi - m xi - b) (-xi) = 0
// DE/Db = sum 2 (yi - m xi - b) (-1) = 0
//
// sum (yi - m xi - b) (xi) = 0
// sum (yi - m xi - b) = 0
//
// sum (xi yi - m xi^2 - b xi) = 0
// sum (yi - m xi - b) = 0
//
// m sum(xi^2) + b sum(xi) = sum(xi yi)
// m sum(xi) + b N = sum(yi)
//
// [ sum(xi^2) sum(xi) ] [ m ] = [ sum(xi yi) ]
// [ sum(xi) N ] [ b ] = [ sum(yi) ]
//
// [ m ] = [ sum(xi^2) sum(xi) ]^-1 [ sum(xi yi) ]
// [ b ] [ sum(xi) N ] [ sum(yi) ]
//
// = [ N -sum(xi) ] [ sum(xi yi) ] * 1/D
// [ -sum(xi) sum(xi^2)] [ sum(yi) ]
//
// where
//
// D = N sum(xi^2) - sum(xi)^2.
//
// So
//
// N sum(xi yi) - sum(xi) sum(yi)
// m = --------------------------------
// D
//
// -sum(xi)sum(xi yi) + sum(xi^2) sum(yi)
// b = ----------------------------------------
// D
typedef struct _stats2_linreg_ols_state_t {
unsigned long long count;
double sumx;

View file

@ -2,27 +2,38 @@
! BUGFIXES !
* --ofmt ignored in put. perhaps best to reglobalize.
rid of ctx.statx; make a mlr_globals_t which is the same.
================================================================
FEATURES
!! quantiles !!
-> be sure to include p99/p50 example (with then-chaining) in mlrwik
!! mode
!! reminder pgr legend is broken !!
http://en.wikipedia.org/wiki/Order_statistic_tree
!! dkvp as generalization of nidx. restructure mlrwik to emphasize this.
tightly integrate 'mlr label'. maybe rename 'mlr label' to 'mlr name' or
some such. perhaps entirely coalesce nidx&dkvp in the code & the docs;
presumably with a different name. something about "header with data" or
"key with value"?? lower-cased only rather than making it an acronym?
!! use 1-(|l2|/|l1|)^2 as pca quality metric? verify against r2 in munch plots.
! sub function. e.g. "300ms" -> "300"
! ordered cut (a la reorder). either a new command (yeck) or cut option (e.g. cut -o)
* stats1 mode: lhmsi & then sort. what about "1"=="1.0"?
* stats1 mode: lhmsi & then sort. what about "1"=="1.0"? doc this, or impl option
w/ temporary sscanf & reformat @ maxlen
* mod op (either c-like, or sane) and put into wikidoc if so.
* linreg-quality 2nd pass -- code it up in stats2 w/ -m {m} -b {b} -- ?
* RV-coefficient -- ?
================================================================
NEATEN
@ -44,9 +55,11 @@ NEATEN
================================================================
ONLINE HELP
* then-chaining note into mlr online help
* jko mlrdoc & gh/jk/mlr urls into mlr online help
* put/filter: have a categorized function lister -- by string/math or arity, or some such ...
* more about I/O and OFMT options @ online help
================================================================
IMPROVEMENTS
@ -109,7 +122,7 @@ DOC
================================================================
PERF
* try mmap(2) for non-stdin case. 1st experiment w/ catc.c.
* try mmap(2) for non-stdin case. 1st experiment w/ catc.c, & make a cutc.c.
================================================================
DATA
@ -123,6 +136,8 @@ MEM MGMT:
* multi-level frees in stats1/stats2/step hashmaps (data-plane structures)
* _free funcptr/funcs for mappers
* free last rec in streamer?
* look strdups at other lhm*
* look at any other strdups
================================================================
FCNS INCL. STRxTIME
@ -157,7 +172,7 @@ INTERNAL DOCS (e.g. README)
================================================================
HARDER HYGIENE
* eliminate compiler warnings for *.l/*.y/etc.
* eliminate compiler warnings for lemon & its autogenerated code
================================================================
PYTHON