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