KOffice – TDE office suite
You can not select more than 25 topics Topics must start with a letter or number, can include dashes ('-') and can be up to 35 characters long.

1244 lines
35KB

  1. /* This file is part of the KDE project
  2. Copyright (C) 1998-2002 The KSpread Team
  3. www.koffice.org/kspread
  4. Copyright (C) 2005 Tomas Mecir <mecirt@gmail.com>
  5. This library is free software; you can redistribute it and/or
  6. modify it under the terms of the GNU Library General Public
  7. License as published by the Free Software Foundation; either
  8. version 2 of the License.
  9. This library is distributed in the hope that it will be useful,
  10. but WITHOUT ANY WARRANTY; without even the implied warranty of
  11. MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
  12. Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public License
  14. along with this library; see the file COPYING.LIB. If not, write to
  15. the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
  16. * Boston, MA 02110-1301, USA.
  17. */
  18. // built-in statistical functions
  19. #include "functions.h"
  20. #include "valuecalc.h"
  21. #include "valueconverter.h"
  22. // needed for MODE
  23. #include <tqmap.h>
  24. using namespace KSpread;
  25. // prototypes (sorted!)
  26. Value func_arrang (valVector args, ValueCalc *calc, FuncExtra *);
  27. Value func_average (valVector args, ValueCalc *calc, FuncExtra *);
  28. Value func_averagea (valVector args, ValueCalc *calc, FuncExtra *);
  29. Value func_avedev (valVector args, ValueCalc *calc, FuncExtra *);
  30. Value func_betadist (valVector args, ValueCalc *calc, FuncExtra *);
  31. Value func_bino (valVector args, ValueCalc *calc, FuncExtra *);
  32. Value func_chidist (valVector args, ValueCalc *calc, FuncExtra *);
  33. Value func_combin (valVector args, ValueCalc *calc, FuncExtra *);
  34. Value func_confidence (valVector args, ValueCalc *calc, FuncExtra *);
  35. Value func_correl_pop (valVector args, ValueCalc *calc, FuncExtra *);
  36. Value func_covar (valVector args, ValueCalc *calc, FuncExtra *);
  37. Value func_devsq (valVector args, ValueCalc *calc, FuncExtra *);
  38. Value func_devsqa (valVector args, ValueCalc *calc, FuncExtra *);
  39. Value func_expondist (valVector args, ValueCalc *calc, FuncExtra *);
  40. Value func_fdist (valVector args, ValueCalc *calc, FuncExtra *);
  41. Value func_fisher (valVector args, ValueCalc *calc, FuncExtra *);
  42. Value func_fisherinv (valVector args, ValueCalc *calc, FuncExtra *);
  43. Value func_gammadist (valVector args, ValueCalc *calc, FuncExtra *);
  44. Value func_gammaln (valVector args, ValueCalc *calc, FuncExtra *);
  45. Value func_gauss (valVector args, ValueCalc *calc, FuncExtra *);
  46. Value func_geomean (valVector args, ValueCalc *calc, FuncExtra *);
  47. Value func_harmean (valVector args, ValueCalc *calc, FuncExtra *);
  48. Value func_hypgeomdist (valVector args, ValueCalc *calc, FuncExtra *);
  49. Value func_kurtosis_est (valVector args, ValueCalc *calc, FuncExtra *);
  50. Value func_kurtosis_pop (valVector args, ValueCalc *calc, FuncExtra *);
  51. Value func_large (valVector args, ValueCalc *calc, FuncExtra *);
  52. Value func_loginv (valVector args, ValueCalc *calc, FuncExtra *);
  53. Value func_lognormdist (valVector args, ValueCalc *calc, FuncExtra *);
  54. Value func_median (valVector args, ValueCalc *calc, FuncExtra *);
  55. Value func_mode (valVector args, ValueCalc *calc, FuncExtra *);
  56. Value func_negbinomdist (valVector args, ValueCalc *calc, FuncExtra *);
  57. Value func_normdist (valVector args, ValueCalc *calc, FuncExtra *);
  58. Value func_norminv (valVector args, ValueCalc *calc, FuncExtra *);
  59. Value func_normsinv (valVector args, ValueCalc *calc, FuncExtra *);
  60. Value func_phi (valVector args, ValueCalc *calc, FuncExtra *);
  61. Value func_poisson (valVector args, ValueCalc *calc, FuncExtra *);
  62. Value func_skew_est (valVector args, ValueCalc *calc, FuncExtra *);
  63. Value func_skew_pop (valVector args, ValueCalc *calc, FuncExtra *);
  64. Value func_small (valVector args, ValueCalc *calc, FuncExtra *);
  65. Value func_standardize (valVector args, ValueCalc *calc, FuncExtra *);
  66. Value func_stddev (valVector args, ValueCalc *calc, FuncExtra *);
  67. Value func_stddeva (valVector args, ValueCalc *calc, FuncExtra *);
  68. Value func_stddevp (valVector args, ValueCalc *calc, FuncExtra *);
  69. Value func_stddevpa (valVector args, ValueCalc *calc, FuncExtra *);
  70. Value func_stdnormdist (valVector args, ValueCalc *calc, FuncExtra *);
  71. Value func_sumproduct (valVector args, ValueCalc *calc, FuncExtra *);
  72. Value func_sumx2py2 (valVector args, ValueCalc *calc, FuncExtra *);
  73. Value func_sumx2my2 (valVector args, ValueCalc *calc, FuncExtra *);
  74. Value func_sumxmy2 (valVector args, ValueCalc *calc, FuncExtra *);
  75. Value func_tdist (valVector args, ValueCalc *calc, FuncExtra *);
  76. Value func_variance (valVector args, ValueCalc *calc, FuncExtra *);
  77. Value func_variancea (valVector args, ValueCalc *calc, FuncExtra *);
  78. Value func_variancep (valVector args, ValueCalc *calc, FuncExtra *);
  79. Value func_variancepa (valVector args, ValueCalc *calc, FuncExtra *);
  80. Value func_weibull (valVector args, ValueCalc *calc, FuncExtra *);
  81. typedef TQValueList<double> List;
  82. // registers all statistical functions
  83. void RegisterStatisticalFunctions()
  84. {
  85. FunctionRepository* repo = FunctionRepository::self();
  86. Function *f;
  87. f = new Function ("AVEDEV", func_avedev);
  88. f->setParamCount (1, -1);
  89. f->setAcceptArray ();
  90. repo->add (f);
  91. f = new Function ("AVERAGE", func_average);
  92. f->setParamCount (1, -1);
  93. f->setAcceptArray ();
  94. repo->add (f);
  95. f = new Function ("AVERAGEA", func_averagea);
  96. f->setParamCount (1, -1);
  97. f->setAcceptArray ();
  98. repo->add (f);
  99. f = new Function ("BETADIST", func_betadist);
  100. f->setParamCount (3, 5);
  101. repo->add (f);
  102. f = new Function ("BINO", func_bino);
  103. f->setParamCount (3);
  104. repo->add (f);
  105. f = new Function ("CHIDIST", func_chidist);
  106. f->setParamCount (2);
  107. repo->add (f);
  108. f = new Function ("COMBIN", func_combin);
  109. f->setParamCount (2);
  110. repo->add (f);
  111. f = new Function ("CONFIDENCE", func_confidence);
  112. f->setParamCount (3);
  113. repo->add (f);
  114. f = new Function ("CORREL", func_correl_pop);
  115. f->setParamCount (2);
  116. f->setAcceptArray ();
  117. repo->add (f);
  118. f = new Function ("COVAR", func_covar);
  119. f->setParamCount (2);
  120. f->setAcceptArray ();
  121. repo->add (f);
  122. f = new Function ("DEVSQ", func_devsq);
  123. f->setParamCount (1, -1);
  124. f->setAcceptArray ();
  125. repo->add (f);
  126. f = new Function ("DEVSQA", func_devsqa);
  127. f->setParamCount (1, -1);
  128. f->setAcceptArray ();
  129. repo->add (f);
  130. f = new Function ("EXPONDIST", func_expondist);
  131. f->setParamCount (3);
  132. repo->add (f);
  133. f = new Function ("FDIST", func_fdist);
  134. f->setParamCount (3);
  135. repo->add (f);
  136. f = new Function ("FISHER", func_fisher);
  137. repo->add (f);
  138. f = new Function ("FISHERINV", func_fisherinv);
  139. repo->add (f);
  140. f = new Function ("GAMMADIST", func_gammadist);
  141. f->setParamCount (4);
  142. repo->add (f);
  143. f = new Function ("GAMMALN", func_gammaln);
  144. repo->add (f);
  145. f = new Function ("GAUSS", func_gauss);
  146. repo->add (f);
  147. f = new Function ("GEOMEAN", func_geomean);
  148. f->setParamCount (1, -1);
  149. f->setAcceptArray ();
  150. repo->add (f);
  151. f = new Function ("HARMEAN", func_harmean);
  152. f->setParamCount (1, -1);
  153. f->setAcceptArray ();
  154. repo->add (f);
  155. f = new Function ("HYPGEOMDIST", func_hypgeomdist);
  156. f->setParamCount (4);
  157. repo->add (f);
  158. f = new Function ("INVBINO", func_bino); // same as BINO, for 1.4 compat
  159. repo->add (f);
  160. f = new Function ("KURT", func_kurtosis_est);
  161. f->setParamCount (1, -1);
  162. f->setAcceptArray ();
  163. repo->add (f);
  164. f = new Function ("KURTP", func_kurtosis_pop);
  165. f->setParamCount (1, -1);
  166. f->setAcceptArray ();
  167. repo->add (f);
  168. f = new Function ("LARGE", func_large);
  169. f->setParamCount (2);
  170. f->setAcceptArray ();
  171. repo->add (f);
  172. f = new Function ("LOGINV", func_loginv);
  173. f->setParamCount (3);
  174. repo->add (f);
  175. f = new Function ("LOGNORMDIST", func_lognormdist);
  176. f->setParamCount (3);
  177. repo->add (f);
  178. f = new Function ("MEDIAN", func_median);
  179. f->setParamCount (1, -1);
  180. f->setAcceptArray ();
  181. repo->add (f);
  182. f = new Function ("MODE", func_mode);
  183. f->setParamCount (1, -1);
  184. f->setAcceptArray ();
  185. repo->add (f);
  186. f = new Function ("NEGBINOMDIST", func_negbinomdist);
  187. f->setParamCount (3);
  188. repo->add (f);
  189. f = new Function ("NORMDIST", func_normdist);
  190. f->setParamCount (4);
  191. repo->add (f);
  192. f = new Function ("NORMINV", func_norminv);
  193. f->setParamCount (3);
  194. repo->add (f);
  195. f = new Function ("NORMSDIST", func_stdnormdist);
  196. repo->add (f);
  197. f = new Function ("NORMSINV", func_normsinv);
  198. repo->add (f);
  199. f = new Function ("PEARSON", func_correl_pop);
  200. f->setParamCount (2);
  201. f->setAcceptArray ();
  202. repo->add (f);
  203. f = new Function ("PERMUT", func_arrang);
  204. f->setParamCount (2);
  205. repo->add (f);
  206. f = new Function ("PHI", func_phi);
  207. repo->add (f);
  208. f = new Function ("POISSON", func_poisson);
  209. f->setParamCount (3);
  210. repo->add (f);
  211. f = new Function ("SKEW", func_skew_est);
  212. f->setParamCount (1, -1);
  213. f->setAcceptArray ();
  214. repo->add (f);
  215. f = new Function ("SKEWP", func_skew_pop);
  216. f->setParamCount (1, -1);
  217. f->setAcceptArray ();
  218. repo->add (f);
  219. f = new Function ("SMALL", func_small);
  220. f->setParamCount (2);
  221. f->setAcceptArray ();
  222. repo->add (f);
  223. f = new Function ("STANDARDIZE", func_standardize);
  224. f->setParamCount (3);
  225. repo->add (f);
  226. f = new Function ("STDEV", func_stddev);
  227. f->setParamCount (1, -1);
  228. f->setAcceptArray ();
  229. repo->add (f);
  230. f = new Function ("STDEVA", func_stddeva);
  231. f->setParamCount (1, -1);
  232. f->setAcceptArray ();
  233. repo->add (f);
  234. f = new Function ("STDEVP", func_stddevp);
  235. f->setParamCount (1, -1);
  236. f->setAcceptArray ();
  237. repo->add (f);
  238. f = new Function ("STDEVPA", func_stddevpa);
  239. f->setParamCount (1, -1);
  240. f->setAcceptArray ();
  241. repo->add (f);
  242. f = new Function ("SUM2XMY", func_sumxmy2);
  243. f->setParamCount (2);
  244. f->setAcceptArray ();
  245. repo->add (f);
  246. f = new Function ("SUMPRODUCT", func_sumproduct);
  247. f->setParamCount (2);
  248. f->setAcceptArray ();
  249. repo->add (f);
  250. f = new Function ("SUMX2PY2", func_sumx2py2);
  251. f->setParamCount (2);
  252. f->setAcceptArray ();
  253. repo->add (f);
  254. f = new Function ("SUMX2MY2", func_sumx2my2);
  255. f->setParamCount (2);
  256. f->setAcceptArray ();
  257. repo->add (f);
  258. f = new Function ("TDIST", func_tdist);
  259. f->setParamCount (3);
  260. repo->add (f);
  261. f = new Function ("VARIANCE", func_variance);
  262. f->setParamCount (1, -1);
  263. f->setAcceptArray ();
  264. repo->add (f);
  265. f = new Function ("VAR", func_variance);
  266. f->setParamCount (1, -1);
  267. f->setAcceptArray ();
  268. repo->add (f);
  269. f = new Function ("VARP", func_variancep);
  270. f->setParamCount (1, -1);
  271. f->setAcceptArray ();
  272. repo->add (f);
  273. f = new Function ("VARA", func_variancea);
  274. f->setParamCount (1, -1);
  275. f->setAcceptArray ();
  276. repo->add (f);
  277. f = new Function ("VARPA", func_variancepa);
  278. f->setParamCount (1, -1);
  279. f->setAcceptArray ();
  280. repo->add (f);
  281. f = new Function ("WEIBULL", func_weibull);
  282. f->setParamCount (4);
  283. repo->add (f);
  284. }
  285. // array-walk functions used in this file
  286. void awSkew (ValueCalc *c, Value &res, Value val, Value p)
  287. {
  288. Value avg = p.element (0, 0);
  289. Value stdev = p.element (1, 0);
  290. // (val - avg) / stddev
  291. Value d = c->div (c->sub (val, avg), stdev);
  292. // res += d*d*d
  293. res = c->add (res, c->mul (d, c->mul (d, d)));
  294. }
  295. void awSumInv (ValueCalc *c, Value &res, Value val, Value)
  296. {
  297. // res += 1/value
  298. res = c->add (res, c->div (1.0, val));
  299. }
  300. void awAveDev (ValueCalc *c, Value &res, Value val,
  301. Value p)
  302. {
  303. // res += abs (val - p)
  304. res = c->add (res, c->abs (c->sub (val, p)));
  305. }
  306. void awKurtosis (ValueCalc *c, Value &res, Value val,
  307. Value p)
  308. {
  309. Value avg = p.element (0, 0);
  310. Value stdev = p.element (1, 0);
  311. //d = (val - avg ) / stdev
  312. Value d = c->div (c->sub (val, avg), stdev);
  313. // res += d^4
  314. res = c->add (res, c->pow (d, 4));
  315. }
  316. Value func_skew_est (valVector args, ValueCalc *calc, FuncExtra *)
  317. {
  318. int number = calc->count (args);
  319. Value avg = calc->avg (args);
  320. if (number < 3)
  321. return Value::errorVALUE();
  322. Value res = calc->stddev (args, avg);
  323. if (res.isZero())
  324. return Value::errorVALUE();
  325. Value params (2, 1);
  326. params.setElement (0, 0, avg);
  327. params.setElement (1, 0, res);
  328. Value tskew;
  329. calc->arrayWalk (args, tskew, awSkew, params);
  330. // ((tskew * number) / (number-1)) / (number-2)
  331. return calc->div (calc->div (calc->mul (tskew, number), number-1), number-2);
  332. }
  333. Value func_skew_pop (valVector args, ValueCalc *calc, FuncExtra *)
  334. {
  335. int number = calc->count (args);
  336. Value avg = calc->avg (args);
  337. if (number < 1)
  338. return Value::errorVALUE();
  339. Value res = calc->stddevP (args, avg);
  340. if (res.isZero())
  341. return Value::errorVALUE();
  342. Value params (2, 1);
  343. params.setElement (0, 0, avg);
  344. params.setElement (1, 0, res);
  345. Value tskew;
  346. calc->arrayWalk (args, tskew, awSkew, params);
  347. // tskew / number
  348. return calc->div (tskew, number);
  349. }
  350. class ContentSheet : public TQMap<double, int> {};
  351. void func_mode_helper (Value range, ValueCalc *calc, ContentSheet &sh)
  352. {
  353. if (!range.isArray())
  354. {
  355. double d = calc->conv()->asFloat (range).asFloat();
  356. sh[d]++;
  357. return;
  358. }
  359. for (unsigned int row = 0; row < range.rows(); ++row)
  360. for (unsigned int col = 0; col < range.columns(); ++col) {
  361. Value v = range.element (col, row);
  362. if (v.isArray())
  363. func_mode_helper (v, calc, sh);
  364. else {
  365. double d = calc->conv()->asFloat (v).asFloat();
  366. sh[d]++;
  367. }
  368. }
  369. }
  370. Value func_mode (valVector args, ValueCalc *calc, FuncExtra *)
  371. {
  372. // does NOT support anything other than doubles !!!
  373. ContentSheet sh;
  374. for (unsigned int i = 0; i < args.count(); ++i)
  375. func_mode_helper (args[i], calc, sh);
  376. // retrieve value with max.count
  377. int maxcount = 0;
  378. double max = 0.0;
  379. ContentSheet::iterator it;
  380. for (it = sh.begin(); it != sh.end(); ++it)
  381. if (it.data() > maxcount) {
  382. max = it.key();
  383. maxcount = it.data();
  384. }
  385. return Value (max);
  386. }
  387. Value func_covar_helper (Value range1, Value range2,
  388. ValueCalc *calc, Value avg1, Value avg2)
  389. {
  390. // two arrays -> cannot use arrayWalk
  391. if ((!range1.isArray()) && (!range2.isArray()))
  392. // (v1-E1)*(v2-E2)
  393. return calc->mul (calc->sub (range1, avg1), calc->sub (range2, avg2));
  394. int rows = range1.rows();
  395. int cols = range1.columns();
  396. int rows2 = range2.rows();
  397. int cols2 = range2.columns();
  398. if ((rows != rows2) || (cols != cols2))
  399. return Value::errorVALUE();
  400. Value result = 0.0;
  401. for (int row = 0; row < rows; ++row)
  402. for (int col = 0; col < cols; ++col) {
  403. Value v1 = range1.element (col, row);
  404. Value v2 = range2.element (col, row);
  405. if (v1.isArray() || v2.isArray())
  406. result = calc->add (result,
  407. func_covar_helper (v1, v2, calc, avg1, avg2));
  408. else
  409. // result += (v1-E1)*(v2-E2)
  410. result = calc->add (result, calc->mul (calc->sub (v1, avg1),
  411. calc->sub (v2, avg2)));
  412. }
  413. return result;
  414. }
  415. Value func_covar (valVector args, ValueCalc *calc, FuncExtra *)
  416. {
  417. Value avg1 = calc->avg (args[0]);
  418. Value avg2 = calc->avg (args[1]);
  419. int number = calc->count (args[0]);
  420. int number2 = calc->count (args[1]);
  421. if (number2 <= 0 || number2 != number)
  422. return Value::errorVALUE();
  423. Value covar = func_covar_helper (args[0], args[1], calc, avg1, avg2);
  424. return calc->div (covar, number);
  425. }
  426. Value func_correl_pop (valVector args, ValueCalc *calc, FuncExtra *)
  427. {
  428. Value covar = func_covar (args, calc, 0);
  429. Value stdevp1 = calc->stddevP (args[0]);
  430. Value stdevp2 = calc->stddevP (args[1]);
  431. if (calc->isZero (stdevp1) || calc->isZero (stdevp2))
  432. return Value::errorDIV0();
  433. // covar / (stdevp1 * stdevp2)
  434. return calc->div (covar, calc->mul (stdevp1, stdevp2));
  435. }
  436. void func_array_helper (Value range, ValueCalc *calc,
  437. List &array, int &number)
  438. {
  439. if (!range.isArray())
  440. {
  441. array << calc->conv()->asFloat (range).asFloat();
  442. ++number;
  443. return;
  444. }
  445. for (unsigned int row = 0; row < range.rows(); ++row)
  446. for (unsigned int col = 0; col < range.columns(); ++col) {
  447. Value v = range.element (col, row);
  448. if (v.isArray ())
  449. func_array_helper (v, calc, array, number);
  450. else {
  451. array << calc->conv()->asFloat (v).asFloat();
  452. ++number;
  453. }
  454. }
  455. }
  456. Value func_large (valVector args, ValueCalc *calc, FuncExtra *)
  457. {
  458. // does NOT support anything other than doubles !!!
  459. int k = calc->conv()->asInteger (args[1]).asInteger();
  460. if ( k < 1 )
  461. return false;
  462. List array;
  463. int number = 1;
  464. func_array_helper (args[0], calc, array, number);
  465. if ( k > number )
  466. return Value::errorVALUE();
  467. qHeapSort (array);
  468. double d = *array.at (number - k - 1);
  469. return Value (d);
  470. }
  471. Value func_small (valVector args, ValueCalc *calc, FuncExtra *)
  472. {
  473. // does NOT support anything other than doubles !!!
  474. int k = calc->conv()->asInteger (args[1]).asInteger();
  475. if ( k < 1 )
  476. return false;
  477. List array;
  478. int number = 1;
  479. func_array_helper (args[0], calc, array, number);
  480. if ( k > number )
  481. return Value::errorVALUE();
  482. qHeapSort (array);
  483. double d = *array.at (k - 1);
  484. return Value (d);
  485. }
  486. Value func_geomean (valVector args, ValueCalc *calc, FuncExtra *)
  487. {
  488. Value count = calc->count (args);
  489. Value prod = calc->product (args, 1.0);
  490. if (calc->isZero (count))
  491. return Value::errorDIV0();
  492. return calc->pow (prod, calc->div (1.0, count));
  493. }
  494. Value func_harmean (valVector args, ValueCalc *calc, FuncExtra *)
  495. {
  496. Value count = calc->count (args);
  497. if (calc->isZero (count))
  498. return Value::errorDIV0();
  499. Value suminv;
  500. calc->arrayWalk (args, suminv, awSumInv, 0);
  501. return calc->div (suminv, count);
  502. }
  503. Value func_loginv (valVector args, ValueCalc *calc, FuncExtra *)
  504. {
  505. Value p = args[0];
  506. Value m = args[1];
  507. Value s = args[2];
  508. if (calc->lower (p, 0) || calc->greater (p, 1))
  509. return Value::errorVALUE();
  510. if (!calc->greater (s, 0))
  511. return Value::errorVALUE();
  512. Value result = 0.0;
  513. if (calc->equal (p, 1)) //p==1
  514. result = Value::errorVALUE();
  515. else if (calc->greater (p, 0)) { //p>0
  516. Value gaussInv = calc->gaussinv (p);
  517. // exp (gaussInv * s + m)
  518. result = calc->exp (calc->add (calc->mul (s, gaussInv), m));
  519. }
  520. return result;
  521. }
  522. Value func_devsq (valVector args, ValueCalc *calc, FuncExtra *)
  523. {
  524. Value res;
  525. calc->arrayWalk (args, res, calc->awFunc ("devsq"), calc->avg (args, false));
  526. return res;
  527. }
  528. Value func_devsqa (valVector args, ValueCalc *calc, FuncExtra *)
  529. {
  530. Value res;
  531. calc->arrayWalk (args, res, calc->awFunc ("devsqa"), calc->avg (args));
  532. return res;
  533. }
  534. Value func_kurtosis_est (valVector args, ValueCalc *calc, FuncExtra *)
  535. {
  536. int count = calc->count (args);
  537. if (count < 4)
  538. return Value::errorVALUE();
  539. Value avg = calc->avg (args);
  540. Value devsq;
  541. calc->arrayWalk (args, devsq, calc->awFunc ("devsqa"), avg);
  542. if (devsq.isZero ())
  543. return Value::errorDIV0();
  544. Value params (2, 1);
  545. params.setElement (0, 0, avg);
  546. params.setElement (1, 0, devsq);
  547. Value x4;
  548. calc->arrayWalk (args, x4, awKurtosis, params);
  549. double den = (double) (count - 2) * (count - 3);
  550. double nth = (double) count * (count + 1) / ((count - 1) * den);
  551. double t = 3.0 * (count - 1) * (count - 1) / den;
  552. // res = x4 * nth - t
  553. return calc->sub (calc->mul (x4, nth), t);
  554. }
  555. Value func_kurtosis_pop (valVector args, ValueCalc *calc, FuncExtra *)
  556. {
  557. int count = calc->count (args);
  558. if (count < 4)
  559. return Value::errorVALUE();
  560. Value avg = calc->avg (args);
  561. Value devsq;
  562. calc->arrayWalk (args, devsq, calc->awFunc ("devsqa"), avg);
  563. if (devsq.isZero ())
  564. return Value::errorDIV0();
  565. Value params (2, 1);
  566. params.setElement (0, 0, avg);
  567. params.setElement (1, 0, devsq);
  568. Value x4;
  569. calc->arrayWalk (args, x4, awKurtosis, params);
  570. // x4 / count - 3
  571. return calc->sub (calc->div (x4, count), 3);
  572. }
  573. Value func_standardize (valVector args, ValueCalc *calc, FuncExtra *)
  574. {
  575. Value x = args[0];
  576. Value m = args[1];
  577. Value s = args[2];
  578. if (!calc->greater (s, 0)) // s must be >0
  579. return Value::errorVALUE();
  580. // (x - m) / s
  581. return calc->div (calc->sub (x, m), s);
  582. }
  583. Value func_hypgeomdist (valVector args, ValueCalc *calc, FuncExtra *)
  584. {
  585. int x = calc->conv()->asInteger (args[0]).asInteger();
  586. int n = calc->conv()->asInteger (args[1]).asInteger();
  587. int M = calc->conv()->asInteger (args[2]).asInteger();
  588. int N = calc->conv()->asInteger (args[3]).asInteger();
  589. if ( x < 0 || n < 0 || M < 0 || N < 0 )
  590. return Value::errorVALUE();
  591. if ( x > M || n > N )
  592. return Value::errorVALUE();
  593. Value d1 = calc->combin (M, x);
  594. Value d2 = calc->combin (N - M, n - x);
  595. Value d3 = calc->combin (N, n);
  596. // d1 * d2 / d3
  597. return calc->div (calc->mul (d1, d2), d3);
  598. }
  599. Value func_negbinomdist (valVector args, ValueCalc *calc, FuncExtra *)
  600. {
  601. int x = calc->conv()->asInteger (args[0]).asInteger();
  602. int r = calc->conv()->asInteger (args[1]).asInteger();
  603. Value p = args[2];
  604. if ((x + r - 1) <= 0)
  605. return Value::errorVALUE();
  606. if (calc->lower (p, 0) || calc->greater (p, 1))
  607. return Value::errorVALUE();
  608. Value d1 = calc->combin (x + r - 1, r - 1);
  609. // d2 = pow (p, r) * pow (1 - p, x)
  610. Value d2 = calc->mul (calc->pow (p, r),
  611. calc->pow (calc->sub (1, p), x));
  612. return calc->mul (d1, d2);
  613. }
  614. // Function: permut
  615. Value func_arrang (valVector args, ValueCalc *calc, FuncExtra *)
  616. {
  617. Value n = args[0];
  618. Value m = args[1];
  619. if (calc->lower (n, m)) // problem if n<m
  620. return Value::errorVALUE();
  621. if (calc->lower (m, 0)) // problem if m<0 (n>=m so that's okay)
  622. return Value::errorVALUE();
  623. // fact(n) / (fact(n-m)
  624. return calc->fact (n, calc->sub (n, m));
  625. }
  626. // Function: average
  627. Value func_average (valVector args, ValueCalc *calc, FuncExtra *)
  628. {
  629. return calc->avg (args, false);
  630. }
  631. // Function: averagea
  632. Value func_averagea (valVector args, ValueCalc *calc, FuncExtra *)
  633. {
  634. return calc->avg (args);
  635. }
  636. // Function: avedev
  637. Value func_avedev (valVector args, ValueCalc *calc, FuncExtra *)
  638. {
  639. Value result;
  640. calc->arrayWalk (args, result, awAveDev, calc->avg (args));
  641. return result;
  642. }
  643. // Function: median
  644. Value func_median (valVector args, ValueCalc *calc, FuncExtra *)
  645. {
  646. // does NOT support anything other than doubles !!!
  647. List array;
  648. int number = 1;
  649. for (unsigned int i = 0; i < args.count(); ++i)
  650. func_array_helper (args[i], calc, array, number);
  651. qHeapSort (array);
  652. double d = *array.at (number / 2 + number % 2);
  653. return Value (d);
  654. }
  655. // Function: variance
  656. Value func_variance (valVector args, ValueCalc *calc, FuncExtra *)
  657. {
  658. int count = calc->count (args, false);
  659. if (count < 2)
  660. return Value::errorVALUE();
  661. Value result = func_devsq (args, calc, 0);
  662. return calc->div (result, count-1);
  663. }
  664. // Function: vara
  665. Value func_variancea (valVector args, ValueCalc *calc, FuncExtra *)
  666. {
  667. int count = calc->count (args);
  668. if (count < 2)
  669. return Value::errorVALUE();
  670. Value result = func_devsqa (args, calc, 0);
  671. return calc->div (result, count-1);
  672. }
  673. // Function: varp
  674. Value func_variancep (valVector args, ValueCalc *calc, FuncExtra *)
  675. {
  676. int count = calc->count (args, false);
  677. if (count == 0)
  678. return Value::errorVALUE();
  679. Value result = func_devsq (args, calc, 0);
  680. return calc->div (result, count);
  681. }
  682. // Function: varpa
  683. Value func_variancepa (valVector args, ValueCalc *calc, FuncExtra *)
  684. {
  685. int count = calc->count (args);
  686. if (count == 0)
  687. return Value::errorVALUE();
  688. Value result = func_devsqa (args, calc, 0);
  689. return calc->div (result, count);
  690. }
  691. // Function: stddev
  692. Value func_stddev (valVector args, ValueCalc *calc, FuncExtra *)
  693. {
  694. return calc->stddev (args, false);
  695. }
  696. // Function: stddeva
  697. Value func_stddeva (valVector args, ValueCalc *calc, FuncExtra *)
  698. {
  699. return calc->stddev (args);
  700. }
  701. // Function: stddevp
  702. Value func_stddevp (valVector args, ValueCalc *calc, FuncExtra *)
  703. {
  704. return calc->stddevP (args, false);
  705. }
  706. // Function: stddevpa
  707. Value func_stddevpa (valVector args, ValueCalc *calc, FuncExtra *)
  708. {
  709. return calc->stddevP (args);
  710. }
  711. // Function: combin
  712. Value func_combin (valVector args, ValueCalc *calc, FuncExtra *)
  713. {
  714. return calc->combin (args[0], args[1]);
  715. }
  716. // Function: bino
  717. Value func_bino (valVector args, ValueCalc *calc, FuncExtra *)
  718. {
  719. Value n = args[0];
  720. Value m = args[1];
  721. Value comb = calc->combin (n, m);
  722. Value prob = args[2];
  723. if (calc->lower (prob,0) || calc->greater (prob,1))
  724. return Value::errorVALUE();
  725. // result = comb * pow (prob, m) * pow (1 - prob, n - m)
  726. Value pow1 = calc->pow (prob, m);
  727. Value pow2 = calc->pow (calc->sub (1, prob), calc->sub (n, m));
  728. return calc->mul (comb, calc->mul (pow1, pow2));
  729. }
  730. // Function: phi
  731. Value func_phi (valVector args, ValueCalc *calc, FuncExtra *)
  732. //distribution function for a standard normal distribution
  733. {
  734. return calc->phi (args[0]);
  735. }
  736. // Function: gauss
  737. Value func_gauss (valVector args, ValueCalc *calc, FuncExtra *)
  738. {
  739. //returns the integral values of the standard normal cumulative distribution
  740. return calc->gauss (args[0]);
  741. }
  742. // Function: gammadist
  743. Value func_gammadist (valVector args, ValueCalc *calc, FuncExtra *)
  744. {
  745. Value x = args[0];
  746. Value alpha = args[1];
  747. Value beta = args[2];
  748. int kum = calc->conv()->asInteger (args[3]).asInteger(); // 0 or 1
  749. Value result;
  750. if (calc->lower (x, 0.0) || (!calc->greater (alpha, 0.0)) ||
  751. (!calc->greater (beta, 0.0)))
  752. return Value::errorVALUE();
  753. if (kum == 0) { //density
  754. Value G = calc->GetGamma (alpha);
  755. // result = pow (x, alpha - 1.0) / exp (x / beta) / pow (beta, alpha) / G
  756. Value pow1 = calc->pow (x, calc->sub (alpha, 1.0));
  757. Value pow2 = calc->exp (calc->div (x, beta));
  758. Value pow3 = calc->pow (beta, alpha);
  759. result = calc->div (calc->div (calc->div (pow1, pow2), pow3), G);
  760. }
  761. else
  762. result = calc->GetGammaDist (x, alpha, beta);
  763. return Value (result);
  764. }
  765. // Function: betadist
  766. Value func_betadist (valVector args, ValueCalc *calc, FuncExtra *)
  767. {
  768. Value x = args[0];
  769. Value alpha = args[1];
  770. Value beta = args[2];
  771. Value fA = 0.0;
  772. Value fB = 1.0;
  773. if (args.count() > 3) fA = args[3];
  774. if (args.count() == 5) fB = args[4];
  775. //x < fA || x > fB || fA == fB || alpha <= 0.0 || beta <= 0.0
  776. if (calc->lower (x, fA) || calc->greater (x, fB) || calc->equal (fA, fB) ||
  777. (!calc->greater (alpha, 0.0)) || (!calc->greater (beta, 0.0)))
  778. return Value::errorVALUE();
  779. // xx = (x - fA) / (fB - fA) // scaling
  780. Value xx = calc->div (calc->sub (x, fA), calc->sub (fB, fA));
  781. return calc->GetBeta (xx, alpha, beta);
  782. }
  783. // Function: fisher
  784. Value func_fisher (valVector args, ValueCalc *calc, FuncExtra *) {
  785. //returns the Fisher transformation for x
  786. // 0.5 * ln ((1.0 + fVal) / (1.0 - fVal))
  787. Value fVal = args[0];
  788. Value num = calc->div (calc->add (fVal, 1.0), calc->sub (1.0, fVal));
  789. return calc->mul (calc->ln (num), 0.5);
  790. }
  791. // Function: fisherinv
  792. Value func_fisherinv (valVector args, ValueCalc *calc, FuncExtra *) {
  793. //returns the inverse of the Fisher transformation for x
  794. Value fVal = args[0];
  795. // (exp (2.0 * fVal) - 1.0) / (exp (2.0 * fVal) + 1.0)
  796. Value ex = calc->exp (calc->mul (fVal, 2.0));
  797. return calc->div (calc->sub (ex, 1.0), calc->add (ex, 1.0));
  798. }
  799. // Function: normdist
  800. Value func_normdist (valVector args, ValueCalc *calc, FuncExtra *) {
  801. //returns the normal cumulative distribution
  802. Value x = args[0];
  803. Value mue = args[1];
  804. Value sigma = args[2];
  805. Value k = args[3];
  806. if (!calc->greater (sigma, 0.0))
  807. return Value::errorVALUE();
  808. // (x - mue) / sigma
  809. Value Y = calc->div (calc->sub (x, mue), sigma);
  810. if (calc->isZero (k)) // density
  811. return calc->div (calc->phi (Y), sigma);
  812. else // distribution
  813. return calc->add (calc->gauss (Y), 0.5);
  814. }
  815. // Function: lognormdist
  816. Value func_lognormdist (valVector args, ValueCalc *calc, FuncExtra *) {
  817. //returns the cumulative lognormal distribution
  818. Value x = args[0];
  819. Value mue = args[1];
  820. Value sigma = args[2];
  821. if (!calc->greater (sigma, 0.0) || (!calc->greater (x, 0.0)))
  822. return Value::errorVALUE();
  823. // (ln(x) - mue) / sigma
  824. Value Y = calc->div (calc->sub (calc->ln (x), mue), sigma);
  825. return calc->add (calc->gauss (Y), 0.5);
  826. }
  827. // Function: normsdist
  828. Value func_stdnormdist (valVector args, ValueCalc *calc, FuncExtra *)
  829. {
  830. //returns the cumulative lognormal distribution, mue=0, sigma=1
  831. return calc->add (calc->gauss (args[0]), 0.5);
  832. }
  833. // Function: expondist
  834. Value func_expondist (valVector args, ValueCalc *calc, FuncExtra *) {
  835. //returns the exponential distribution
  836. Value x = args[0];
  837. Value lambda = args[1];
  838. Value kum = args[2];
  839. Value result = 0.0;
  840. if (!calc->greater (lambda, 0.0))
  841. return Value::errorVALUE();
  842. // ex = exp (-lambda * x)
  843. Value ex = calc->exp (calc->mul (calc->mul (lambda, -1), x));
  844. if (calc->isZero (kum)) { //density
  845. if (!calc->lower (x, 0.0))
  846. // lambda * ex
  847. result = calc->mul (lambda, ex);
  848. }
  849. else { //distribution
  850. if (calc->greater (x, 0.0))
  851. // 1.0 - ex
  852. result = calc->sub (1.0, ex);
  853. }
  854. return result;
  855. }
  856. // Function: weibull
  857. Value func_weibull (valVector args, ValueCalc *calc, FuncExtra *) {
  858. //returns the Weibull distribution
  859. Value x = args[0];
  860. Value alpha = args[1];
  861. Value beta = args[2];
  862. Value kum = args[3];
  863. Value result;
  864. if ((!calc->greater (alpha, 0.0)) || (!calc->greater (beta, 0.0)) ||
  865. calc->lower (x, 0.0))
  866. return Value::errorVALUE();
  867. // ex = exp (-pow (x / beta, alpha))
  868. Value ex;
  869. ex = calc->exp (calc->mul (calc->pow (calc->div (x, beta), alpha), -1));
  870. if (calc->isZero (kum)) // density
  871. {
  872. // result = alpha / pow(beta,alpha) * pow(x,alpha-1.0) * ex
  873. result = calc->div (alpha, calc->pow (beta, alpha));
  874. result = calc->mul (result, calc->mul (calc->pow (x,
  875. calc->sub (alpha, 1)), ex));
  876. }
  877. else // distribution
  878. result = calc->sub (1.0, ex);
  879. return result;
  880. }
  881. // Function: normsinv
  882. Value func_normsinv (valVector args, ValueCalc *calc, FuncExtra *) {
  883. //returns the inverse of the standard normal cumulative distribution
  884. Value x = args[0];
  885. if (!(calc->greater (x, 0.0) && calc->lower (x, 1.0)))
  886. return Value::errorVALUE();
  887. return calc->gaussinv (x);
  888. }
  889. // Function: norminv
  890. Value func_norminv (valVector args, ValueCalc *calc, FuncExtra *) {
  891. //returns the inverse of the normal cumulative distribution
  892. Value x = args[0];
  893. Value mue = args[1];
  894. Value sigma = args[2];
  895. if (!calc->greater (sigma, 0.0))
  896. return Value::errorVALUE();
  897. if (!(calc->greater (x, 0.0) && calc->lower (x, 1.0)))
  898. return Value::errorVALUE();
  899. // gaussinv (x)*sigma + mue
  900. return calc->add (calc->mul (calc->gaussinv (x), sigma), mue);
  901. }
  902. // Function: gammaln
  903. Value func_gammaln (valVector args, ValueCalc *calc, FuncExtra *) {
  904. //returns the natural logarithm of the gamma function
  905. if (calc->greater (args[0], 0.0))
  906. return calc->GetLogGamma (args[0]);
  907. return Value::errorVALUE();
  908. }
  909. // Function: poisson
  910. Value func_poisson (valVector args, ValueCalc *calc, FuncExtra *) {
  911. //returns the Poisson distribution
  912. Value x = args[0];
  913. Value lambda = args[1];
  914. Value kum = args[2];
  915. // lambda < 0.0 || x < 0.0
  916. if (calc->lower (lambda, 0.0) || calc->lower (x, 0.0))
  917. return Value::errorVALUE();
  918. Value result;
  919. // ex = exp (-lambda)
  920. Value ex = calc->exp (calc->mul (lambda, -1));
  921. if (calc->isZero (kum)) { // density
  922. if (calc->isZero (lambda))
  923. result = 0;
  924. else
  925. // ex * pow (lambda, x) / fact (x)
  926. result = calc->div (calc->mul (ex, calc->pow (lambda, x)), calc->fact (x));
  927. }
  928. else { // distribution
  929. if (calc->isZero (lambda))
  930. result = 1;
  931. else
  932. {
  933. result = 1.0;
  934. Value fFak = 1.0;
  935. unsigned long nEnd = calc->conv()->asInteger (x).asInteger();
  936. for (unsigned long i = 1; i <= nEnd; i++)
  937. {
  938. // fFak *= i
  939. fFak = calc->mul (fFak, i);
  940. // result += pow (lambda, i) / fFak
  941. result = calc->add (result, calc->div (calc->pow (lambda, i), fFak));
  942. }
  943. result = calc->mul (result, ex);
  944. }
  945. }
  946. return result;
  947. }
  948. // Function: confidence
  949. Value func_confidence (valVector args, ValueCalc *calc, FuncExtra *) {
  950. //returns the confidence interval for a population mean
  951. Value alpha = args[0];
  952. Value sigma = args[1];
  953. Value n = args[2];
  954. // sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1
  955. if ((!calc->greater (sigma, 0.0)) || (!calc->greater (alpha, 0.0)) ||
  956. (!calc->lower (alpha, 1.0)) || calc->lower (n, 1))
  957. return Value::errorVALUE();
  958. // g = gaussinv (1.0 - alpha / 2.0)
  959. Value g = calc->gaussinv (calc->sub (1.0, calc->div (alpha, 2.0)));
  960. // g * sigma / sqrt (n)
  961. return calc->div (calc->mul (g, sigma), calc->sqrt (n));
  962. }
  963. // Function: tdist
  964. Value func_tdist (valVector args, ValueCalc *calc, FuncExtra *) {
  965. //returns the t-distribution
  966. Value T = args[0];
  967. Value fDF = args[1];
  968. int flag = calc->conv()->asInteger (args[2]).asInteger();
  969. if (calc->lower (fDF, 1) || calc->lower (T, 0.0) || (flag != 1 && flag != 2))
  970. return Value::errorVALUE();
  971. // arg = fDF / (fDF + T * T)
  972. Value arg = calc->div (fDF, calc->add (fDF, calc->sqr (T)));
  973. Value R;
  974. R = calc->mul (calc->GetBeta (arg, calc->div (fDF, 2.0), 0.5), 0.5);
  975. if (flag == 1)
  976. return R;
  977. return calc->mul (R, 2); // flag is 2 here
  978. }
  979. // Function: fdist
  980. Value func_fdist (valVector args, ValueCalc *calc, FuncExtra *) {
  981. //returns the f-distribution
  982. Value x = args[0];
  983. Value fF1 = args[1];
  984. Value fF2 = args[2];
  985. // x < 0.0 || fF1 < 1 || fF2 < 1 || fF1 >= 1.0E10 || fF2 >= 1.0E10
  986. if (calc->lower (x, 0.0) || calc->lower (fF1, 1) || calc->lower (fF2, 1) ||
  987. (!calc->lower (fF1, 1.0E10)) || (!calc->lower (fF2, 1.0E10)))
  988. return Value::errorVALUE();
  989. // arg = fF2 / (fF2 + fF1 * x)
  990. Value arg = calc->div (fF2, calc->add (fF2, calc->mul (fF1, x)));
  991. // alpha = fF2/2.0
  992. Value alpha = calc->div (fF2, 2.0);
  993. // beta = fF1/2.0
  994. Value beta = calc->div (fF1, 2.0);
  995. return calc->GetBeta (arg, alpha, beta);
  996. }
  997. // Function: chidist
  998. Value func_chidist (valVector args, ValueCalc *calc, FuncExtra *) {
  999. //returns the chi-distribution
  1000. Value fChi = args[0];
  1001. Value fDF = args[1];
  1002. // fDF < 1 || fDF >= 1.0E5 || fChi < 0.0
  1003. if (calc->lower (fDF, 1) || (!calc->lower (fDF, 1.0E5)) ||
  1004. calc->lower (fChi, 0.0))
  1005. return Value::errorVALUE();
  1006. // 1.0 - GetGammaDist (fChi / 2.0, fDF / 2.0, 1.0)
  1007. return calc->sub (1.0, calc->GetGammaDist (calc->div (fChi, 2.0),
  1008. calc->div (fDF, 2.0), 1.0));
  1009. }
  1010. // two-array-walk functions used in the two-sum functions
  1011. void tawSumproduct (ValueCalc *c, Value &res, Value v1,
  1012. Value v2) {
  1013. // res += v1*v2
  1014. res = c->add (res, c->mul (v1, v2));
  1015. }
  1016. void tawSumx2py2 (ValueCalc *c, Value &res, Value v1,
  1017. Value v2) {
  1018. // res += sqr(v1)+sqr(v2)
  1019. res = c->add (res, c->add (c->sqr (v1), c->sqr (v2)));
  1020. }
  1021. void tawSumx2my2 (ValueCalc *c, Value &res, Value v1,
  1022. Value v2) {
  1023. // res += sqr(v1)-sqr(v2)
  1024. res = c->add (res, c->sub (c->sqr (v1), c->sqr (v2)));
  1025. }
  1026. void tawSumxmy2 (ValueCalc *c, Value &res, Value v1,
  1027. Value v2) {
  1028. // res += sqr(v1-v2)
  1029. res = c->add (res, c->sqr (c->sub (v1, v2)));
  1030. }
  1031. // Function: sumproduct
  1032. Value func_sumproduct (valVector args, ValueCalc *calc, FuncExtra *)
  1033. {
  1034. Value result;
  1035. calc->twoArrayWalk (args[0], args[1], result, tawSumproduct);
  1036. return result;
  1037. }
  1038. // Function: sumx2py2
  1039. Value func_sumx2py2 (valVector args, ValueCalc *calc, FuncExtra *)
  1040. {
  1041. Value result;
  1042. calc->twoArrayWalk (args[0], args[1], result, tawSumx2py2);
  1043. return result;
  1044. }
  1045. // Function: sumx2my2
  1046. Value func_sumx2my2 (valVector args, ValueCalc *calc, FuncExtra *)
  1047. {
  1048. Value result;
  1049. calc->twoArrayWalk (args[0], args[1], result, tawSumx2my2);
  1050. return result;
  1051. }
  1052. // Function: sum2xmy
  1053. Value func_sumxmy2 (valVector args, ValueCalc *calc, FuncExtra *)
  1054. {
  1055. Value result;
  1056. calc->twoArrayWalk (args[0], args[1], result, tawSumxmy2);
  1057. return result;
  1058. }