kspread

kspread_functions_statistical.cc

00001 /* This file is part of the KDE project
00002    Copyright (C) 1998-2002 The KSpread Team
00003                            www.koffice.org/kspread
00004    Copyright (C) 2005 Tomas Mecir <mecirt@gmail.com>
00005 
00006    This library is free software; you can redistribute it and/or
00007    modify it under the terms of the GNU Library General Public
00008    License as published by the Free Software Foundation; either
00009    version 2 of the License.
00010 
00011    This library is distributed in the hope that it will be useful,
00012    but WITHOUT ANY WARRANTY; without even the implied warranty of
00013    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00014    Library General Public License for more details.
00015 
00016    You should have received a copy of the GNU Library General Public License
00017    along with this library; see the file COPYING.LIB.  If not, write to
00018    the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
00019  * Boston, MA 02110-1301, USA.
00020 */
00021 
00022 // built-in statistical functions
00023 
00024 #include "functions.h"
00025 #include "valuecalc.h"
00026 #include "valueconverter.h"
00027 
00028 // needed for MODE
00029 #include <qmap.h>
00030 
00031 using namespace KSpread;
00032 
00033 // prototypes (sorted!)
00034 Value func_arrang (valVector args, ValueCalc *calc, FuncExtra *);
00035 Value func_average (valVector args, ValueCalc *calc, FuncExtra *);
00036 Value func_averagea (valVector args, ValueCalc *calc, FuncExtra *);
00037 Value func_avedev (valVector args, ValueCalc *calc, FuncExtra *);
00038 Value func_betadist (valVector args, ValueCalc *calc, FuncExtra *);
00039 Value func_bino (valVector args, ValueCalc *calc, FuncExtra *);
00040 Value func_chidist (valVector args, ValueCalc *calc, FuncExtra *);
00041 Value func_combin (valVector args, ValueCalc *calc, FuncExtra *);
00042 Value func_confidence (valVector args, ValueCalc *calc, FuncExtra *);
00043 Value func_correl_pop (valVector args, ValueCalc *calc, FuncExtra *);
00044 Value func_covar (valVector args, ValueCalc *calc, FuncExtra *);
00045 Value func_devsq (valVector args, ValueCalc *calc, FuncExtra *);
00046 Value func_devsqa (valVector args, ValueCalc *calc, FuncExtra *);
00047 Value func_expondist (valVector args, ValueCalc *calc, FuncExtra *);
00048 Value func_fdist (valVector args, ValueCalc *calc, FuncExtra *);
00049 Value func_fisher (valVector args, ValueCalc *calc, FuncExtra *);
00050 Value func_fisherinv (valVector args, ValueCalc *calc, FuncExtra *);
00051 Value func_gammadist (valVector args, ValueCalc *calc, FuncExtra *);
00052 Value func_gammaln (valVector args, ValueCalc *calc, FuncExtra *);
00053 Value func_gauss (valVector args, ValueCalc *calc, FuncExtra *);
00054 Value func_geomean (valVector args, ValueCalc *calc, FuncExtra *);
00055 Value func_harmean (valVector args, ValueCalc *calc, FuncExtra *);
00056 Value func_hypgeomdist (valVector args, ValueCalc *calc, FuncExtra *);
00057 Value func_kurtosis_est (valVector args, ValueCalc *calc, FuncExtra *);
00058 Value func_kurtosis_pop (valVector args, ValueCalc *calc, FuncExtra *);
00059 Value func_large (valVector args, ValueCalc *calc, FuncExtra *);
00060 Value func_loginv (valVector args, ValueCalc *calc, FuncExtra *);
00061 Value func_lognormdist (valVector args, ValueCalc *calc, FuncExtra *);
00062 Value func_median (valVector args, ValueCalc *calc, FuncExtra *);
00063 Value func_mode (valVector args, ValueCalc *calc, FuncExtra *);
00064 Value func_negbinomdist (valVector args, ValueCalc *calc, FuncExtra *);
00065 Value func_normdist (valVector args, ValueCalc *calc, FuncExtra *);
00066 Value func_norminv (valVector args, ValueCalc *calc, FuncExtra *);
00067 Value func_normsinv (valVector args, ValueCalc *calc, FuncExtra *);
00068 Value func_phi (valVector args, ValueCalc *calc, FuncExtra *);
00069 Value func_poisson (valVector args, ValueCalc *calc, FuncExtra *);
00070 Value func_skew_est (valVector args, ValueCalc *calc, FuncExtra *);
00071 Value func_skew_pop (valVector args, ValueCalc *calc, FuncExtra *);
00072 Value func_small (valVector args, ValueCalc *calc, FuncExtra *);
00073 Value func_standardize (valVector args, ValueCalc *calc, FuncExtra *);
00074 Value func_stddev (valVector args, ValueCalc *calc, FuncExtra *);
00075 Value func_stddeva (valVector args, ValueCalc *calc, FuncExtra *);
00076 Value func_stddevp (valVector args, ValueCalc *calc, FuncExtra *);
00077 Value func_stddevpa (valVector args, ValueCalc *calc, FuncExtra *);
00078 Value func_stdnormdist (valVector args, ValueCalc *calc, FuncExtra *);
00079 Value func_sumproduct (valVector args, ValueCalc *calc, FuncExtra *);
00080 Value func_sumx2py2 (valVector args, ValueCalc *calc, FuncExtra *);
00081 Value func_sumx2my2 (valVector args, ValueCalc *calc, FuncExtra *);
00082 Value func_sumxmy2 (valVector args, ValueCalc *calc, FuncExtra *);
00083 Value func_tdist (valVector args, ValueCalc *calc, FuncExtra *);
00084 Value func_variance (valVector args, ValueCalc *calc, FuncExtra *);
00085 Value func_variancea (valVector args, ValueCalc *calc, FuncExtra *);
00086 Value func_variancep (valVector args, ValueCalc *calc, FuncExtra *);
00087 Value func_variancepa (valVector args, ValueCalc *calc, FuncExtra *);
00088 Value func_weibull (valVector args, ValueCalc *calc, FuncExtra *);
00089 
00090 typedef QValueList<double> List;
00091 
00092 // registers all statistical functions
00093 void RegisterStatisticalFunctions()
00094 {
00095   FunctionRepository* repo = FunctionRepository::self();
00096   Function *f;
00097 
00098   f = new Function ("AVEDEV", func_avedev);
00099   f->setParamCount (1, -1);
00100   f->setAcceptArray ();
00101   repo->add (f);
00102   f = new Function ("AVERAGE", func_average);
00103   f->setParamCount (1, -1);
00104   f->setAcceptArray ();
00105   repo->add (f);
00106   f = new Function ("AVERAGEA", func_averagea);
00107   f->setParamCount (1, -1);
00108   f->setAcceptArray ();
00109   repo->add (f);
00110   f = new Function ("BETADIST", func_betadist);
00111   f->setParamCount (3, 5);
00112   repo->add (f);
00113   f = new Function ("BINO", func_bino);
00114   f->setParamCount (3);
00115   repo->add (f);
00116   f = new Function ("CHIDIST", func_chidist);
00117   f->setParamCount (2);
00118   repo->add (f);
00119   f = new Function ("COMBIN", func_combin);
00120   f->setParamCount (2);
00121   repo->add (f);
00122   f = new Function ("CONFIDENCE", func_confidence);
00123   f->setParamCount (3);
00124   repo->add (f);
00125   f = new Function ("CORREL", func_correl_pop);
00126   f->setParamCount (2);
00127   f->setAcceptArray ();
00128   repo->add (f);
00129   f = new Function ("COVAR", func_covar);
00130   f->setParamCount (2);
00131   f->setAcceptArray ();
00132   repo->add (f);
00133   f = new Function ("DEVSQ", func_devsq);
00134   f->setParamCount (1, -1);
00135   f->setAcceptArray ();
00136   repo->add (f);
00137   f = new Function ("DEVSQA", func_devsqa);
00138   f->setParamCount (1, -1);
00139   f->setAcceptArray ();
00140   repo->add (f);
00141   f = new Function ("EXPONDIST", func_expondist);
00142   f->setParamCount (3);
00143   repo->add (f);
00144   f = new Function ("FDIST", func_fdist);
00145   f->setParamCount (3);
00146   repo->add (f);
00147   f = new Function ("FISHER", func_fisher);
00148   repo->add (f);
00149   f = new Function ("FISHERINV", func_fisherinv);
00150   repo->add (f);
00151   f = new Function ("GAMMADIST", func_gammadist);
00152   f->setParamCount (4);
00153   repo->add (f);
00154   f = new Function ("GAMMALN", func_gammaln);
00155   repo->add (f);
00156   f = new Function ("GAUSS", func_gauss);
00157   repo->add (f);
00158   f = new Function ("GEOMEAN", func_geomean);
00159   f->setParamCount (1, -1);
00160   f->setAcceptArray ();
00161   repo->add (f);
00162   f = new Function ("HARMEAN", func_harmean);
00163   f->setParamCount (1, -1);
00164   f->setAcceptArray ();
00165   repo->add (f);
00166   f = new Function ("HYPGEOMDIST", func_hypgeomdist);
00167   f->setParamCount (4);
00168   repo->add (f);
00169   f = new Function ("INVBINO", func_bino);  // same as BINO, for 1.4 compat
00170   repo->add (f);
00171   f = new Function ("KURT", func_kurtosis_est);
00172   f->setParamCount (1, -1);
00173   f->setAcceptArray ();
00174   repo->add (f);
00175   f = new Function ("KURTP", func_kurtosis_pop);
00176   f->setParamCount (1, -1);
00177   f->setAcceptArray ();
00178   repo->add (f);
00179   f = new Function ("LARGE", func_large);
00180   f->setParamCount (2);
00181   f->setAcceptArray ();
00182   repo->add (f);
00183   f = new Function ("LOGINV", func_loginv);
00184   f->setParamCount (3);
00185   repo->add (f);
00186   f = new Function ("LOGNORMDIST", func_lognormdist);
00187   f->setParamCount (3);
00188   repo->add (f);
00189   f = new Function ("MEDIAN", func_median);
00190   f->setParamCount (1, -1);
00191   f->setAcceptArray ();
00192   repo->add (f);
00193   f = new Function ("MODE", func_mode);
00194   f->setParamCount (1, -1);
00195   f->setAcceptArray ();
00196   repo->add (f);
00197   f = new Function ("NEGBINOMDIST", func_negbinomdist);
00198   f->setParamCount (3);
00199   repo->add (f);
00200   f = new Function ("NORMDIST", func_normdist);
00201   f->setParamCount (4);
00202   repo->add (f);
00203   f = new Function ("NORMINV", func_norminv);
00204   f->setParamCount (3);
00205   repo->add (f);
00206   f = new Function ("NORMSDIST", func_stdnormdist);
00207   repo->add (f);
00208   f = new Function ("NORMSINV", func_normsinv);
00209   repo->add (f);
00210   f = new Function ("PEARSON", func_correl_pop);
00211   f->setParamCount (2);
00212   f->setAcceptArray ();
00213   repo->add (f);
00214   f = new Function ("PERMUT", func_arrang);
00215   f->setParamCount (2);
00216   repo->add (f);
00217   f = new Function ("PHI", func_phi);
00218   repo->add (f);
00219   f = new Function ("POISSON", func_poisson);
00220   f->setParamCount (3);
00221   repo->add (f);
00222   f = new Function ("SKEW", func_skew_est);
00223   f->setParamCount (1, -1);
00224   f->setAcceptArray ();
00225   repo->add (f);
00226   f = new Function ("SKEWP", func_skew_pop);
00227   f->setParamCount (1, -1);
00228   f->setAcceptArray ();
00229   repo->add (f);
00230   f = new Function ("SMALL", func_small);
00231   f->setParamCount (2);
00232   f->setAcceptArray ();
00233   repo->add (f);
00234   f = new Function ("STANDARDIZE", func_standardize);
00235   f->setParamCount (3);
00236   repo->add (f);
00237   f = new Function ("STDEV", func_stddev);
00238   f->setParamCount (1, -1);
00239   f->setAcceptArray ();
00240   repo->add (f);
00241   f = new Function ("STDEVA", func_stddeva);
00242   f->setParamCount (1, -1);
00243   f->setAcceptArray ();
00244   repo->add (f);
00245   f = new Function ("STDEVP", func_stddevp);
00246   f->setParamCount (1, -1);
00247   f->setAcceptArray ();
00248   repo->add (f);
00249   f = new Function ("STDEVPA", func_stddevpa);
00250   f->setParamCount (1, -1);
00251   f->setAcceptArray ();
00252   repo->add (f);
00253   f = new Function ("SUM2XMY", func_sumxmy2);
00254   f->setParamCount (2);
00255   f->setAcceptArray ();
00256   repo->add (f);
00257   f = new Function ("SUMPRODUCT", func_sumproduct);
00258   f->setParamCount (2);
00259   f->setAcceptArray ();
00260   repo->add (f);
00261   f = new Function ("SUMX2PY2", func_sumx2py2);
00262   f->setParamCount (2);
00263   f->setAcceptArray ();
00264   repo->add (f);
00265   f = new Function ("SUMX2MY2", func_sumx2my2);
00266   f->setParamCount (2);
00267   f->setAcceptArray ();
00268   repo->add (f);
00269   f = new Function ("TDIST", func_tdist);
00270   f->setParamCount (3);
00271   repo->add (f);
00272   f = new Function ("VARIANCE", func_variance);
00273   f->setParamCount (1, -1);
00274   f->setAcceptArray ();
00275   repo->add (f);
00276   f = new Function ("VAR", func_variance);
00277   f->setParamCount (1, -1);
00278   f->setAcceptArray ();
00279   repo->add (f);
00280   f = new Function ("VARP", func_variancep);
00281   f->setParamCount (1, -1);
00282   f->setAcceptArray ();
00283   repo->add (f);
00284   f = new Function ("VARA", func_variancea);
00285   f->setParamCount (1, -1);
00286   f->setAcceptArray ();
00287   repo->add (f);
00288   f = new Function ("VARPA", func_variancepa);
00289   f->setParamCount (1, -1);
00290   f->setAcceptArray ();
00291   repo->add (f);
00292   f = new Function ("WEIBULL", func_weibull);
00293   f->setParamCount (4);
00294   repo->add (f);
00295 }
00296 
00297 // array-walk functions used in this file
00298 
00299 void awSkew (ValueCalc *c, Value &res, Value val, Value p)
00300 {
00301   Value avg = p.element (0, 0);
00302   Value stdev = p.element (1, 0);
00303   // (val - avg) / stddev
00304   Value d = c->div (c->sub (val, avg), stdev);
00305   // res += d*d*d
00306   res = c->add (res, c->mul (d, c->mul (d, d)));
00307 }
00308 
00309 void awSumInv (ValueCalc *c, Value &res, Value val, Value)
00310 {
00311   // res += 1/value
00312   res = c->add (res, c->div (1.0, val));
00313 }
00314 
00315 void awAveDev (ValueCalc *c, Value &res, Value val,
00316     Value p)
00317 {
00318   // res += abs (val - p)
00319   res = c->add (res, c->abs (c->sub (val, p)));
00320 }
00321 
00322 void awKurtosis (ValueCalc *c, Value &res, Value val,
00323     Value p)
00324 {
00325   Value avg = p.element (0, 0);
00326   Value stdev = p.element (1, 0);
00327   //d = (val - avg ) / stdev
00328   Value d = c->div (c->sub (val, avg), stdev);
00329   // res += d^4
00330   res = c->add (res, c->pow (d, 4));
00331 }
00332 
00333 
00334 Value func_skew_est (valVector args, ValueCalc *calc, FuncExtra *)
00335 {
00336   int number = calc->count (args);
00337   Value avg = calc->avg (args);
00338   if (number < 3)
00339     return Value::errorVALUE();
00340 
00341   Value res = calc->stddev (args, avg);
00342   if (res.isZero())
00343     return Value::errorVALUE();
00344   
00345   Value params (2, 1);
00346   params.setElement (0, 0, avg);
00347   params.setElement (1, 0, res);
00348   Value tskew;
00349   calc->arrayWalk (args, tskew, awSkew, params);
00350 
00351   // ((tskew * number) / (number-1)) / (number-2)
00352   return calc->div (calc->div (calc->mul (tskew, number), number-1), number-2);
00353 }
00354 
00355 Value func_skew_pop (valVector args, ValueCalc *calc, FuncExtra *)
00356 {
00357   int number = calc->count (args);
00358   Value avg = calc->avg (args);
00359   if (number < 1)
00360     return Value::errorVALUE();
00361   
00362   Value res = calc->stddevP (args, avg);
00363   if (res.isZero())
00364     return Value::errorVALUE();
00365   
00366   Value params (2, 1);
00367   params.setElement (0, 0, avg);
00368   params.setElement (1, 0, res);
00369   Value tskew;
00370   calc->arrayWalk (args, tskew, awSkew, params);
00371   
00372   // tskew / number
00373   return calc->div (tskew, number);
00374 }
00375 
00376 class ContentSheet : public QMap<double, int> {};
00377 
00378 void func_mode_helper (Value range, ValueCalc *calc, ContentSheet &sh)
00379 {
00380   if (!range.isArray())
00381   {
00382     double d = calc->conv()->asFloat (range).asFloat();
00383     sh[d]++;
00384     return;
00385   }
00386   
00387   for (unsigned int row = 0; row < range.rows(); ++row)
00388     for (unsigned int col = 0; col < range.columns(); ++col) {
00389       Value v = range.element (col, row);
00390       if (v.isArray())
00391         func_mode_helper (v, calc, sh);
00392       else {
00393         double d = calc->conv()->asFloat (v).asFloat();
00394         sh[d]++;
00395       }
00396     }
00397 }
00398 
00399 Value func_mode (valVector args, ValueCalc *calc, FuncExtra *)
00400 {
00401   // does NOT support anything other than doubles !!!
00402   ContentSheet sh;
00403   for (unsigned int i = 0; i < args.count(); ++i)
00404     func_mode_helper (args[i], calc, sh);
00405   
00406   // retrieve value with max.count
00407   int maxcount = 0;
00408   double max = 0.0;
00409   ContentSheet::iterator it;
00410   for (it = sh.begin(); it != sh.end(); ++it)
00411     if (it.data() > maxcount) {
00412       max = it.key();
00413       maxcount = it.data();
00414     }
00415   return Value (max);
00416 }
00417 
00418 Value func_covar_helper (Value range1, Value range2,
00419     ValueCalc *calc, Value avg1, Value avg2)
00420 {
00421   // two arrays -> cannot use arrayWalk
00422   if ((!range1.isArray()) && (!range2.isArray()))
00423     // (v1-E1)*(v2-E2)
00424     return calc->mul (calc->sub (range1, avg1), calc->sub (range2, avg2));
00425   
00426   int rows = range1.rows();
00427   int cols = range1.columns();
00428   int rows2 = range2.rows();
00429   int cols2 = range2.columns();
00430   if ((rows != rows2) || (cols != cols2))
00431     return Value::errorVALUE();
00432     
00433   Value result = 0.0;
00434   for (int row = 0; row < rows; ++row)
00435     for (int col = 0; col < cols; ++col) {
00436       Value v1 = range1.element (col, row);
00437       Value v2 = range2.element (col, row);
00438       if (v1.isArray() || v2.isArray())
00439         result = calc->add (result,
00440             func_covar_helper (v1, v2, calc, avg1, avg2));
00441       else
00442         // result += (v1-E1)*(v2-E2)
00443         result = calc->add (result, calc->mul (calc->sub (v1, avg1),
00444             calc->sub (v2, avg2)));
00445     }
00446 
00447   return result;
00448 }
00449 
00450 Value func_covar (valVector args, ValueCalc *calc, FuncExtra *)
00451 {
00452   Value avg1 = calc->avg (args[0]);
00453   Value avg2 = calc->avg (args[1]);
00454   int number = calc->count (args[0]);
00455   int number2 = calc->count (args[1]);
00456   
00457   if (number2 <= 0 || number2 != number)
00458     return Value::errorVALUE();
00459 
00460   Value covar = func_covar_helper (args[0], args[1], calc, avg1, avg2);
00461   return calc->div (covar, number);
00462 }
00463 
00464 Value func_correl_pop (valVector args, ValueCalc *calc, FuncExtra *)
00465 {
00466   Value covar = func_covar (args, calc, 0);
00467   Value stdevp1 = calc->stddevP (args[0]);
00468   Value stdevp2 = calc->stddevP (args[1]);
00469   
00470   if (calc->isZero (stdevp1) || calc->isZero (stdevp2))
00471     return Value::errorDIV0();
00472 
00473   // covar / (stdevp1 * stdevp2)
00474   return calc->div (covar, calc->mul (stdevp1, stdevp2));
00475 }
00476 
00477 void func_array_helper (Value range, ValueCalc *calc,
00478     List &array, int &number)
00479 {
00480   if (!range.isArray())
00481   {
00482     array << calc->conv()->asFloat (range).asFloat();
00483     ++number;
00484     return;
00485   }
00486   
00487   for (unsigned int row = 0; row < range.rows(); ++row)
00488     for (unsigned int col = 0; col < range.columns(); ++col) {
00489       Value v = range.element (col, row);
00490       if (v.isArray ())
00491         func_array_helper (v, calc, array, number);
00492       else {
00493         array << calc->conv()->asFloat (v).asFloat();
00494         ++number;
00495       }
00496     }
00497 }
00498 
00499 Value func_large (valVector args, ValueCalc *calc, FuncExtra *)
00500 {
00501   // does NOT support anything other than doubles !!!
00502   int k = calc->conv()->asInteger (args[1]).asInteger();
00503   if ( k < 1 )
00504     return false;
00505 
00506   List array;
00507   int number = 1;
00508   
00509   func_array_helper (args[0], calc, array, number);
00510 
00511   if ( k > number )
00512     return Value::errorVALUE();
00513 
00514   qHeapSort (array);
00515   double d = *array.at (number - k - 1);
00516   return Value (d);
00517 }
00518 
00519 Value func_small (valVector args, ValueCalc *calc, FuncExtra *)
00520 {
00521   // does NOT support anything other than doubles !!!
00522   int k = calc->conv()->asInteger (args[1]).asInteger();
00523   if ( k < 1 )
00524     return false;
00525 
00526   List array;
00527   int number = 1;
00528   
00529   func_array_helper (args[0], calc, array, number);
00530 
00531   if ( k > number )
00532     return Value::errorVALUE();
00533 
00534   qHeapSort (array);
00535   double d = *array.at (k - 1);
00536   return Value (d);
00537 }
00538 
00539 Value func_geomean (valVector args, ValueCalc *calc, FuncExtra *)
00540 {
00541   Value count = calc->count (args);
00542   Value prod = calc->product (args, 1.0);
00543   if (calc->isZero (count))
00544     return Value::errorDIV0();
00545   return calc->pow (prod, calc->div (1.0, count));
00546 }
00547 
00548 Value func_harmean (valVector args, ValueCalc *calc, FuncExtra *)
00549 {
00550   Value count = calc->count (args);
00551   if (calc->isZero (count))
00552     return Value::errorDIV0();
00553   Value suminv;
00554   calc->arrayWalk (args, suminv, awSumInv, 0);
00555   return calc->div (suminv, count);
00556 }
00557 
00558 Value func_loginv (valVector args, ValueCalc *calc, FuncExtra *)
00559 {
00560   Value p = args[0];
00561   Value m = args[1];
00562   Value s = args[2];
00563 
00564   if (calc->lower (p, 0) || calc->greater (p, 1))
00565     return Value::errorVALUE();
00566   
00567   if (!calc->greater (s, 0))
00568     return Value::errorVALUE();
00569 
00570   Value result = 0.0;
00571   if (calc->equal (p, 1))   //p==1
00572     result = Value::errorVALUE();
00573   else if (calc->greater (p, 0)) {   //p>0
00574     Value gaussInv = calc->gaussinv (p);
00575     // exp (gaussInv * s + m)
00576     result = calc->exp (calc->add (calc->mul (s, gaussInv), m));
00577   }
00578 
00579   return result;
00580 }
00581 
00582 Value func_devsq (valVector args, ValueCalc *calc, FuncExtra *)
00583 {
00584   Value res;
00585   calc->arrayWalk (args, res, calc->awFunc ("devsq"), calc->avg (args, false));
00586   return res;
00587 }
00588 
00589 Value func_devsqa (valVector args, ValueCalc *calc, FuncExtra *)
00590 {
00591   Value res;
00592   calc->arrayWalk (args, res, calc->awFunc ("devsqa"), calc->avg (args));
00593   return res;
00594 }
00595 
00596 Value func_kurtosis_est (valVector args, ValueCalc *calc, FuncExtra *)
00597 {
00598   int count = calc->count (args);
00599   if (count < 4)
00600     return Value::errorVALUE();
00601 
00602   Value avg = calc->avg (args);
00603   Value devsq;
00604   calc->arrayWalk (args, devsq, calc->awFunc ("devsqa"), avg);
00605 
00606   if (devsq.isZero ())
00607     return Value::errorDIV0();
00608 
00609   Value params (2, 1);
00610   params.setElement (0, 0, avg);
00611   params.setElement (1, 0, devsq);
00612   Value x4;
00613   calc->arrayWalk (args, x4, awKurtosis, params);
00614 
00615   double den = (double) (count - 2) * (count - 3);
00616   double nth = (double) count * (count + 1) / ((count - 1) * den);
00617   double t = 3.0 * (count - 1) * (count - 1) / den;
00618 
00619   // res = x4 * nth - t
00620   return calc->sub (calc->mul (x4, nth), t);
00621 }
00622 
00623 Value func_kurtosis_pop (valVector args, ValueCalc *calc, FuncExtra *)
00624 {
00625   int count = calc->count (args);
00626   if (count < 4)
00627     return Value::errorVALUE();
00628 
00629   Value avg = calc->avg (args);
00630   Value devsq;
00631   calc->arrayWalk (args, devsq, calc->awFunc ("devsqa"), avg);
00632 
00633   if (devsq.isZero ())
00634     return Value::errorDIV0();
00635 
00636   Value params (2, 1);
00637   params.setElement (0, 0, avg);
00638   params.setElement (1, 0, devsq);
00639   Value x4;
00640   calc->arrayWalk (args, x4, awKurtosis, params);
00641   
00642   // x4 / count - 3
00643   return calc->sub (calc->div (x4, count), 3);
00644 }
00645 
00646 Value func_standardize (valVector args, ValueCalc *calc, FuncExtra *)
00647 {
00648   Value x = args[0];
00649   Value m = args[1];
00650   Value s = args[2];
00651 
00652   if (!calc->greater (s, 0))  // s must be >0
00653     return Value::errorVALUE();
00654 
00655   // (x - m) / s
00656   return calc->div (calc->sub (x, m), s);
00657 }
00658 
00659 Value func_hypgeomdist (valVector args, ValueCalc *calc, FuncExtra *)
00660 {
00661   int x = calc->conv()->asInteger (args[0]).asInteger();
00662   int n = calc->conv()->asInteger (args[1]).asInteger();
00663   int M = calc->conv()->asInteger (args[2]).asInteger();
00664   int N = calc->conv()->asInteger (args[3]).asInteger();
00665 
00666   if ( x < 0 || n < 0 || M < 0 || N < 0 )
00667     return Value::errorVALUE();
00668 
00669   if ( x > M || n > N )
00670     return Value::errorVALUE();
00671 
00672   Value d1 = calc->combin (M, x);
00673   Value d2 = calc->combin (N - M, n - x);
00674   Value d3 = calc->combin (N, n);
00675 
00676   // d1 * d2 / d3
00677   return calc->div (calc->mul (d1, d2), d3);
00678 }
00679 
00680 Value func_negbinomdist (valVector args, ValueCalc *calc, FuncExtra *)
00681 {
00682   int x = calc->conv()->asInteger (args[0]).asInteger();
00683   int r = calc->conv()->asInteger (args[1]).asInteger();
00684   Value p = args[2];
00685 
00686   if ((x + r - 1) <= 0)
00687     return Value::errorVALUE();
00688   if (calc->lower (p, 0) || calc->greater (p, 1))
00689     return Value::errorVALUE();
00690 
00691   Value d1 = calc->combin (x + r - 1, r - 1);
00692   // d2 = pow (p, r) * pow (1 - p, x)
00693   Value d2 = calc->mul (calc->pow (p, r),
00694       calc->pow (calc->sub (1, p), x));
00695 
00696   return calc->mul (d1, d2);
00697 }
00698 
00699 // Function: permut
00700 Value func_arrang (valVector args, ValueCalc *calc, FuncExtra *)
00701 { 
00702   Value n = args[0];
00703   Value m = args[1];
00704   if (calc->lower (n, m))  // problem if n<m
00705     return Value::errorVALUE();
00706 
00707   if (calc->lower (m, 0))  // problem if m<0  (n>=m so that's okay)
00708     return Value::errorVALUE();
00709 
00710   // fact(n) / (fact(n-m)
00711   return calc->fact (n, calc->sub (n, m));
00712 }
00713 
00714 // Function: average
00715 Value func_average (valVector args, ValueCalc *calc, FuncExtra *)
00716 {
00717   return calc->avg (args, false);
00718 }
00719 
00720 // Function: averagea
00721 Value func_averagea (valVector args, ValueCalc *calc, FuncExtra *)
00722 {
00723   return calc->avg (args);
00724 }
00725 
00726 // Function: avedev
00727 Value func_avedev (valVector args, ValueCalc *calc, FuncExtra *)
00728 {
00729   Value result;
00730   calc->arrayWalk (args, result, awAveDev, calc->avg (args));
00731   return result;
00732 }
00733 
00734 // Function: median
00735 Value func_median (valVector args, ValueCalc *calc, FuncExtra *)
00736 {
00737   // does NOT support anything other than doubles !!!
00738   List array;
00739   int number = 1;
00740   
00741   for (unsigned int i = 0; i < args.count(); ++i)
00742     func_array_helper (args[i], calc, array, number);
00743 
00744   qHeapSort (array);
00745   double d = *array.at (number / 2 + number % 2);
00746   return Value (d);
00747 }
00748 
00749 // Function: variance
00750 Value func_variance (valVector args, ValueCalc *calc, FuncExtra *)
00751 {
00752   int count = calc->count (args, false);
00753   if (count < 2)
00754     return Value::errorVALUE();
00755 
00756   Value result = func_devsq (args, calc, 0);
00757   return calc->div (result, count-1);
00758 }
00759 
00760 // Function: vara
00761 Value func_variancea (valVector args, ValueCalc *calc, FuncExtra *)
00762 {
00763   int count = calc->count (args);
00764   if (count < 2)
00765     return Value::errorVALUE();
00766 
00767   Value result = func_devsqa (args, calc, 0);
00768   return calc->div (result, count-1);
00769 }
00770 
00771 // Function: varp
00772 Value func_variancep (valVector args, ValueCalc *calc, FuncExtra *)
00773 {
00774   int count = calc->count (args, false);
00775   if (count == 0)
00776     return Value::errorVALUE();
00777 
00778   Value result = func_devsq (args, calc, 0);
00779   return calc->div (result, count);
00780 }
00781 
00782 // Function: varpa
00783 Value func_variancepa (valVector args, ValueCalc *calc, FuncExtra *)
00784 {
00785   int count = calc->count (args);
00786   if (count == 0)
00787     return Value::errorVALUE();
00788 
00789   Value result = func_devsqa (args, calc, 0);
00790   return calc->div (result, count);
00791 }
00792 
00793 // Function: stddev
00794 Value func_stddev (valVector args, ValueCalc *calc, FuncExtra *)
00795 {
00796   return calc->stddev (args, false);
00797 }
00798 
00799 // Function: stddeva
00800 Value func_stddeva (valVector args, ValueCalc *calc, FuncExtra *)
00801 {
00802   return calc->stddev (args);
00803 }
00804 
00805 // Function: stddevp
00806 Value func_stddevp (valVector args, ValueCalc *calc, FuncExtra *)
00807 {
00808   return calc->stddevP (args, false);
00809 }
00810 
00811 // Function: stddevpa
00812 Value func_stddevpa (valVector args, ValueCalc *calc, FuncExtra *)
00813 {
00814   return calc->stddevP (args);
00815 }
00816 
00817 // Function: combin
00818 Value func_combin (valVector args, ValueCalc *calc, FuncExtra *)
00819 {
00820   return calc->combin (args[0], args[1]);
00821 }
00822 
00823 // Function: bino
00824 Value func_bino (valVector args, ValueCalc *calc, FuncExtra *)
00825 {
00826   Value n = args[0];
00827   Value m = args[1];
00828   Value comb = calc->combin (n, m);
00829   Value prob = args[2];
00830   
00831   if (calc->lower (prob,0) || calc->greater (prob,1))
00832     return Value::errorVALUE();
00833     
00834   // result = comb * pow (prob, m) * pow (1 - prob, n - m)
00835   Value pow1 = calc->pow (prob, m);
00836   Value pow2 = calc->pow (calc->sub (1, prob), calc->sub (n, m));
00837   return calc->mul (comb, calc->mul (pow1, pow2));
00838 }
00839 
00840 // Function: phi
00841 Value func_phi (valVector args, ValueCalc *calc, FuncExtra *)
00842 //distribution function for a standard normal distribution
00843 {
00844   return calc->phi (args[0]);
00845 }
00846 
00847 // Function: gauss
00848 Value func_gauss (valVector args, ValueCalc *calc, FuncExtra *)
00849 {
00850   //returns the integral values of the standard normal cumulative distribution
00851   return calc->gauss (args[0]);
00852 }
00853 
00854 // Function: gammadist
00855 Value func_gammadist (valVector args, ValueCalc *calc, FuncExtra *)
00856 {
00857   Value x = args[0];
00858   Value alpha = args[1];
00859   Value beta = args[2];
00860   int kum = calc->conv()->asInteger (args[3]).asInteger();  // 0 or 1
00861   
00862   Value result;
00863 
00864   if (calc->lower (x, 0.0) || (!calc->greater (alpha, 0.0)) ||
00865       (!calc->greater (beta, 0.0)))
00866     return Value::errorVALUE();
00867   
00868   if (kum == 0) {  //density
00869     Value G = calc->GetGamma (alpha);
00870     // result = pow (x, alpha - 1.0) / exp (x / beta) / pow (beta, alpha) / G
00871     Value pow1 = calc->pow (x, calc->sub (alpha, 1.0));
00872     Value pow2 = calc->exp (calc->div (x, beta));
00873     Value pow3 = calc->pow (beta, alpha);
00874     result = calc->div (calc->div (calc->div (pow1, pow2), pow3), G);
00875   }
00876   else
00877     result = calc->GetGammaDist (x, alpha, beta);
00878 
00879   return Value (result);
00880 }
00881 
00882 // Function: betadist
00883 Value func_betadist (valVector args, ValueCalc *calc, FuncExtra *)
00884 {
00885   Value x = args[0];
00886   Value alpha = args[1];
00887   Value beta = args[2];
00888   
00889   Value fA = 0.0;
00890   Value fB = 1.0;
00891   if (args.count() > 3) fA = args[3];
00892   if (args.count() == 5) fB = args[4];
00893 
00894   //x < fA || x > fB || fA == fB || alpha <= 0.0 || beta <= 0.0
00895   if (calc->lower (x, fA) || calc->greater (x, fB) || calc->equal (fA, fB) ||
00896       (!calc->greater (alpha, 0.0)) || (!calc->greater (beta, 0.0)))
00897     return Value::errorVALUE();
00898   
00899   // xx = (x - fA) / (fB - fA)  // scaling
00900   Value xx = calc->div (calc->sub (x, fA), calc->sub (fB, fA));
00901 
00902   return calc->GetBeta (xx, alpha, beta);
00903 }
00904 
00905 // Function: fisher
00906 Value func_fisher (valVector args, ValueCalc *calc, FuncExtra *) {
00907   //returns the Fisher transformation for x
00908 
00909   // 0.5 * ln ((1.0 + fVal) / (1.0 - fVal))
00910   Value fVal = args[0];
00911   Value num = calc->div (calc->add (fVal, 1.0), calc->sub (1.0, fVal));
00912   return calc->mul (calc->ln (num), 0.5);
00913 }
00914 
00915 // Function: fisherinv
00916 Value func_fisherinv (valVector args, ValueCalc *calc, FuncExtra *) {
00917   //returns the inverse of the Fisher transformation for x
00918   
00919   Value fVal = args[0];
00920   // (exp (2.0 * fVal) - 1.0) / (exp (2.0 * fVal) + 1.0)
00921   Value ex = calc->exp (calc->mul (fVal, 2.0));
00922   return calc->div (calc->sub (ex, 1.0), calc->add (ex, 1.0));
00923 }
00924 
00925 // Function: normdist
00926 Value func_normdist (valVector args, ValueCalc *calc, FuncExtra *) {
00927   //returns the normal cumulative distribution
00928   Value x = args[0];
00929   Value mue = args[1];
00930   Value sigma = args[2];
00931   Value k = args[3];
00932 
00933   if (!calc->greater (sigma, 0.0))
00934     return Value::errorVALUE();
00935 
00936   // (x - mue) / sigma
00937   Value Y = calc->div (calc->sub (x, mue), sigma);
00938   if (calc->isZero (k))   // density
00939     return calc->div (calc->phi (Y), sigma);
00940   else          // distribution
00941     return calc->add (calc->gauss (Y), 0.5);
00942 }
00943 
00944 // Function: lognormdist
00945 Value func_lognormdist (valVector args, ValueCalc *calc, FuncExtra *) {
00946   //returns the cumulative lognormal distribution
00947   Value x = args[0];
00948   Value mue = args[1];
00949   Value sigma = args[2];
00950 
00951   if (!calc->greater (sigma, 0.0) || (!calc->greater (x, 0.0)))
00952     return Value::errorVALUE();
00953   
00954   // (ln(x) - mue) / sigma
00955   Value Y = calc->div (calc->sub (calc->ln (x), mue), sigma);
00956   return calc->add (calc->gauss (Y), 0.5);
00957 }
00958 
00959 // Function: normsdist
00960 Value func_stdnormdist (valVector args, ValueCalc *calc, FuncExtra *)
00961 {
00962   //returns the cumulative lognormal distribution, mue=0, sigma=1
00963   return calc->add (calc->gauss (args[0]), 0.5);
00964 }
00965 
00966 // Function: expondist
00967 Value func_expondist (valVector args, ValueCalc *calc, FuncExtra *) {
00968   //returns the exponential distribution
00969   Value x = args[0];
00970   Value lambda = args[1];
00971   Value kum = args[2];
00972 
00973   Value result = 0.0;
00974 
00975   if (!calc->greater (lambda, 0.0))
00976     return Value::errorVALUE();
00977   
00978   // ex = exp (-lambda * x)
00979   Value ex = calc->exp (calc->mul (calc->mul (lambda, -1), x));
00980   if (calc->isZero (kum)) {  //density
00981     if (!calc->lower (x, 0.0))
00982       // lambda * ex
00983       result = calc->mul (lambda, ex);
00984   }
00985   else {  //distribution
00986     if (calc->greater (x, 0.0))
00987       // 1.0 - ex
00988       result = calc->sub (1.0, ex);
00989   }
00990   return result;
00991 }
00992 
00993 // Function: weibull
00994 Value func_weibull (valVector args, ValueCalc *calc, FuncExtra *) {
00995   //returns the Weibull distribution
00996 
00997   Value x = args[0];
00998   Value alpha = args[1];
00999   Value beta = args[2];
01000   Value kum = args[3];
01001 
01002   Value result;
01003 
01004   if ((!calc->greater (alpha, 0.0)) || (!calc->greater (beta, 0.0)) ||
01005       calc->lower (x, 0.0))
01006     return Value::errorVALUE();
01007   
01008   // ex = exp (-pow (x / beta, alpha))
01009   Value ex;
01010   ex = calc->exp (calc->mul (calc->pow (calc->div (x, beta), alpha), -1));
01011   if (calc->isZero (kum))    // density
01012   {
01013     // result = alpha / pow(beta,alpha) * pow(x,alpha-1.0) * ex
01014     result = calc->div (alpha, calc->pow (beta, alpha));
01015     result = calc->mul (result, calc->mul (calc->pow (x,
01016         calc->sub (alpha, 1)), ex));
01017   }
01018   else    // distribution
01019     result = calc->sub (1.0, ex);
01020 
01021   return result;
01022 }
01023 
01024 // Function: normsinv
01025 Value func_normsinv (valVector args, ValueCalc *calc, FuncExtra *) {
01026   //returns the inverse of the standard normal cumulative distribution
01027   
01028   Value x = args[0];
01029   if (!(calc->greater (x, 0.0) && calc->lower (x, 1.0)))
01030     return Value::errorVALUE();
01031   
01032   return calc->gaussinv (x);
01033 }
01034 
01035 // Function: norminv
01036 Value func_norminv (valVector args, ValueCalc *calc, FuncExtra *) {
01037   //returns the inverse of the normal cumulative distribution
01038   Value x = args[0];
01039   Value mue = args[1];
01040   Value sigma = args[2];
01041 
01042   if (!calc->greater (sigma, 0.0))
01043     return Value::errorVALUE();
01044   if (!(calc->greater (x, 0.0) && calc->lower (x, 1.0)))
01045     return Value::errorVALUE();
01046 
01047   // gaussinv (x)*sigma + mue
01048   return calc->add (calc->mul (calc->gaussinv (x), sigma), mue);
01049 }
01050 
01051 // Function: gammaln
01052 Value func_gammaln (valVector args, ValueCalc *calc, FuncExtra *) {
01053   //returns the natural logarithm of the gamma function
01054   
01055   if (calc->greater (args[0], 0.0))
01056     return calc->GetLogGamma (args[0]);
01057   return Value::errorVALUE();
01058 }
01059 
01060 // Function: poisson
01061 Value func_poisson (valVector args, ValueCalc *calc, FuncExtra *) {
01062   //returns the Poisson distribution
01063 
01064   Value x = args[0];
01065   Value lambda = args[1];
01066   Value kum = args[2];
01067 
01068   // lambda < 0.0 || x < 0.0
01069   if (calc->lower (lambda, 0.0) || calc->lower (x, 0.0))
01070     return Value::errorVALUE();
01071   
01072   Value result;
01073 
01074   // ex = exp (-lambda)
01075   Value ex = calc->exp (calc->mul (lambda, -1));
01076 
01077   if (calc->isZero (kum)) {   // density
01078     if (calc->isZero (lambda))
01079       result = 0;
01080     else
01081       // ex * pow (lambda, x) / fact (x)
01082     result = calc->div (calc->mul (ex, calc->pow (lambda, x)), calc->fact (x));
01083   }
01084   else {   // distribution
01085     if (calc->isZero (lambda))
01086       result = 1;
01087     else
01088     {
01089       result = 1.0;
01090       Value fFak = 1.0;
01091       unsigned long nEnd = calc->conv()->asInteger (x).asInteger();
01092       for (unsigned long i = 1; i <= nEnd; i++)
01093       {
01094         // fFak *= i
01095         fFak = calc->mul (fFak, i);
01096         // result += pow (lambda, i) / fFak
01097         result = calc->add (result, calc->div (calc->pow (lambda, i), fFak));
01098       }
01099       result = calc->mul (result, ex);
01100     }
01101   }
01102 
01103   return result;
01104 }
01105 
01106 // Function: confidence
01107 Value func_confidence (valVector args, ValueCalc *calc, FuncExtra *) {
01108   //returns the confidence interval for a population mean
01109   Value alpha = args[0];
01110   Value sigma = args[1];
01111   Value n = args[2];
01112 
01113   // sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1
01114   if ((!calc->greater (sigma, 0.0)) || (!calc->greater (alpha, 0.0)) ||
01115       (!calc->lower (alpha, 1.0)) || calc->lower (n, 1))
01116     return Value::errorVALUE();
01117     
01118   // g = gaussinv (1.0 - alpha / 2.0)
01119   Value g = calc->gaussinv (calc->sub (1.0, calc->div (alpha, 2.0)));
01120   // g * sigma / sqrt (n)
01121   return calc->div (calc->mul (g, sigma), calc->sqrt (n));
01122 }
01123 
01124 // Function: tdist
01125 Value func_tdist (valVector args, ValueCalc *calc, FuncExtra *) {
01126   //returns the t-distribution
01127   
01128   Value T = args[0];
01129   Value fDF = args[1];
01130   int flag = calc->conv()->asInteger (args[2]).asInteger();
01131 
01132   if (calc->lower (fDF, 1) || calc->lower (T, 0.0) || (flag != 1 && flag != 2))
01133     return Value::errorVALUE();
01134 
01135   // arg = fDF / (fDF + T * T)
01136   Value arg = calc->div (fDF, calc->add (fDF, calc->sqr (T)));
01137   
01138   Value R;
01139   R = calc->mul (calc->GetBeta (arg, calc->div (fDF, 2.0), 0.5), 0.5);
01140   
01141   if (flag == 1)
01142     return R;
01143   return calc->mul (R, 2);   // flag is 2 here
01144 }
01145 
01146 // Function: fdist
01147 Value func_fdist (valVector args, ValueCalc *calc, FuncExtra *) {
01148   //returns the f-distribution
01149 
01150   Value x = args[0];
01151   Value fF1 = args[1];
01152   Value fF2 = args[2];
01153 
01154   // x < 0.0 || fF1 < 1 || fF2 < 1 || fF1 >= 1.0E10 || fF2 >= 1.0E10
01155   if (calc->lower (x, 0.0) || calc->lower (fF1, 1) || calc->lower (fF2, 1) ||
01156       (!calc->lower (fF1, 1.0E10)) || (!calc->lower (fF2, 1.0E10)))
01157     return Value::errorVALUE();
01158 
01159   // arg = fF2 / (fF2 + fF1 * x)
01160   Value arg = calc->div (fF2, calc->add (fF2, calc->mul (fF1, x)));
01161   // alpha = fF2/2.0
01162   Value alpha = calc->div (fF2, 2.0);
01163   // beta = fF1/2.0
01164   Value beta = calc->div (fF1, 2.0);
01165   return calc->GetBeta (arg, alpha, beta);
01166 }
01167 
01168 // Function: chidist
01169 Value func_chidist (valVector args, ValueCalc *calc, FuncExtra *) {
01170   //returns the chi-distribution
01171   
01172   Value fChi = args[0];
01173   Value fDF = args[1];
01174 
01175   // fDF < 1 || fDF >= 1.0E5 || fChi < 0.0
01176   if (calc->lower (fDF, 1) || (!calc->lower (fDF, 1.0E5)) ||
01177       calc->lower (fChi, 0.0))
01178     return Value::errorVALUE();
01179 
01180   // 1.0 - GetGammaDist (fChi / 2.0, fDF / 2.0, 1.0)
01181   return calc->sub (1.0, calc->GetGammaDist (calc->div (fChi, 2.0),
01182       calc->div (fDF, 2.0), 1.0));
01183 }
01184 
01185 
01186 // two-array-walk functions used in the two-sum functions
01187 
01188 void tawSumproduct (ValueCalc *c, Value &res, Value v1,
01189     Value v2) {
01190   // res += v1*v2
01191   res = c->add (res, c->mul (v1, v2));  
01192 }
01193 
01194 void tawSumx2py2 (ValueCalc *c, Value &res, Value v1,
01195     Value v2) {
01196   // res += sqr(v1)+sqr(v2)
01197   res = c->add (res, c->add (c->sqr (v1), c->sqr (v2)));
01198 }
01199 
01200 void tawSumx2my2 (ValueCalc *c, Value &res, Value v1,
01201     Value v2) {
01202   // res += sqr(v1)-sqr(v2)
01203   res = c->add (res, c->sub (c->sqr (v1), c->sqr (v2)));
01204 }
01205 
01206 void tawSumxmy2 (ValueCalc *c, Value &res, Value v1,
01207     Value v2) {
01208   // res += sqr(v1-v2)
01209   res = c->add (res, c->sqr (c->sub (v1, v2)));
01210   
01211 }
01212 
01213 // Function: sumproduct
01214 Value func_sumproduct (valVector args, ValueCalc *calc, FuncExtra *)
01215 {
01216   Value result;
01217   calc->twoArrayWalk (args[0], args[1], result, tawSumproduct);
01218   return result;
01219 }
01220 
01221 // Function: sumx2py2
01222 Value func_sumx2py2 (valVector args, ValueCalc *calc, FuncExtra *)
01223 {
01224   Value result;
01225   calc->twoArrayWalk (args[0], args[1], result, tawSumx2py2);
01226   return result;
01227 }
01228 
01229 // Function: sumx2my2
01230 Value func_sumx2my2 (valVector args, ValueCalc *calc, FuncExtra *)
01231 {
01232   Value result;
01233   calc->twoArrayWalk (args[0], args[1], result, tawSumx2my2);
01234   return result;
01235 }
01236 
01237 // Function: sum2xmy
01238 Value func_sumxmy2 (valVector args, ValueCalc *calc, FuncExtra *)
01239 {
01240   Value result;
01241   calc->twoArrayWalk (args[0], args[1], result, tawSumxmy2);
01242   return result;
01243 }
KDE Home | KDE Accessibility Home | Description of Access Keys