1 /*
2 Gri - A language for scientific graphics programming
3 Copyright (C) 2008 Daniel Kelley
4
5 This program is free software; you can redistribute it and/or modify
6 it under the terms of the GNU General Public License as published by
7 the Free Software Foundation; version 3 of the License, or
8 (at your option) any later version.
9
10 This program is distributed in the hope that it will be useful,
11 but WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
13 GNU General Public License for more details.
14
15 You should have received a copy of the GNU General Public License along
16 with this program; if not, write to the Free Software Foundation, Inc.,
17 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
18 */
19
20 #include <string>
21 #include <algorithm> // for sort
22 #include <vector> // part of STL
23 #include <math.h>
24 #include "gr.hh"
25 #include "private.hh"
26
27
28 extern char _grTempString[];
29
30 static double array_at_i(const std::vector<double>& x,
31 double idouble,
32 int n);
33
34 void
moment(double * data,int n,double * ave,double * adev,double * sdev,double * svar,double * skew,double * kurt)35 moment(double *data,
36 int n,
37 double *ave,
38 double *adev,
39 double *sdev,
40 double *svar,
41 double *skew,
42 double *kurt)
43 {
44 if (n < 2)
45 *ave = *adev = *sdev = *svar = *skew = *kurt = 0.0;
46 else {
47 int ngood = 0, i;
48 double s = 0.0, p;
49 for (i = 0; i < n; i++)
50 if (!gr_missing((double) (*(data + i)))) {
51 s += *(data + i);
52 ngood++;
53 }
54 *ave = s / ngood;
55 *adev = (*svar) = (*skew) = (*kurt) = 0.0;
56 for (i = 0; i < n; i++) {
57 if (!gr_missing((double) (*(data + i)))) {
58 *adev += fabs(s = *(data + i) - (*ave));
59 *svar += (p = s * s);
60 *skew += (p *= s);
61 *kurt += (p *= s);
62 }
63 }
64 *adev /= ngood;
65 *svar /= (ngood - 1);
66 *sdev = sqrt(*svar);
67 if (*svar) {
68 *skew /= (ngood * (*svar) * (*sdev));
69 *kurt = (*kurt) / (ngood * (*svar) * (*svar))/* -3.0*/;
70 }
71 }
72 }
73
74 // calculate q1, q2 = median, and q3 for n data in x
75 void
histogram_stats(const double * x,unsigned int n,double * q1,double * q2,double * q3)76 histogram_stats(const double* x,
77 unsigned int n,
78 double *q1,
79 double *q2,
80 double *q3)
81 {
82 //void sort(vector<double>::iterator, vector<double>::iterator);
83 if (n < 2)
84 *q1 = *q2 = *q3 = 0.0;
85 else {
86 unsigned int ngood = 0;
87 std::vector<double> xcopy;
88 for (unsigned int i = 0; i < n; i++)
89 if (!gr_missing(*(x + i)))
90 xcopy.push_back(*(x + i));
91 ngood = xcopy.size();
92 std::sort(xcopy.begin(), xcopy.end());
93 *q1 = array_at_i(xcopy, 0.25 * (ngood - 1), ngood);
94 *q2 = array_at_i(xcopy, 0.50 * (ngood - 1), ngood);
95 *q3 = array_at_i(xcopy, 0.75 * (ngood - 1), ngood);
96 }
97 PUT_VAR("..q1..", *q1);
98 PUT_VAR("..q2..", *q2);
99 PUT_VAR("..q3..", *q3);
100 }
101
102 static double
array_at_i(const std::vector<double> & x,double idouble,int n)103 array_at_i(const std::vector<double>& x, double idouble, int n)
104 {
105 int i = int(floor(idouble));
106 if (i < 0)
107 return x[0];
108 else if (i >= n)
109 return x[n - 1];
110 else {
111 double r = idouble - i;
112 return (x[i + 1] * r + x[i] * (1.0 - r));
113 }
114 }
115