File: stats.cc

package info (click to toggle)
gri 2.4.2-1
  • links: PTS
  • area: main
  • in suites: potato
  • size: 4,540 kB
  • ctags: 1,966
  • sloc: cpp: 32,542; lisp: 3,243; perl: 806; makefile: 548; sh: 253
file content (80 lines) | stat: -rw-r--r-- 1,999 bytes parent folder | download
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
#include <string>
#include <algorithm>		// for sort
#include <vector>		// part of STL
#include <math.h>
#include "gr.hh"
#include "private.hh"


extern char     _grTempString[];


static double    array_at_i(const double *x, double idouble, int n);

void
moment(double *data, int n, double *ave, double *adev, double *sdev, double *svar, double *skew, double *kurt)
{
    if (n < 2)
	*ave = *adev = *sdev = *svar = *skew = *kurt = 0.0;
    else {
	int             ngood = 0, i;
	double          s = 0.0, p;
	for (i = 0; i < n; i++)
	    if (!gr_missing((double) (*(data + i)))) {
		s += *(data + i);
		ngood++;
	    }
	*ave = s / ngood;
	*adev = (*svar) = (*skew) = (*kurt) = 0.0;
	for (i = 0; i < n; i++) {
	    if (!gr_missing((double) (*(data + i)))) {
		*adev += fabs(s = *(data + i) - (*ave));
		*svar += (p = s * s);
		*skew += (p *= s);
		*kurt += (p *= s);
	    }
	}
	*adev /= ngood;
	*svar /= (ngood - 1);
	*sdev = sqrt(*svar);
	if (*svar) {
	    *skew /= (ngood * (*svar) * (*sdev));
	    *kurt = (*kurt) / (ngood * (*svar) * (*svar)) - 3.0;
	}
    }
}

// calculate q1, q2 = median, and q3 for n data in x
void
histogram_stats(const double *x, unsigned int n, double *q1, double *q2, double *q3)
{
    //void sort(vector<double>::iterator, vector<double>::iterator);
    if (n < 2)
	*q1 = *q2 = *q3 = 0.0;
    else {
	unsigned int ngood = 0;
	vector<double> xcopy;
	for (unsigned int i = 0; i < n; i++)
	    if (!gr_missing(*(x + i)))
		xcopy.push_back(x[i]);
	ngood = xcopy.size();
	sort(xcopy.begin(), xcopy.end());
	*q1 = array_at_i(xcopy.begin(), 0.25 * (ngood - 1), ngood);
	*q2 = array_at_i(xcopy.begin(), 0.50 * (ngood - 1), ngood);
	*q3 = array_at_i(xcopy.begin(), 0.75 * (ngood - 1), ngood);
    }
}

static double
array_at_i(const double *x, double idouble, int n)
{
    int i = int(floor(idouble));
    if (i < 0)
	return *x;
    else if (i >= n)
	return *(x + n - 1);
    else {
	double           r = idouble - i;
	return (*(x + i + 1) * r + *(x + i) * (1.0 - r));
    }
}