1LoadFunctionLibrary("libv3/UtilityFunctions.bf");
2
3math.Infinity = 1/0;
4
5/** @module math */
6
7/**
8* Converts a float to integer by rounding
9* @name math.Int
10* @param float
11*/
12lfunction math.Int (float) {
13  return (float + 0.5)$1;
14}
15
16/**
17* Computes small-sample AIC
18* @name math.GetIC
19* @param logl
20* @param params
21* @param samples
22*/
23lfunction math.GetIC(logl,params,samples) {
24  return -2*logl + 2*samples/(samples-params-1)*params;
25}
26
27/**
28* Computes the LRT and p-value
29* @name math.DoLRT
30* @param lognull
31* @param logalt
32* @param df
33*/
34lfunction math.DoLRT (lognull,logalt,df) {
35  lrt = 2 * (logalt - lognull);
36  return {
37           utility.getGlobalValue("terms.LRT"): lrt__,
38           utility.getGlobalValue("terms.p_value") : 1-CChi2 (lrt__, df__)
39         }
40}
41
42/**
43* Computes sum of 1xN vector
44* @name math.Sum
45* @param {Matrix} _data_vector - 1xN vector
46*/
47lfunction math.Sum(_data_vector) {
48
49  return +_data_vector;
50
51}
52
53/**
54* Returns median average of 1xN vector
55* @name math.Median
56* @param {Matrix} _data_vector - 1xN vector
57*/
58lfunction math.Median(_data_vector) {
59
60  count = utility.Array1D (_data_vector);
61     // this will also let you handle other non 1xN vectors
62
63  // sort tmp_data_vector
64  _tmp_data_vector = _data_vector % 0;
65
66  if (count%2) {
67    median = _tmp_data_vector[count$2];
68    // $ is integer division, so you don't end up taking index 1.5
69  } else {
70    counter = count$2-1;
71    median = (_tmp_data_vector[counter]+_tmp_data_vector[counter+1])/2;
72  }
73
74  return median;
75
76}
77
78/**
79* Returns mean average of 1xN vector
80* @name math.Mean
81* @param {Matrix} _data_vector - 1xN vector
82*/
83lfunction math.Mean(_data_vector) {
84
85  return math.Sum (_data_vector) / utility.Array1D (_data_vector);
86
87}
88
89/**
90* Returns kurtosis of 1xN vector
91* @name math.Kurtosis
92* @param {Matrix} _data_vector - 1xN vector
93*/
94lfunction math.Kurtosis(_data_vector) {
95
96  count = utility.Array1D (_data_vector);
97  mean = math.Mean(_data_vector);
98
99  /*
100  for (_k = 0; _k < count; _k = _k+1) {
101    diff += Abs(_data_vector[_k] - mean)^4;
102  }*/
103
104  // it is MUCH faster to do "iterator ops" like this
105  diff = + _data_vector ["(_MATRIX_ELEMENT_VALUE_ - mean)^4"];
106  _moment = diff/count;
107
108  return _moment/math.Variance(_data_vector)^2;
109
110}
111
112/**
113* Returns skewness of 1xN vector
114* @name math.Skewness
115* @param {Matrix} _data_vector - 1xN vector
116*/
117lfunction math.Skewness(_data_vector) {
118
119  count = utility.Array1D (_data_vector);
120  mean = math.Mean(_data_vector);
121  diff = + _data_vector ["(_MATRIX_ELEMENT_VALUE_ - mean)^3"];
122
123  _moment = diff/count;
124  return _moment/math.Std(_data_vector)^3;
125
126}
127
128/**
129* Returns variance of 1xN vector
130* @name math.Variance
131* @param {Matrix} _data_vector - 1xN vector
132*/
133lfunction math.Variance(_data_vector) {
134
135  count = utility.Array1D (_data_vector);
136  mean = math.Mean(_data_vector);
137
138  diff = + _data_vector ["(_MATRIX_ELEMENT_VALUE_ - mean)^2"];
139  return diff/count;
140
141}
142
143/**
144* Returns standard deviation of 1xN vector
145* @name math.Std
146* @param {Matrix} _data_vector - 1xN vector
147*/
148lfunction math.Std(_data_vector) {
149  return Sqrt (math.Variance (_data_vector));
150};
151
152/**
153* Returns suite of statistical information for a 1xN vector
154* @name math.GatherDescriptiveStats
155* @param {Matrix} _data_vector - 1xN vector
156*/
157lfunction math.GatherDescriptiveStats(_data_vector) {
158
159
160  _count = utility.Array1D (_data_vector);
161  _sorted_data_vector = _data_vector % 0;
162
163  _variance = math.Variance(_data_vector);
164  _std = Sqrt(_variance);
165  _sum = math.Sum(_data_vector);
166
167  _nonnegatives = utility.Filter (_data_vector, "_value_", "_value_ >= 0");
168
169  _dstats = {};
170  _dstats[utility.getGlobalValue("terms.math.count")] = _count;
171  _dstats[utility.getGlobalValue("terms.math.mean")] = math.Mean(_data_vector);
172  _dstats[utility.getGlobalValue("terms.math.median")] = math.Median(_data_vector);
173  _dstats[utility.getGlobalValue("terms.math.minimum")] = _sorted_data_vector[0];
174  _dstats[utility.getGlobalValue("terms.math.maximum")] = _sorted_data_vector[_count - 1];
175  _dstats[utility.getGlobalValue("terms.math._2.5")] = _sorted_data_vector[(_count*0.025)$1];
176  _dstats[utility.getGlobalValue("terms.math._97.5")] = _sorted_data_vector[Min((_count*0.975+0.5)$1, _count-1)];
177  _dstats[utility.getGlobalValue("terms.math.sum")] = _sum;
178  //_dstats["Mean"] = _sum/_count;
179  _dstats[utility.getGlobalValue("terms.math.stddev")] = _std;
180  _dstats[utility.getGlobalValue("terms.math.variance")] = _variance;
181  _dstats[utility.getGlobalValue("terms.math.cov")] = _std*_count/_sum;
182  _dstats[utility.getGlobalValue("terms.math.skewness")] = math.Skewness(_data_vector);
183  _dstats[utility.getGlobalValue("terms.math.kurtosis")] = math.Kurtosis(_data_vector);
184  _dstats [utility.getGlobalValue("terms.math.square_sum")] = math.Sum(_data_vector*2);
185  _dstats [utility.getGlobalValue("terms.math.non_negative")]= Abs(_nonnegatives);
186
187  return _dstats;
188
189}
190
191/**
192* Returns Holm-Bonferroni corrected p-values
193* @name math.HolmBonferroniCorrection
194* @param {Dict} ps - a key -> p-value (or null, for not done)
195
196*/
197lfunction math.HolmBonferroniCorrection(ps) {
198  tested  = utility.Filter (ps, "_value_", "None!=_value_");
199  count   = utility.Array1D (tested);
200  names   = utility.Keys (tested);
201  indexed = {count, 2};
202  for (k = 0; k < count; k+=1) {
203    indexed [k][0] = tested[names[k]];
204    indexed [k][1] = k;
205  }
206  indexed = indexed % 0;
207
208
209  for (k = 0; k < count; k+=1) {
210    indexed[k][0] = indexed[k][0] * (count - k);
211  }
212
213  corrected = {};
214
215  for (k = 0; k < count; k+=1) {
216    corrected[names[indexed[k][1]]] = Min (1,indexed[k][0]);
217  }
218  utility.Extend (corrected, ps);
219  return corrected;
220};
221
222
223/**
224* Returns Benjamini-Hochberg corrected q-values (FDR)
225* @name math.BenjaminiHochbergFDR
226* @param {Dict} ps - a key -> p-value (or null, for not done)
227
228*/
229
230lfunction math.BenjaminiHochbergFDR(ps) {
231  tested  = utility.Filter (ps, "_value_", "None!=_value_");
232  count   = utility.Array1D (tested);
233  names   = utility.Keys (tested);
234  indexed = {count, 2};
235  for (k = 0; k < count; k+=1) {
236    indexed [k][0] = tested[names[k]];
237    indexed [k][1] = k;
238  }
239  indexed = indexed % 0;
240
241
242  for (k = 0; k < count; k+=1) {
243    indexed[k][0] = indexed[k][0] * count/(k+1);
244  }
245
246  corrected = {};
247
248  for (k = 0; k < count; k+=1) {
249    corrected[names[indexed[k][1]]] = Min (1,indexed[k][0]);
250  }
251
252  utility.Extend (corrected, ps);
253  return corrected;
254};
255
256
257/**
258* Returns the range normalized to the lowest value
259* @name math.minNormalizedRange
260* @param {Matrix || Dictonary} if Dictonary, the values will be used
261* @returns {number}
262
263*/
264lfunction math.minNormalizedRange(object) {
265  if (Type (object) == "AssociativeList") {
266    object = utility.Values(object);
267  }
268  min_value = Min(object, 0);
269  if (min_value == 0) {
270    return ^"math.Infinity";
271  }
272  return (Max(object, 0) - min_value ) / min_value;
273};
274
275
276