File indexing completed on 2024-05-12 16:35:35
0001 /* This file is part of the KDE project 0002 Copyright (C) 1998-2002 The KSpread Team <calligra-devel@kde.org> 0003 Copyright (C) 2005 Tomas Mecir <mecirt@gmail.com> 0004 Copyright (C) 2007 Sascha Pfau <MrPeacock@gmail.com> 0005 0006 This library is free software; you can redistribute it and/or 0007 modify it under the terms of the GNU Library General Public 0008 License as published by the Free Software Foundation; only 0009 version 2 of the License. 0010 0011 This library is distributed in the hope that it will be useful, 0012 but WITHOUT ANY WARRANTY; without even the implied warranty of 0013 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 0014 Library General Public License for more details. 0015 0016 You should have received a copy of the GNU Library General Public License 0017 along with this library; see the file COPYING.LIB. If not, write to 0018 the Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, 0019 Boston, MA 02110-1301, USA. 0020 */ 0021 0022 // built-in statistical functions 0023 #include "StatisticalModule.h" 0024 0025 #include "Function.h" 0026 #include "FunctionModuleRegistry.h" 0027 #include "ValueCalc.h" 0028 #include "ValueConverter.h" 0029 #include "SheetsDebug.h" 0030 0031 #include <Formula.h> 0032 0033 // needed for MODE 0034 #include <QList> 0035 #include <QMap> 0036 0037 #include <algorithm> 0038 0039 using namespace Calligra::Sheets; 0040 0041 // prototypes (sorted!) 0042 Value func_arrang(valVector args, ValueCalc *calc, FuncExtra *); 0043 Value func_average(valVector args, ValueCalc *calc, FuncExtra *); 0044 Value func_averagea(valVector args, ValueCalc *calc, FuncExtra *); 0045 Value func_averageif(valVector args, ValueCalc *calc, FuncExtra *); 0046 Value func_averageifs(valVector args, ValueCalc *calc, FuncExtra *e); 0047 Value func_avedev(valVector args, ValueCalc *calc, FuncExtra *); 0048 Value func_betadist(valVector args, ValueCalc *calc, FuncExtra *); 0049 Value func_betainv(valVector args, ValueCalc *calc, FuncExtra *); 0050 Value func_bino(valVector args, ValueCalc *calc, FuncExtra *); 0051 Value func_binomdist(valVector args, ValueCalc *calc, FuncExtra *); 0052 Value func_chidist(valVector args, ValueCalc *calc, FuncExtra *); 0053 Value func_combin(valVector args, ValueCalc *calc, FuncExtra *); 0054 Value func_combina(valVector args, ValueCalc *calc, FuncExtra *); 0055 Value func_confidence(valVector args, ValueCalc *calc, FuncExtra *); 0056 Value func_correl_pop(valVector args, ValueCalc *calc, FuncExtra *); 0057 Value func_covar(valVector args, ValueCalc *calc, FuncExtra *); 0058 Value func_devsq(valVector args, ValueCalc *calc, FuncExtra *); 0059 Value func_devsqa(valVector args, ValueCalc *calc, FuncExtra *); 0060 Value func_expondist(valVector args, ValueCalc *calc, FuncExtra *); 0061 Value func_fdist(valVector args, ValueCalc *calc, FuncExtra *); 0062 Value func_finv(valVector args, ValueCalc *calc, FuncExtra *); 0063 Value func_fisher(valVector args, ValueCalc *calc, FuncExtra *); 0064 Value func_fisherinv(valVector args, ValueCalc *calc, FuncExtra *); 0065 Value func_frequency(valVector args, ValueCalc *calc, FuncExtra *); 0066 Value func_ftest(valVector args, ValueCalc *calc, FuncExtra *); 0067 Value func_gammadist(valVector args, ValueCalc *calc, FuncExtra *); 0068 Value func_gammainv(valVector args, ValueCalc *calc, FuncExtra *); 0069 Value func_gammaln(valVector args, ValueCalc *calc, FuncExtra *); 0070 Value func_gauss(valVector args, ValueCalc *calc, FuncExtra *); 0071 Value func_growth(valVector args, ValueCalc *calc, FuncExtra *); 0072 Value func_geomean(valVector args, ValueCalc *calc, FuncExtra *); 0073 Value func_harmean(valVector args, ValueCalc *calc, FuncExtra *); 0074 Value func_hypgeomdist(valVector args, ValueCalc *calc, FuncExtra *); 0075 Value func_intercept(valVector args, ValueCalc *calc, FuncExtra *); 0076 Value func_kurtosis_est(valVector args, ValueCalc *calc, FuncExtra *); 0077 Value func_kurtosis_pop(valVector args, ValueCalc *calc, FuncExtra *); 0078 Value func_large(valVector args, ValueCalc *calc, FuncExtra *); 0079 Value func_legacychidist(valVector args, ValueCalc *calc, FuncExtra *); 0080 Value func_legacychiinv(valVector args, ValueCalc *calc, FuncExtra *); 0081 Value func_legacyfdist(valVector args, ValueCalc *calc, FuncExtra *); 0082 Value func_legacyfinv(valVector args, ValueCalc *calc, FuncExtra *); 0083 Value func_loginv(valVector args, ValueCalc *calc, FuncExtra *); 0084 Value func_lognormdist(valVector args, ValueCalc *calc, FuncExtra *); 0085 Value func_median(valVector args, ValueCalc *calc, FuncExtra *); 0086 Value func_mode(valVector args, ValueCalc *calc, FuncExtra *); 0087 Value func_negbinomdist(valVector args, ValueCalc *calc, FuncExtra *); 0088 Value func_normdist(valVector args, ValueCalc *calc, FuncExtra *); 0089 Value func_norminv(valVector args, ValueCalc *calc, FuncExtra *); 0090 Value func_normsinv(valVector args, ValueCalc *calc, FuncExtra *); 0091 Value func_percentile(valVector args, ValueCalc *calc, FuncExtra *); 0092 Value func_permutationa(valVector args, ValueCalc *calc, FuncExtra *); 0093 Value func_phi(valVector args, ValueCalc *calc, FuncExtra *); 0094 Value func_poisson(valVector args, ValueCalc *calc, FuncExtra *); 0095 Value func_rank(valVector args, ValueCalc *calc, FuncExtra *); 0096 Value func_rsq(valVector args, ValueCalc *calc, FuncExtra *); 0097 Value func_quartile(valVector args, ValueCalc *calc, FuncExtra *); 0098 Value func_skew_est(valVector args, ValueCalc *calc, FuncExtra *); 0099 Value func_skew_pop(valVector args, ValueCalc *calc, FuncExtra *); 0100 Value func_slope(valVector args, ValueCalc *calc, FuncExtra *); 0101 Value func_small(valVector args, ValueCalc *calc, FuncExtra *); 0102 Value func_standardize(valVector args, ValueCalc *calc, FuncExtra *); 0103 Value func_stddev(valVector args, ValueCalc *calc, FuncExtra *); 0104 Value func_stddeva(valVector args, ValueCalc *calc, FuncExtra *); 0105 Value func_stddevp(valVector args, ValueCalc *calc, FuncExtra *); 0106 Value func_stddevpa(valVector args, ValueCalc *calc, FuncExtra *); 0107 Value func_stdnormdist(valVector args, ValueCalc *calc, FuncExtra *); 0108 Value func_steyx(valVector args, ValueCalc *calc, FuncExtra *); 0109 Value func_sumproduct(valVector args, ValueCalc *calc, FuncExtra *); 0110 Value func_sumx2py2(valVector args, ValueCalc *calc, FuncExtra *); 0111 Value func_sumx2my2(valVector args, ValueCalc *calc, FuncExtra *); 0112 Value func_sumxmy2(valVector args, ValueCalc *calc, FuncExtra *); 0113 Value func_tdist(valVector args, ValueCalc *calc, FuncExtra *); 0114 Value func_tinv(valVector args, ValueCalc *calc, FuncExtra *); 0115 Value func_trend(valVector args, ValueCalc *calc, FuncExtra *); 0116 Value func_trimmean(valVector args, ValueCalc *calc, FuncExtra *); 0117 Value func_ttest(valVector args, ValueCalc *calc, FuncExtra *); 0118 Value func_variance(valVector args, ValueCalc *calc, FuncExtra *); 0119 Value func_variancea(valVector args, ValueCalc *calc, FuncExtra *); 0120 Value func_variancep(valVector args, ValueCalc *calc, FuncExtra *); 0121 Value func_variancepa(valVector args, ValueCalc *calc, FuncExtra *); 0122 Value func_weibull(valVector args, ValueCalc *calc, FuncExtra *); 0123 Value func_ztest(valVector args, ValueCalc *calc, FuncExtra *); 0124 0125 typedef QList<double> List; 0126 0127 0128 CALLIGRA_SHEETS_EXPORT_FUNCTION_MODULE("kspreadstatisticalmodule.json", StatisticalModule) 0129 0130 0131 StatisticalModule::StatisticalModule(QObject* parent, const QVariantList&) 0132 : FunctionModule(parent) 0133 { 0134 Function *f; 0135 0136 f = new Function("AVEDEV", func_avedev); 0137 f->setParamCount(1, -1); 0138 f->setAcceptArray(); 0139 add(f); 0140 f = new Function("AVERAGE", func_average); 0141 f->setParamCount(1, -1); 0142 f->setAcceptArray(); 0143 add(f); 0144 f = new Function("AVERAGEA", func_averagea); 0145 f->setParamCount(1, -1); 0146 f->setAcceptArray(); 0147 add(f); 0148 f = new Function("AVERAGEIF", func_averageif); 0149 f->setParamCount(2, 3); 0150 f->setAcceptArray(); 0151 f->setNeedsExtra(true); 0152 add(f); 0153 f = new Function("AVERAGEIFS", func_averageifs); 0154 f->setParamCount(3, -1); 0155 f->setAcceptArray(); 0156 f->setNeedsExtra(true); 0157 add(f); 0158 f = new Function("BETADIST", func_betadist); 0159 f->setParamCount(3, 6); 0160 add(f); 0161 f = new Function("BETAINV", func_betainv); 0162 f->setParamCount(3, 5); 0163 add(f); 0164 f = new Function("BINO", func_bino); 0165 f->setParamCount(3); 0166 add(f); 0167 f = new Function("BINOMDIST", func_binomdist); 0168 f->setParamCount(4); 0169 add(f); 0170 f = new Function("CHIDIST", func_chidist); 0171 f->setParamCount(2); 0172 add(f); 0173 f = new Function("COMBIN", func_combin); 0174 f->setParamCount(2); 0175 add(f); 0176 f = new Function("COMBINA", func_combina); 0177 f->setParamCount(2); 0178 add(f); 0179 f = new Function("CONFIDENCE", func_confidence); 0180 f->setParamCount(3); 0181 add(f); 0182 f = new Function("CORREL", func_correl_pop); 0183 f->setParamCount(2); 0184 f->setAcceptArray(); 0185 add(f); 0186 f = new Function("COVAR", func_covar); 0187 f->setParamCount(2); 0188 f->setAcceptArray(); 0189 add(f); 0190 f = new Function("DEVSQ", func_devsq); 0191 f->setParamCount(1, -1); 0192 f->setAcceptArray(); 0193 add(f); 0194 f = new Function("DEVSQA", func_devsqa); 0195 f->setParamCount(1, -1); 0196 f->setAcceptArray(); 0197 add(f); 0198 f = new Function("EXPONDIST", func_expondist); 0199 f->setParamCount(3); 0200 add(f); 0201 f = new Function("FDIST", func_fdist); 0202 f->setParamCount(3, 4); 0203 add(f); 0204 f = new Function("FINV", func_finv); 0205 f->setParamCount(3); 0206 add(f); 0207 f = new Function("FISHER", func_fisher); 0208 add(f); 0209 f = new Function("FISHERINV", func_fisherinv); 0210 add(f); 0211 f = new Function("FREQUENCY", func_frequency); 0212 f->setParamCount(2); 0213 f->setAcceptArray(); 0214 add(f); 0215 f = new Function("FTEST", func_ftest); 0216 f->setParamCount(2); 0217 f->setAcceptArray(); 0218 add(f); 0219 f = new Function("GAMMADIST", func_gammadist); 0220 f->setParamCount(4); 0221 add(f); 0222 f = new Function("GAMMAINV", func_gammainv); 0223 f->setParamCount(3); 0224 add(f); 0225 f = new Function("GAMMALN", func_gammaln); 0226 add(f); 0227 f = new Function("GAUSS", func_gauss); 0228 add(f); 0229 f = new Function("GROWTH", func_growth); 0230 f->setParamCount(1, 4); 0231 f->setAcceptArray(); 0232 add(f); 0233 f = new Function("GEOMEAN", func_geomean); 0234 f->setParamCount(1, -1); 0235 f->setAcceptArray(); 0236 add(f); 0237 f = new Function("HARMEAN", func_harmean); 0238 f->setParamCount(1, -1); 0239 f->setAcceptArray(); 0240 add(f); 0241 f = new Function("HYPGEOMDIST", func_hypgeomdist); 0242 f->setParamCount(4, 5); 0243 add(f); 0244 f = new Function("INTERCEPT", func_intercept); 0245 f->setParamCount(2); 0246 f->setAcceptArray(); 0247 add(f); 0248 f = new Function("INVBINO", func_bino); // same as BINO, for 1.4 compat 0249 add(f); 0250 f = new Function("KURT", func_kurtosis_est); 0251 f->setParamCount(1, -1); 0252 f->setAcceptArray(); 0253 add(f); 0254 f = new Function("KURTP", func_kurtosis_pop); 0255 f->setParamCount(1, -1); 0256 f->setAcceptArray(); 0257 add(f); 0258 f = new Function("LARGE", func_large); 0259 f->setParamCount(2); 0260 f->setAcceptArray(); 0261 add(f); 0262 f = new Function("LEGACYCHIDIST", func_legacychidist); 0263 f->setParamCount(2); 0264 add(f); 0265 f = new Function("LEGACYCHIINV", func_legacychiinv); 0266 f->setParamCount(2); 0267 add(f); 0268 f = new Function("LEGACYFDIST", func_legacyfdist); 0269 f->setParamCount(3); 0270 add(f); 0271 f = new Function("LEGACYFINV", func_legacyfinv); 0272 f->setParamCount(3); 0273 add(f); 0274 // this is meant to be a copy of the function NORMSDIST. 0275 // required for OpenFormula compliance. 0276 f = new Function("LEGACYNORMSDIST", func_stdnormdist); 0277 add(f); 0278 // this is meant to be a copy of the function NORMSINV. 0279 // required for OpenFormula compliance. 0280 f = new Function("LEGACYNORMSINV", func_normsinv); 0281 add(f); 0282 f = new Function("LOGINV", func_loginv); 0283 f->setParamCount(1, 3); 0284 add(f); 0285 f = new Function("LOGNORMDIST", func_lognormdist); 0286 f->setParamCount(1, 4); 0287 add(f); 0288 f = new Function("MEDIAN", func_median); 0289 f->setParamCount(1, -1); 0290 f->setAcceptArray(); 0291 add(f); 0292 f = new Function("MODE", func_mode); 0293 f->setParamCount(1, -1); 0294 f->setAcceptArray(); 0295 add(f); 0296 f = new Function("NEGBINOMDIST", func_negbinomdist); 0297 f->setParamCount(3); 0298 add(f); 0299 f = new Function("NORMDIST", func_normdist); 0300 f->setParamCount(4); 0301 add(f); 0302 f = new Function("NORMINV", func_norminv); 0303 f->setParamCount(3); 0304 add(f); 0305 f = new Function("NORMSDIST", func_stdnormdist); 0306 add(f); 0307 f = new Function("NORMSINV", func_normsinv); 0308 add(f); 0309 f = new Function("PEARSON", func_correl_pop); 0310 f->setParamCount(2); 0311 f->setAcceptArray(); 0312 add(f); 0313 f = new Function("PERCENTILE", func_percentile); 0314 f->setParamCount(2); 0315 f->setAcceptArray(); 0316 add(f); 0317 f = new Function("PERMUT", func_arrang); 0318 f->setParamCount(2); 0319 add(f); 0320 f = new Function("PERMUTATIONA", func_permutationa); 0321 f->setParamCount(2); 0322 add(f); 0323 f = new Function("PHI", func_phi); 0324 add(f); 0325 f = new Function("POISSON", func_poisson); 0326 f->setParamCount(3); 0327 add(f); 0328 f = new Function("RANK", func_rank); 0329 f->setParamCount(2, 3); 0330 f->setAcceptArray(); 0331 add(f); 0332 f = new Function("RSQ", func_rsq); 0333 f->setParamCount(2); 0334 f->setAcceptArray(); 0335 add(f); 0336 f = new Function("QUARTILE", func_quartile); 0337 f->setParamCount(2); 0338 f->setAcceptArray(); 0339 add(f); 0340 f = new Function("SKEW", func_skew_est); 0341 f->setParamCount(1, -1); 0342 f->setAcceptArray(); 0343 add(f); 0344 f = new Function("SKEWP", func_skew_pop); 0345 f->setParamCount(1, -1); 0346 f->setAcceptArray(); 0347 add(f); 0348 f = new Function("SLOPE", func_slope); 0349 f->setParamCount(2); 0350 f->setAcceptArray(); 0351 add(f); 0352 f = new Function("SMALL", func_small); 0353 f->setParamCount(2); 0354 f->setAcceptArray(); 0355 add(f); 0356 f = new Function("STANDARDIZE", func_standardize); 0357 f->setParamCount(3); 0358 add(f); 0359 f = new Function("STDEV", func_stddev); 0360 f->setParamCount(1, -1); 0361 f->setAcceptArray(); 0362 add(f); 0363 f = new Function("STDEVA", func_stddeva); 0364 f->setParamCount(1, -1); 0365 f->setAcceptArray(); 0366 add(f); 0367 f = new Function("STDEVP", func_stddevp); 0368 f->setParamCount(1, -1); 0369 f->setAcceptArray(); 0370 add(f); 0371 f = new Function("STDEVPA", func_stddevpa); 0372 f->setParamCount(1, -1); 0373 f->setAcceptArray(); 0374 add(f); 0375 f = new Function("STEYX", func_steyx); 0376 f->setParamCount(2); 0377 f->setAcceptArray(); 0378 add(f); 0379 f = new Function("SUM2XMY", func_sumxmy2); // deprecated, use SUMXMY2 0380 f->setParamCount(2); 0381 f->setAcceptArray(); 0382 add(f); 0383 f = new Function("SUMXMY2", func_sumxmy2); 0384 f->setParamCount(2); 0385 f->setAcceptArray(); 0386 add(f); 0387 f = new Function("SUMPRODUCT", func_sumproduct); 0388 f->setParamCount(2); 0389 f->setAcceptArray(); 0390 add(f); 0391 f = new Function("SUMX2PY2", func_sumx2py2); 0392 f->setParamCount(2); 0393 f->setAcceptArray(); 0394 add(f); 0395 f = new Function("SUMX2MY2", func_sumx2my2); 0396 f->setParamCount(2); 0397 f->setAcceptArray(); 0398 add(f); 0399 f = new Function("TDIST", func_tdist); 0400 f->setParamCount(3); 0401 add(f); 0402 f = new Function("TINV", func_tinv); 0403 f->setParamCount(2); 0404 add(f); 0405 f = new Function("TREND", func_trend); 0406 f->setParamCount(1, 4); 0407 f->setAcceptArray(); 0408 add(f); 0409 f = new Function("TRIMMEAN", func_trimmean); 0410 f->setParamCount(2); 0411 f->setAcceptArray(); 0412 add(f); 0413 f = new Function("TTEST", func_ttest); 0414 f->setParamCount(4); 0415 f->setAcceptArray(); 0416 add(f); 0417 f = new Function("VARIANCE", func_variance); 0418 f->setParamCount(1, -1); 0419 f->setAcceptArray(); 0420 add(f); 0421 f = new Function("VAR", func_variance); 0422 f->setParamCount(1, -1); 0423 f->setAcceptArray(); 0424 add(f); 0425 f = new Function("VARP", func_variancep); 0426 f->setParamCount(1, -1); 0427 f->setAcceptArray(); 0428 add(f); 0429 f = new Function("VARA", func_variancea); 0430 f->setParamCount(1, -1); 0431 f->setAcceptArray(); 0432 add(f); 0433 f = new Function("VARPA", func_variancepa); 0434 f->setParamCount(1, -1); 0435 f->setAcceptArray(); 0436 add(f); 0437 f = new Function("WEIBULL", func_weibull); 0438 f->setParamCount(4); 0439 add(f); 0440 f = new Function("ZTEST", func_ztest); 0441 f->setParamCount(2, 3); 0442 f->setAcceptArray(); 0443 add(f); 0444 } 0445 0446 QString StatisticalModule::descriptionFileName() const 0447 { 0448 return QString("statistical.xml"); 0449 } 0450 0451 0452 /////////////////////////////////////////////////////////// 0453 // 0454 // helper functions 0455 // 0456 /////////////////////////////////////////////////////////// 0457 0458 // 0459 // helper: array_helper 0460 // 0461 void func_array_helper(Value range, ValueCalc *calc, 0462 List &array, int &number) 0463 { 0464 if (!range.isArray()) { 0465 array << numToDouble(calc->conv()->toFloat(range)); 0466 ++number; 0467 return; 0468 } 0469 0470 for (unsigned int row = 0; row < range.rows(); ++row) 0471 for (unsigned int col = 0; col < range.columns(); ++col) { 0472 Value v = range.element(col, row); 0473 if (v.isArray()) 0474 func_array_helper(v, calc, array, number); 0475 else { 0476 array << numToDouble(calc->conv()->toFloat(v)); 0477 ++number; 0478 } 0479 } 0480 } 0481 0482 0483 // 0484 // helper: covar_helper 0485 // 0486 Value func_covar_helper(Value range1, Value range2, 0487 ValueCalc *calc, Value avg1, Value avg2) 0488 { 0489 // two arrays -> cannot use arrayWalk 0490 if ((!range1.isArray()) && (!range2.isArray())) 0491 // (v1-E1)*(v2-E2) 0492 return calc->mul(calc->sub(range1, avg1), calc->sub(range2, avg2)); 0493 0494 int rows = range1.rows(); 0495 int cols = range1.columns(); 0496 int rows2 = range2.rows(); 0497 int cols2 = range2.columns(); 0498 if ((rows != rows2) || (cols != cols2)) 0499 return Value::errorVALUE(); 0500 0501 Value result(0.0); 0502 for (int row = 0; row < rows; ++row) 0503 for (int col = 0; col < cols; ++col) { 0504 Value v1 = range1.element(col, row); 0505 Value v2 = range2.element(col, row); 0506 if (v1.isArray() || v2.isArray()) 0507 result = calc->add(result, 0508 func_covar_helper(v1, v2, calc, avg1, avg2)); 0509 else 0510 // result += (v1-E1)*(v2-E2) 0511 result = calc->add(result, calc->mul(calc->sub(v1, avg1), 0512 calc->sub(v2, avg2))); 0513 } 0514 0515 return result; 0516 } 0517 0518 0519 /* 0520 // helper: GetValue - Returns result of a formula. 0521 static double GetValue(const QString& formula, const double x) 0522 { 0523 Formula f; 0524 QString expr = formula; 0525 if (expr[0] != '=') 0526 expr.prepend('='); 0527 0528 expr.replace(QString("x"), QString::number(x, 'g', 12)); 0529 0530 //debugSheets<<"expression"<<expr; 0531 f.setExpression(expr); 0532 Value result = f.eval(); 0533 0534 return result.asFloat(); 0535 } 0536 */ 0537 0538 /** 0539 * Helper-class to inverse iterate over a value to determinate 0540 * an unknown value. 0541 */ 0542 class InverseIterator 0543 { 0544 public: 0545 InverseIterator(FunctionPtr ptr, const valVector &args, ValueCalc *calc) 0546 : m_caller(FunctionCaller(ptr, args, calc)) {} 0547 Value exec(double unknown, double x0, double x1, bool& convergenceError); 0548 private: 0549 FunctionCaller m_caller; 0550 inline double getValue(Value arg) { 0551 valVector args = m_caller.m_args; 0552 args.prepend(arg); 0553 return m_caller.exec(args).asFloat(); 0554 } 0555 inline double getValue(double arg) { 0556 return getValue(Value(arg)); 0557 } 0558 }; 0559 0560 Value InverseIterator::exec(double unknown, double x0, double x1, bool& convergenceError) 0561 // static Value InverseIterator(const double unknown, FunctionPtr ptr, double x0, double x1, bool& convergenceError) 0562 { 0563 convergenceError = false; // reset error flag 0564 double eps = 1.0E-7; // define Epsilon 0565 0566 debugSheets << "searching for " << unknown << " in interval x0=" << x0 << " x1=" << x1; 0567 0568 if (x0 > x1) 0569 debugSheets << "InverseIterator: wrong interval"; 0570 0571 double f0 = unknown - getValue(x0); 0572 double f1 = unknown - getValue(x1); 0573 0574 debugSheets << " f(" << x0 << ") =" << f0; 0575 debugSheets << " f(" << x1 << ") =" << f1; 0576 0577 double xs; 0578 int i; 0579 for (i = 0; i < 1000 && f0*f1 > 0.0; ++i) { 0580 if (fabs(f0) <= fabs(f1)) { 0581 xs = x0; 0582 x0 += 2.0 * (x0 - x1); 0583 if (x0 < 0.0) 0584 x0 = 0.0; 0585 x1 = xs; 0586 f1 = f0; 0587 f0 = unknown - getValue(x0); 0588 } else { 0589 xs = x1; 0590 x1 += 2.0 * (x1 - x0); 0591 x0 = xs; 0592 f0 = f1; 0593 f1 = unknown - getValue(x1); 0594 } 0595 } 0596 0597 if (f0 == 0.0) 0598 return Value(x0); 0599 if (f1 == 0.0) 0600 return Value(x1); 0601 0602 // simple iteration 0603 //debugSheets<<"simple iteration f0="<<f0<<" f1="<<f1; 0604 0605 double x00 = x0; 0606 double x11 = x1; 0607 double fs = 0.0; 0608 for (i = 0; i < 100; ++i) { 0609 xs = 0.5 * (x0 + x1); 0610 if (fabs(f1 - f0) >= eps) { 0611 fs = unknown - getValue(xs); 0612 if (f0*fs <= 0.0) { 0613 x1 = xs; 0614 f1 = fs; 0615 } else { 0616 x0 = xs; 0617 f0 = fs; 0618 } 0619 } else { 0620 //add one step of regula falsi to improve precision 0621 if (x0 != x1) { 0622 double regxs = (f1 - f0) / (x1 - x0); 0623 if (regxs != 0.0) { 0624 double regx = x1 - f1 / regxs; 0625 if (regx >= x00 && regx <= x11) { 0626 double regfs = unknown - getValue(regx); 0627 if (fabs(regfs) < fabs(fs)) 0628 xs = regx; 0629 } 0630 } 0631 } 0632 return Value(xs); 0633 } 0634 // debugSheets<<"probe no. "<<i<<" : "<<xs<<" error diff ="<<fs; 0635 } 0636 0637 // error no convergence - set flag 0638 convergenceError = true; 0639 return Value(0.0); 0640 } 0641 0642 // 0643 // helper: mode_helper 0644 // 0645 class ContentSheet : public QMap<double, int> {}; 0646 0647 void func_mode_helper(Value range, ValueCalc *calc, ContentSheet &sh) 0648 { 0649 if (!range.isArray()) { 0650 double d = numToDouble(calc->conv()->toFloat(range)); 0651 ++sh[d]; 0652 return; 0653 } 0654 0655 for (unsigned int row = 0; row < range.rows(); ++row) 0656 for (unsigned int col = 0; col < range.columns(); ++col) { 0657 Value v = range.element(col, row); 0658 if (v.isArray()) 0659 func_mode_helper(v, calc, sh); 0660 else { 0661 double d = numToDouble(calc->conv()->toFloat(v)); 0662 ++sh[d]; 0663 } 0664 } 0665 } 0666 0667 /////////////////////////////////////////////////////////// 0668 // 0669 // array-walk functions used in this file 0670 // 0671 /////////////////////////////////////////////////////////// 0672 0673 void awSkew(ValueCalc *c, Value &res, Value val, Value p) 0674 { 0675 Value avg = p.element(0, 0); 0676 Value stdev = p.element(1, 0); 0677 // (val - avg) / stddev 0678 Value d = c->div(c->sub(val, avg), stdev); 0679 // res += d*d*d 0680 res = c->add(res, c->mul(d, c->mul(d, d))); 0681 } 0682 0683 void awSumInv(ValueCalc *c, Value &res, Value val, Value) 0684 { 0685 // res += 1/value 0686 res = c->add(res, c->div(Value(1.0), val)); 0687 } 0688 0689 void awAveDev(ValueCalc *c, Value &res, Value val, 0690 Value p) 0691 { 0692 // res += abs (val - p) 0693 res = c->add(res, c->abs(c->sub(val, p))); 0694 } 0695 0696 void awKurtosis(ValueCalc *c, Value &res, Value val, 0697 Value p) 0698 { 0699 Value avg = p.element(0, 0); 0700 Value stdev = p.element(1, 0); 0701 //d = (val - avg) / stdev 0702 Value d = c->div(c->sub(val, avg), stdev); 0703 // res += d^4 0704 res = c->add(res, c->pow(d, 4)); 0705 } 0706 0707 /////////////////////////////////////////////////////////// 0708 // 0709 // two-array-walk functions used in the two-sum functions 0710 // 0711 /////////////////////////////////////////////////////////// 0712 0713 void tawSumproduct(ValueCalc *c, Value &res, Value v1, 0714 Value v2) 0715 { 0716 // res += v1*v2 0717 res = c->add(res, c->mul(v1, v2)); 0718 } 0719 0720 void tawSumx2py2(ValueCalc *c, Value &res, Value v1, 0721 Value v2) 0722 { 0723 // res += sqr(v1)+sqr(v2) 0724 res = c->add(res, c->add(c->sqr(v1), c->sqr(v2))); 0725 } 0726 0727 void tawSumx2my2(ValueCalc *c, Value &res, Value v1, 0728 Value v2) 0729 { 0730 // res += sqr(v1)-sqr(v2) 0731 res = c->add(res, c->sub(c->sqr(v1), c->sqr(v2))); 0732 } 0733 0734 void tawSumxmy2(ValueCalc *c, Value &res, Value v1, 0735 Value v2) 0736 { 0737 // res += sqr(v1-v2) 0738 res = c->add(res, c->sqr(c->sub(v1, v2))); 0739 } 0740 0741 void tawSumxmy(ValueCalc *c, Value &res, Value v1, 0742 Value v2) 0743 { 0744 // res += (v1-v2) 0745 res = c->add(res, c->sub(v1, v2)); 0746 } 0747 0748 /////////////////////////////////////////////////////////// 0749 // 0750 // functions used in this file 0751 // 0752 /////////////////////////////////////////////////////////// 0753 0754 // 0755 // Function: permut 0756 // 0757 Value func_arrang(valVector args, ValueCalc *calc, FuncExtra *) 0758 { 0759 Value n = args[0]; 0760 Value m = args[1]; 0761 if (calc->lower(n, m)) // problem if n<m 0762 return Value::errorVALUE(); 0763 0764 if (calc->lower(m, Value(0))) // problem if m<0 (n>=m so that's okay) 0765 return Value::errorVALUE(); 0766 0767 // fact(n) / (fact(n-m) 0768 return calc->fact(n, calc->sub(n, m)); 0769 } 0770 0771 // 0772 // Function: average 0773 // 0774 Value func_average(valVector args, ValueCalc *calc, FuncExtra *) 0775 { 0776 return calc->avg(args, false); 0777 } 0778 0779 // 0780 // Function: averagea 0781 // 0782 Value func_averagea(valVector args, ValueCalc *calc, FuncExtra *) 0783 { 0784 return calc->avg(args); 0785 } 0786 0787 // 0788 // Function: averageif 0789 // 0790 Value func_averageif(valVector args, ValueCalc *calc, FuncExtra *e) 0791 { 0792 Value checkRange = args[0]; 0793 QString condition = calc->conv()->asString(args[1]).asString(); 0794 Condition cond; 0795 calc->getCond(cond, Value(condition)); 0796 0797 if (args.count() == 3) { 0798 Cell avgRangeStart(e->sheet, e->ranges[2].col1, e->ranges[2].row1); 0799 return calc->averageIf(avgRangeStart, checkRange, cond); 0800 } else { 0801 return calc->averageIf(checkRange, cond); 0802 } 0803 } 0804 0805 // 0806 // Function: averageifs 0807 // 0808 Value func_averageifs(valVector args, ValueCalc *calc, FuncExtra *e) 0809 { 0810 int lim = (int) (args.count()-1)/2; 0811 0812 QList<Value> c_Range; 0813 QStringList condition; 0814 QList<Condition> cond; 0815 0816 c_Range.append(args.value(0)); //first element - range to be operated on 0817 0818 for (int i = 1; i < args.count(); i += 2) { 0819 c_Range.append(args[i]); 0820 condition.append(calc->conv()->asString(args[i+1]).asString()); 0821 Condition c; 0822 calc->getCond(c, Value(condition.last())); 0823 cond.append(c); 0824 } 0825 Cell avgRangeStart(e->sheet, e->ranges[2].col1, e->ranges[2].row1); 0826 return calc->averageIfs(avgRangeStart, c_Range, cond, lim); 0827 } 0828 0829 // 0830 // Function: avedev 0831 // 0832 Value func_avedev(valVector args, ValueCalc *calc, FuncExtra *) 0833 { 0834 Value result; 0835 calc->arrayWalk(args, result, awAveDev, calc->avg(args)); 0836 return calc->div(result, calc->count(args)); 0837 } 0838 0839 // 0840 // Function: betadist 0841 // 0842 Value func_betadist(valVector args, ValueCalc *calc, FuncExtra *) 0843 { 0844 Value x = args[0]; 0845 Value alpha = args[1]; 0846 Value beta = args[2]; 0847 0848 // default values parameter 4 - 6 0849 Value fA(0.0); 0850 Value fB(1.0); 0851 bool kum = true; 0852 0853 if (args.count() > 3) fA = args[3]; 0854 if (args.count() > 4) fB = args[4]; 0855 if (args.count() > 5) 0856 kum = calc->conv()->asInteger(args[5]).asInteger(); // 0 or 1 0857 0858 // constraints x < fA || x > fB || fA == fB || alpha <= 0.0 || beta <= 0.0 0859 if (calc->lower(x, fA) || calc->equal(fA, fB) || 0860 (!calc->greater(alpha, 0.0)) || !calc->greater(beta, 0.0)) 0861 return Value(0.0); 0862 0863 // constraints x > b 0864 if (calc->greater(x, fB)) { 0865 if (kum) 0866 return Value(1.0); 0867 else 0868 return Value(0.0); 0869 } 0870 0871 // scale = (x - fA) / (fB - fA) // prescaling 0872 Value scale = calc->div(calc->sub(x, fA), calc->sub(fB, fA)); 0873 0874 if (kum) 0875 return calc->GetBeta(scale, alpha, beta); 0876 else { 0877 Value res = calc->div(calc->mul(calc->GetGamma(alpha), calc->GetGamma(beta)), 0878 calc->GetGamma(calc->add(alpha, beta))); 0879 Value b1 = calc->pow(scale, calc->sub(alpha, Value(1.0))); 0880 Value b2 = calc->pow(calc->sub(Value(1.0), scale), calc->sub(beta, Value(1.0))); 0881 0882 return calc->mul(calc->mul(res, b1), b2); 0883 } 0884 } 0885 0886 // 0887 // Function: betainv 0888 // 0889 Value func_betainv(valVector args, ValueCalc *calc, FuncExtra *) 0890 { 0891 Value p = args[0]; 0892 Value alpha = args[1]; 0893 Value beta = args[2]; 0894 0895 // default values parameter 4 - 6 0896 Value fA(0.0); 0897 Value fB(1.0); 0898 0899 if (args.count() > 3) fA = args[3]; 0900 if (args.count() > 4) fB = args[4]; 0901 0902 Value result; 0903 0904 // constraints 0905 if (calc->lower(alpha, 0.0) || calc->lower(beta, 0.0) || 0906 calc->greater(p, 1.0) || calc->lower(p, 0.0) || calc->equal(fA, fB)) 0907 return Value::errorVALUE(); 0908 0909 bool convergenceError; 0910 0911 result = InverseIterator(func_betadist, valVector() << alpha << beta, calc).exec(p.asFloat(), 0.0, 1.0, convergenceError); 0912 0913 if (convergenceError) 0914 return Value::errorVALUE(); 0915 0916 result = calc->add(fA, calc->mul(result, calc->sub(fB, fA))); 0917 0918 return result; 0919 } 0920 0921 // 0922 // Function: bino 0923 // 0924 // kspread version 0925 // 0926 Value func_bino(valVector args, ValueCalc *calc, FuncExtra *) 0927 { 0928 Value n = args[0]; 0929 Value m = args[1]; 0930 Value comb = calc->combin(n, m); 0931 Value prob = args[2]; 0932 0933 if (calc->lower(prob, Value(0)) || calc->greater(prob, Value(1))) 0934 return Value::errorVALUE(); 0935 0936 // result = comb * pow (prob, m) * pow (1 - prob, n - m) 0937 Value pow1 = calc->pow(prob, m); 0938 Value pow2 = calc->pow(calc->sub(1, prob), calc->sub(n, m)); 0939 return calc->mul(comb, calc->mul(pow1, pow2)); 0940 } 0941 0942 // 0943 // Function: binomdist 0944 // 0945 // binomdist ( x, n, p, bool cumulative ) 0946 // 0947 Value func_binomdist(valVector args, ValueCalc *calc, FuncExtra *) 0948 { 0949 // TODO using approxfloor 0950 double x = calc->conv()->asFloat(args[0]).asFloat(); 0951 double n = calc->conv()->asFloat(args[1]).asFloat(); 0952 double p = calc->conv()->asFloat(args[2]).asFloat(); 0953 bool kum = calc->conv()->asInteger(args[3]).asInteger(); 0954 0955 debugSheets << "x= " << x << " n= " << n << " p= " << p; 0956 0957 // check constraints 0958 if (n < 0.0 || x < 0.0 || x > n || p < 0.0 || p > 1.0) 0959 return Value::errorVALUE(); 0960 0961 double res; 0962 double factor; 0963 double q; 0964 0965 if (kum) { 0966 // calculation of binom-distribution 0967 debugSheets << "calc distribution"; 0968 if (n == x) 0969 res = 1.0; 0970 else { 0971 q = 1.0 - p; 0972 factor = pow(q, n); 0973 if (factor == 0.0) { 0974 factor = pow(p, n); 0975 if (factor == 0.0) 0976 return Value::errorNA(); //SetNoValue(); 0977 else { 0978 res = 1.0 - factor; 0979 unsigned long max = (unsigned long)(n - x) - 1; 0980 for (unsigned long i = 0; i < max && factor > 0.0; ++i) { 0981 factor *= (n - i) / (i + 1) * q / p; 0982 res -= factor; 0983 } 0984 if (res < 0.0) 0985 res = 0.0; 0986 } 0987 } else { 0988 res = factor; 0989 unsigned long max = (unsigned long) x; 0990 for (unsigned long i = 0; i < max && factor > 0.0; ++i) { 0991 factor *= (n - i) / (i + 1) * p / q; 0992 res += factor; 0993 } 0994 } 0995 } 0996 } else { // density 0997 debugSheets << "calc density"; 0998 q = 1.0 - p; 0999 factor = pow(q, n); 1000 if (factor == 0.0) { 1001 factor = pow(p, n); 1002 if (factor == 0.0) 1003 return Value::errorNA(); //SetNoValue(); 1004 else { 1005 unsigned long max = (unsigned long)(n - x); 1006 for (unsigned long i = 0; i < max && factor > 0.0; ++i) 1007 factor *= (n - i) / (i + 1) * q / p; 1008 res = factor; 1009 } 1010 } else { 1011 unsigned long max = (unsigned long) x; 1012 for (unsigned long i = 0; i < max && factor > 0.0; ++i) 1013 factor *= (n - i) / (i + 1) * p / q; 1014 res = factor; 1015 } 1016 } 1017 return Value(res); 1018 } 1019 1020 // 1021 // Function: chidist 1022 // 1023 // returns the chi-distribution 1024 // 1025 Value func_chidist(valVector args, ValueCalc *calc, FuncExtra *) 1026 { 1027 Value fChi = args[0]; 1028 Value fDF = args[1]; 1029 1030 // fDF < 1 || fDF >= 1.0E5 1031 if (calc->lower(fDF, Value(1)) || (!calc->lower(fDF, Value(1.0E5)))) 1032 return Value::errorVALUE(); 1033 // fChi <= 0.0 1034 if (calc->lower(fChi, Value(0.0)) || (!calc->greater(fChi, Value(0.0)))) 1035 return Value(1.0); 1036 1037 // 1.0 - GetGammaDist (fChi / 2.0, fDF / 2.0, 1.0) 1038 return calc->sub(Value(1.0), calc->GetGammaDist(calc->div(fChi, 2.0), 1039 calc->div(fDF, 2.0), Value(1.0))); 1040 } 1041 1042 // 1043 // Function: combin 1044 // 1045 Value func_combin(valVector args, ValueCalc *calc, FuncExtra *) 1046 { 1047 if (calc->lower(args[1], Value(0.0)) || calc->lower(args[1], Value(0.0)) || calc->greater(args[1], args[0])) 1048 return Value::errorNUM(); 1049 1050 return calc->combin(args[0], args[1]); 1051 } 1052 1053 // 1054 // Function: combina 1055 // 1056 Value func_combina(valVector args, ValueCalc *calc, FuncExtra *) 1057 { 1058 if (calc->lower(args[1], Value(0.0)) || calc->lower(args[1], Value(0.0)) || calc->greater(args[1], args[0])) 1059 return Value::errorNUM(); 1060 1061 return calc->combin(calc->sub(calc->add(args[0], args[1]), Value(1.0)), args[1]); 1062 } 1063 1064 // 1065 // Function: confidence 1066 // 1067 // returns the confidence interval for a population mean 1068 // 1069 Value func_confidence(valVector args, ValueCalc *calc, FuncExtra *) 1070 { 1071 Value alpha = args[0]; 1072 Value sigma = args[1]; 1073 Value n = args[2]; 1074 1075 // sigma <= 0.0 || alpha <= 0.0 || alpha >= 1.0 || n < 1 1076 if ((!calc->greater(sigma, Value(0.0))) || (!calc->greater(alpha, Value(0.0))) || 1077 (!calc->lower(alpha, Value(1.0))) || calc->lower(n, Value(1))) 1078 return Value::errorVALUE(); 1079 1080 // g = gaussinv (1.0 - alpha / 2.0) 1081 Value g = calc->gaussinv(calc->sub(Value(1.0), calc->div(alpha, 2.0))); 1082 // g * sigma / sqrt (n) 1083 return calc->div(calc->mul(g, sigma), calc->sqrt(n)); 1084 } 1085 1086 // 1087 // function: correl_pop 1088 // 1089 Value func_correl_pop(valVector args, ValueCalc *calc, FuncExtra *) 1090 { 1091 Value covar = func_covar(args, calc, 0); 1092 Value stdevp1 = calc->stddevP(args[0]); 1093 Value stdevp2 = calc->stddevP(args[1]); 1094 1095 if (calc->isZero(stdevp1) || calc->isZero(stdevp2)) 1096 return Value::errorDIV0(); 1097 1098 // covar / (stdevp1 * stdevp2) 1099 return calc->div(covar, calc->mul(stdevp1, stdevp2)); 1100 } 1101 1102 // 1103 // function: covar 1104 // 1105 Value func_covar(valVector args, ValueCalc *calc, FuncExtra *) 1106 { 1107 Value avg1 = calc->avg(args[0]); 1108 Value avg2 = calc->avg(args[1]); 1109 int number = calc->count(args[0]); 1110 int number2 = calc->count(args[1]); 1111 1112 if (number2 <= 0 || number2 != number) 1113 return Value::errorVALUE(); 1114 1115 Value covar = func_covar_helper(args[0], args[1], calc, avg1, avg2); 1116 return calc->div(covar, number); 1117 } 1118 1119 // 1120 // function: devsq 1121 // 1122 Value func_devsq(valVector args, ValueCalc *calc, FuncExtra *) 1123 { 1124 Value res; 1125 calc->arrayWalk(args, res, calc->awFunc("devsq"), calc->avg(args, false)); 1126 return res; 1127 } 1128 1129 // 1130 // function: devsqa 1131 // 1132 Value func_devsqa(valVector args, ValueCalc *calc, FuncExtra *) 1133 { 1134 Value res; 1135 calc->arrayWalk(args, res, calc->awFunc("devsqa"), calc->avg(args)); 1136 return res; 1137 } 1138 1139 // 1140 // Function: expondist 1141 // 1142 // returns the exponential distribution 1143 // 1144 Value func_expondist(valVector args, ValueCalc *calc, FuncExtra *) 1145 { 1146 Value x = args[0]; 1147 Value lambda = args[1]; 1148 Value kum = args[2]; 1149 1150 Value result(0.0); 1151 1152 if (!calc->greater(lambda, 0.0)) 1153 return Value::errorVALUE(); 1154 1155 // ex = exp (-lambda * x) 1156 Value ex = calc->exp(calc->mul(calc->mul(lambda, -1), x)); 1157 if (calc->isZero(kum)) { //density 1158 if (!calc->lower(x, 0.0)) 1159 // lambda * ex 1160 result = calc->mul(lambda, ex); 1161 } else { //distribution 1162 if (calc->greater(x, 0.0)) 1163 // 1.0 - ex 1164 result = calc->sub(1.0, ex); 1165 } 1166 return result; 1167 } 1168 1169 // 1170 // Function: fdist 1171 // 1172 // returns the f-distribution 1173 // 1174 Value func_fdist(valVector args, ValueCalc *calc, FuncExtra *) 1175 { 1176 Value x = args[0]; 1177 Value fF1 = args[1]; 1178 Value fF2 = args[2]; 1179 1180 bool kum = true; 1181 if (args.count() > 3) 1182 kum = calc->conv()->asInteger(args[3]).asInteger(); 1183 1184 // check constraints 1185 // x < 0.0 || fF1 < 1 || fF2 < 1 || fF1 >= 1.0E10 || fF2 >= 1.0E10 1186 if (calc->lower(x, Value(0.0)) || calc->lower(fF1, Value(1)) || calc->lower(fF2, Value(1)) || 1187 (!calc->lower(fF1, Value(1.0E10))) || (!calc->lower(fF2, Value(1.0E10)))) 1188 return Value::errorVALUE(); 1189 1190 if (kum) { 1191 // arg = fF2 / (fF2 + fF1 * x) 1192 Value arg = calc->div(fF2, calc->add(fF2, calc->mul(fF1, x))); 1193 // alpha = fF2/2.0 1194 Value alpha = calc->div(fF2, 2.0); 1195 // beta = fF1/2.0 1196 Value beta = calc->div(fF1, 2.0); 1197 return calc->sub(Value(1), calc->GetBeta(arg, alpha, beta)); 1198 } else { 1199 // non-cumulative calculation 1200 if (calc->lower(x, Value(0.0))) 1201 return Value(0); 1202 else { 1203 double res = 0.0; 1204 double f1 = calc->conv()->asFloat(args[1]).asFloat(); 1205 double f2 = calc->conv()->asFloat(args[2]).asFloat(); 1206 double xx = calc->conv()->asFloat(args[0]).asFloat(); 1207 1208 1209 double tmp1 = (f1 / 2) * log(f1) + (f2 / 2) * log(f2); 1210 double tmp2 = calc->GetLogGamma(Value((f1 + f2) / 2)).asFloat(); 1211 double tmp3 = calc->GetLogGamma(Value(f1 / 2)).asFloat(); 1212 double tmp4 = calc->GetLogGamma(Value(f2 / 2)).asFloat(); 1213 1214 res = exp(tmp1 + tmp2 - tmp3 - tmp4) * pow(xx, (f1 / 2) - 1) * pow(f2 + f1 * xx, (-f1 / 2) - (f2 / 2)); 1215 return Value(res); 1216 } 1217 } 1218 } 1219 1220 // 1221 // Function: finv 1222 // 1223 // returns the inverse f-distribution 1224 // 1225 Value func_finv(valVector args, ValueCalc *calc, FuncExtra *) 1226 { 1227 Value p = args[0]; 1228 Value f1 = args[1]; 1229 Value f2 = args[2]; 1230 1231 Value result; 1232 1233 //TODO constraints 1234 // if ( calc->lower(DF, 1.0) || calc->greater(DF, 1.0E5) || 1235 // calc->greater(p, 1.0) || calc->lower(p,0.0) ) 1236 // return Value::errorVALUE(); 1237 1238 bool convergenceError; 1239 1240 result = InverseIterator(func_fdist, valVector() << f1 << f2 << Value(1), calc).exec(p.asFloat(), f1.asFloat() * 0.5, f1.asFloat(), convergenceError); 1241 1242 if (convergenceError) 1243 return Value::errorVALUE(); 1244 1245 return result; 1246 } 1247 1248 // 1249 // Function: fisher 1250 // 1251 // returns the Fisher transformation for x 1252 // 1253 Value func_fisher(valVector args, ValueCalc *calc, FuncExtra *) 1254 { 1255 // 0.5 * ln ((1.0 + fVal) / (1.0 - fVal)) 1256 Value fVal = args[0]; 1257 Value num = calc->div(calc->add(fVal, 1.0), calc->sub(1.0, fVal)); 1258 return calc->mul(calc->ln(num), 0.5); 1259 } 1260 1261 // 1262 // Function: fisherinv 1263 // 1264 // returns the inverse of the Fisher transformation for x 1265 // 1266 Value func_fisherinv(valVector args, ValueCalc *calc, FuncExtra *) 1267 { 1268 Value fVal = args[0]; 1269 // (exp (2.0 * fVal) - 1.0) / (exp (2.0 * fVal) + 1.0) 1270 Value ex = calc->exp(calc->mul(fVal, 2.0)); 1271 return calc->div(calc->sub(ex, 1.0), calc->add(ex, 1.0)); 1272 } 1273 1274 // 1275 // Function: FREQUENCY 1276 // 1277 Value func_frequency(valVector args, ValueCalc*, FuncExtra*) 1278 { 1279 const Value bins = args[1]; 1280 if (bins.columns() != 1) 1281 return Value::errorVALUE(); 1282 1283 // create a data vector 1284 QVector<double> data; 1285 for (uint v = 0; v < args[0].count(); ++v) { 1286 if (args[0].element(v).isNumber()) 1287 data.append(numToDouble(args[0].element(v).asFloat())); 1288 } 1289 1290 // no intervals given? 1291 if (bins.isEmpty()) 1292 return Value(data.count()); 1293 1294 // sort the data 1295 std::stable_sort(data.begin(), data.end()); 1296 1297 Value result(Value::Array); 1298 QVector<double>::ConstIterator begin = data.constBegin(); 1299 QVector<double>::ConstIterator it = data.constBegin(); 1300 for (uint v = 0; v < bins.count(); ++v) { 1301 if (!bins.element(v).isNumber()) 1302 continue; 1303 it = std::upper_bound(begin, data.constEnd(), bins.element(v).asFloat()); 1304 // exact match? 1305 if (*it == numToDouble(bins.element(v).asFloat())) 1306 ++it; 1307 // add the number of values in this interval to the result 1308 result.setElement(0, v, Value(static_cast<qint64>(it - begin))); 1309 begin = it; 1310 } 1311 // the remaining values 1312 result.setElement(0, bins.count(), Value(static_cast<qint64>(data.constEnd() - begin))); 1313 1314 return result; 1315 } 1316 1317 // 1318 // Function: FTEST 1319 // 1320 // TODO - check if parameters are arrays 1321 Value func_ftest(valVector args, ValueCalc *calc, FuncExtra*) 1322 { 1323 const Value matrixA = args[0]; 1324 const Value matrixB = args[1]; 1325 1326 double val = 0.0; // stores current array value 1327 double countA = 0.0; // 1328 double countB = 0.0; // 1329 double sumA = 0.0, sumSqrA = 0.0; 1330 double sumB = 0.0, sumSqrB = 0.0; 1331 1332 // matrixA 1333 for (uint v = 0; v < matrixA.count(); ++v) { 1334 if (matrixA.element(v).isNumber()) { 1335 ++countA; 1336 val = numToDouble(matrixA.element(v).asFloat()); 1337 sumA += val; // add sum 1338 sumSqrA += val * val; // add square 1339 } 1340 } 1341 1342 // matrixB 1343 for (uint v = 0; v < matrixB.count(); ++v) { 1344 if (matrixB.element(v).isNumber()) { 1345 ++countB; 1346 val = numToDouble(matrixB.element(v).asFloat()); 1347 sumB += val; // add sum 1348 sumSqrB += val * val; // add square 1349 } 1350 } 1351 1352 // check constraints 1353 if (countA < 2 || countB < 2) 1354 return Value::errorNA(); 1355 1356 double sA = (sumSqrA - sumA * sumA / countA) / (countA - 1.0); 1357 double sB = (sumSqrB - sumB * sumB / countB) / (countB - 1.0); 1358 1359 if (sA == 0.0 || sB == 0.0) 1360 return Value::errorNA(); 1361 1362 double x, r1, r2; 1363 1364 if (sA > sB) { 1365 x = sA / sB; 1366 r1 = countA - 1.0; 1367 r2 = countB - 1.0; 1368 } else { 1369 x = sB / sA; 1370 r1 = countB - 1.0; 1371 r2 = countA - 1.0; 1372 } 1373 1374 valVector param; 1375 param.append(Value(x)); 1376 param.append(Value(r1)); 1377 param.append(Value(r2)); 1378 1379 return calc->mul(Value(2.0), func_legacyfdist(param, calc, 0)); 1380 } 1381 1382 // 1383 // Function: gammadist 1384 // 1385 Value func_gammadist(valVector args, ValueCalc *calc, FuncExtra *) 1386 { 1387 Value x = args[0]; 1388 Value alpha = args[1]; 1389 Value beta = args[2]; 1390 int kum = calc->conv()->asInteger(args[3]).asInteger(); // 0 or 1 1391 1392 Value result; 1393 1394 if (calc->lower(x, 0.0) || (!calc->greater(alpha, 0.0)) || 1395 (!calc->greater(beta, 0.0))) 1396 return Value::errorVALUE(); 1397 1398 if (kum == 0) { //density 1399 Value G = calc->GetGamma(alpha); 1400 // result = pow (x, alpha - 1.0) / exp (x / beta) / pow (beta, alpha) / G 1401 Value pow1 = calc->pow(x, calc->sub(alpha, 1.0)); 1402 Value pow2 = calc->exp(calc->div(x, beta)); 1403 Value pow3 = calc->pow(beta, alpha); 1404 result = calc->div(calc->div(calc->div(pow1, pow2), pow3), G); 1405 } else 1406 result = calc->GetGammaDist(x, alpha, beta); 1407 1408 return Value(result); 1409 } 1410 1411 // 1412 // Function: gammainv 1413 // 1414 Value func_gammainv(valVector args, ValueCalc *calc, FuncExtra *) 1415 { 1416 Value p = args[0]; 1417 Value alpha = args[1]; 1418 Value beta = args[2]; 1419 1420 Value result; 1421 1422 // constraints 1423 if (calc->lower(alpha, 0.0) || calc->lower(beta, 0.0) || 1424 calc->lower(p, 0.0) || !calc->lower(p, 1.0)) 1425 return Value::errorVALUE(); 1426 1427 bool convergenceError; 1428 Value start = calc->mul(alpha, beta); 1429 1430 result = InverseIterator(func_gammadist, valVector() << alpha << beta << Value(1), calc).exec(p.asFloat(), start.asFloat() * 0.5, start.asFloat(), convergenceError); 1431 1432 if (convergenceError) 1433 return Value::errorVALUE(); 1434 1435 return result; 1436 } 1437 1438 // 1439 // Function: gammaln 1440 // 1441 // returns the natural logarithm of the gamma function 1442 // 1443 Value func_gammaln(valVector args, ValueCalc *calc, FuncExtra *) 1444 { 1445 if (calc->greater(args[0], Value(0.0))) 1446 return calc->GetLogGamma(args[0]); 1447 return Value::errorVALUE(); 1448 } 1449 1450 // 1451 // Function: gauss 1452 // 1453 Value func_gauss(valVector args, ValueCalc *calc, FuncExtra *) 1454 { 1455 //returns the integral values of the standard normal cumulative distribution 1456 return calc->gauss(args[0]); 1457 } 1458 1459 // 1460 // Function: growth 1461 // 1462 // GROWTH ( knownY [; [knownX] [; [newX] [; allowOsset = TRUE() ] ] ] ) 1463 // 1464 Value func_growth(valVector args, ValueCalc *calc, FuncExtra *) 1465 { 1466 debugSheets << "GROWTH"; // Debug 1467 Value known_Y = args[0]; 1468 1469 // default 1470 bool withOffset = true; 1471 1472 if (args.count() > 3) 1473 withOffset = calc->conv()->asInteger(args[3]).asInteger(); 1474 1475 // check constraints 1476 if (known_Y.isEmpty()) { 1477 debugSheets << "known_Y is empty"; 1478 return Value::errorNA(); 1479 } 1480 1481 // check if array known_Y contains only numbers 1482 for (uint i = 0; i < known_Y.count(); ++i) { 1483 if (!known_Y.element(i).isNumber()) { 1484 debugSheets << "count_Y (" << i << ") is non Value"; 1485 return Value::errorNA(); 1486 } 1487 } 1488 1489 uint cols_X, cols_Y; // columns in X and Y-Matrix 1490 uint rows_X, rows_Y; // rows in X and Y-Matrix 1491 1492 // stores count of elements in array 1493 // int count_Y = 0; 1494 // int count_X = 0; 1495 1496 // 1497 int nCase = 0; 1498 1499 // get size Y-Matrix 1500 rows_Y = known_Y.rows(); 1501 cols_Y = known_Y.columns(); 1502 debugSheets << "Y has " << rows_Y << " rows"; 1503 debugSheets << "Y has " << cols_Y << " cols"; 1504 1505 // convert all Value in known_Y into log 1506 for (uint r = 0; r < rows_Y; ++r) 1507 for (uint c = 0; c < cols_Y; ++c) { 1508 debugSheets << "col " << c << " row " << r << " log of Y(" << known_Y.element(c, r) << ") Value=" << calc->log(known_Y.element(c, r)); // Debug 1509 known_Y.setElement(c, r, calc->ln(known_Y.element(c, r))); 1510 } 1511 1512 Value known_X(Value::Array); 1513 1514 // 1515 // knownX is given ... 1516 // 1517 if (args.count() > 1) { 1518 // 1519 // get X-Matrix and print size 1520 // 1521 known_X = args[1]; 1522 1523 rows_X = known_Y.rows(); 1524 cols_X = known_X.columns(); 1525 debugSheets << "X has " << rows_X << " rows"; 1526 debugSheets << "X has " << cols_X << " cols"; 1527 1528 // 1529 // check if array known_X contains only numbers 1530 // 1531 for (uint i = 0; i < known_X.count(); ++i) { 1532 if (!known_X.element(i).isNumber()) { 1533 debugSheets << "count_X (" << i << ") is non Value"; 1534 return Value::errorNA(); 1535 } 1536 } 1537 1538 // 1539 // check for simple regression 1540 // 1541 if (cols_X == cols_Y && rows_X == rows_Y) 1542 nCase = 1; 1543 1544 else if (cols_Y != 1 && rows_Y != 1) { 1545 debugSheets << "Y-Matrix only has one row or column"; 1546 return Value::errorNA(); // TODO which errortype VALUE? 1547 } else if (cols_Y == 1) { 1548 // 1549 // row alignment 1550 // 1551 if (rows_X != rows_Y) { 1552 debugSheets << "--> row aligned"; 1553 debugSheets << "row sizes not equal"; 1554 return Value::errorNA(); 1555 } else { 1556 debugSheets << "--> row aligned"; 1557 nCase = 2; // row alignment 1558 } 1559 } 1560 1561 // 1562 // only column alignment left 1563 // 1564 else if (cols_X != cols_Y) { 1565 debugSheets << "--> col aligned"; 1566 debugSheets << "col sizes not equal"; 1567 return Value::errorNA(); 1568 } else { 1569 debugSheets << "--> col aligned"; 1570 nCase = 3; // column alignment 1571 } 1572 } else 1573 // 1574 // if known_X is empty it has to be set to 1575 // the sequence 1,2,3... n (n number of counts knownY) in one row. 1576 // 1577 { 1578 debugSheets << "fill X-Matrix with 0,1,2,3 .. sequence"; 1579 const int known_Y_count = known_Y.count(); 1580 for (int i = 0; i < known_Y_count; ++i) 1581 known_X.setElement(i, 0, Value(i)); 1582 1583 cols_X = cols_Y; 1584 rows_X = rows_Y; 1585 1586 // simple regression 1587 nCase = 1; 1588 } 1589 1590 Value newX(Value::Array); 1591 uint cols_newX, rows_newX; 1592 1593 if (args.count() < 3) { 1594 debugSheets << "no newX-Matrix --> copy X-Matrix"; 1595 cols_newX = cols_X; 1596 rows_newX = rows_X; 1597 newX = known_X; 1598 } else { 1599 newX = args[2]; 1600 // get dimensions 1601 cols_newX = newX.columns(); 1602 rows_newX = newX.rows(); 1603 1604 if ((nCase == 2 && cols_X != cols_newX) || (nCase == 3 && rows_X != rows_newX)) { 1605 debugSheets << "newX does not fit..."; 1606 return Value::errorNA(); 1607 } 1608 1609 // check if array newX contains only numbers 1610 for (uint i = 0; i < newX.count(); ++i) { 1611 if (!newX.element(i).isNumber()) { 1612 debugSheets << "newX (" << i << ") is non Value"; 1613 return Value::errorNA(); 1614 } 1615 } 1616 } 1617 1618 debugSheets << "known_X = " << known_X; 1619 debugSheets << "newX = " << newX; 1620 debugSheets << "newX has " << rows_newX << " rows"; 1621 debugSheets << "newX has " << cols_newX << " cols"; 1622 1623 // create the resulting matrix 1624 Value res(Value::Array); 1625 1626 // 1627 // simple regression 1628 // 1629 if (nCase == 1) { 1630 debugSheets << "Simple regression detected"; // Debug 1631 1632 double count = 0.0; 1633 double sumX = 0.0; 1634 double sumSqrX = 0.0; 1635 double sumY = 0.0; 1636 double sumSqrY = 0.0; 1637 double sumXY = 0.0; 1638 double valX, valY; 1639 1640 // 1641 // Gehe �ber Matrix Reihen/Spaltenweise 1642 // 1643 for (uint c = 0; c < cols_Y; ++c) { 1644 for (uint r = 0; r < rows_Y; ++r) { 1645 valX = known_X.element(c, r).asFloat(); 1646 valY = known_Y.element(c, r).asFloat(); 1647 sumX += valX; 1648 sumSqrX += valX * valX; 1649 sumY += valY; 1650 sumSqrY += valY * valY; 1651 sumXY += valX * valY; 1652 ++count; 1653 } 1654 } 1655 1656 if (count < 1.0) { 1657 debugSheets << "count less than 1.0"; 1658 return Value::errorNA(); 1659 } else { 1660 double f1 = count * sumXY - sumX * sumY; 1661 double X = count * sumSqrX - sumX * sumX; 1662 double b, m; 1663 1664 if (withOffset) { 1665 // with offset 1666 b = sumY / count - f1 / X * sumX / count; 1667 m = f1 / X; 1668 } else { 1669 // without offset 1670 b = 0.0; 1671 m = sumXY / sumSqrX; 1672 } 1673 1674 // 1675 // Fill result matrix 1676 // 1677 for (uint c = 0; c < cols_newX; ++c) { 1678 for (uint r = 0; r < rows_newX; ++r) { 1679 double result = 0.0; 1680 result = exp(newX.element(c, r).asFloat() * m + b); 1681 debugSheets << "res(" << c << "," << r << ") = " << result; 1682 res.setElement(c, r, Value(result)); 1683 } 1684 } 1685 } 1686 debugSheets << res; 1687 } else { 1688 if (nCase == 2) { 1689 debugSheets << "column alignment"; 1690 } else { 1691 debugSheets << "row alignment"; 1692 } 1693 } 1694 1695 return (res); // return array 1696 } 1697 1698 // 1699 // function: geomean 1700 // 1701 Value func_geomean(valVector args, ValueCalc *calc, FuncExtra *) 1702 { 1703 Value count = Value(calc->count(args)); 1704 Value prod = calc->product(args, Value(1.0)); 1705 if (calc->isZero(count)) 1706 return Value::errorDIV0(); 1707 return calc->pow(prod, calc->div(Value(1.0), count)); 1708 } 1709 1710 // 1711 // function: harmean 1712 // 1713 Value func_harmean(valVector args, ValueCalc *calc, FuncExtra *) 1714 { 1715 Value count(calc->count(args)); 1716 if (calc->isZero(count)) 1717 return Value::errorDIV0(); 1718 Value suminv; 1719 calc->arrayWalk(args, suminv, awSumInv, Value(0)); 1720 return calc->div(count, suminv); 1721 } 1722 1723 // 1724 // function: hypgeomdist 1725 // 1726 Value func_hypgeomdist(valVector args, ValueCalc *calc, FuncExtra *) 1727 { 1728 int x = calc->conv()->asInteger(args[0]).asInteger(); 1729 int n = calc->conv()->asInteger(args[1]).asInteger(); 1730 int M = calc->conv()->asInteger(args[2]).asInteger(); 1731 int N = calc->conv()->asInteger(args[3]).asInteger(); 1732 1733 double res = 0.0; 1734 1735 bool kum = false; 1736 if (args.count() > 4) 1737 kum = calc->conv()->asInteger(args[4]).asInteger(); 1738 1739 if (x < 0 || n < 0 || M < 0 || N < 0) 1740 return Value::errorVALUE(); 1741 1742 if (x > M || n > N) 1743 return Value::errorVALUE(); 1744 1745 if (kum) { 1746 for (int i = 0; i < x + 1; ++i) { 1747 Value d1 = calc->combin(M, i); 1748 Value d2 = calc->combin(N - M, n - i); 1749 Value d3 = calc->combin(N, n); 1750 1751 // d1 * d2 / d3 1752 res += calc->div(calc->mul(d1, d2), d3).asFloat(); 1753 } 1754 } else { 1755 Value d1 = calc->combin(M, x); 1756 Value d2 = calc->combin(N - M, n - x); 1757 Value d3 = calc->combin(N, n); 1758 1759 // d1 * d2 / d3 1760 res = calc->div(calc->mul(d1, d2), d3).asFloat(); 1761 } 1762 1763 return Value(res); 1764 } 1765 1766 // 1767 // Function: INTERCEPT 1768 // 1769 Value func_intercept(valVector args, ValueCalc* calc, FuncExtra*) 1770 { 1771 int numberY = calc->count(args[0]); 1772 int numberX = calc->count(args[1]); 1773 1774 if (numberY < 1 || numberX < 1 || numberY != numberX) 1775 return Value::errorVALUE(); 1776 1777 Value denominator; 1778 Value avgY = calc->avg(args[0]); 1779 Value avgX = calc->avg(args[1]); 1780 Value nominator = func_covar_helper(args[0], args[1], calc, avgY, avgX); 1781 calc->arrayWalk(args[1], denominator, calc->awFunc("devsq"), avgX); 1782 // result = Ey - SLOPE * Ex 1783 return calc->sub(avgY, calc->mul(calc->div(nominator, denominator), avgX)); 1784 } 1785 1786 // 1787 // function: kurtosis_est 1788 // 1789 Value func_kurtosis_est(valVector args, ValueCalc *calc, FuncExtra *) 1790 { 1791 int count = calc->count(args); 1792 if (count < 4) 1793 return Value::errorVALUE(); 1794 1795 Value avg = calc->avg(args); 1796 1797 Value stdev = calc->stddev(args, false); 1798 if (stdev.isZero()) 1799 return Value::errorDIV0(); 1800 1801 Value params(Value::Array); 1802 params.setElement(0, 0, avg); 1803 params.setElement(1, 0, stdev); 1804 Value x4; 1805 calc->arrayWalk(args, x4, awKurtosis, params); 1806 1807 // res = ( n*(n+1)*x4 - 3*(n-1)^3) / ( (n-3)*(n-2)*(n-1) ) 1808 int v1 = count * (count + 1); 1809 int v2 = 3 * (count - 1) * (count - 1) * (count - 1); 1810 int v3 = (count - 3) * (count - 2) * (count - 1); 1811 return calc->div(calc->sub(calc->mul(x4, v1), v2), v3); 1812 } 1813 1814 // 1815 // function: kurtosis_pop 1816 // 1817 Value func_kurtosis_pop(valVector args, ValueCalc *calc, FuncExtra *) 1818 { 1819 int count = calc->count(args); 1820 if (count < 4) 1821 return Value::errorVALUE(); 1822 1823 Value avg = calc->avg(args); 1824 Value stdev = calc->stddev(args, false); 1825 if (stdev.isZero()) 1826 return Value::errorDIV0(); 1827 1828 Value params(Value::Array); 1829 params.setElement(0, 0, avg); 1830 params.setElement(1, 0, stdev); 1831 Value x4; 1832 calc->arrayWalk(args, x4, awKurtosis, params); 1833 1834 // x4 / count - 3 1835 return calc->sub(calc->div(x4, count), 3); 1836 } 1837 1838 // 1839 // function: large 1840 // 1841 Value func_large(valVector args, ValueCalc *calc, FuncExtra *) 1842 { 1843 // does NOT support anything other than doubles !!! 1844 int k = calc->conv()->asInteger(args[1]).asInteger(); 1845 if (k < 1) 1846 return Value::errorVALUE(); 1847 1848 List array; 1849 int number = 1; 1850 1851 func_array_helper(args[0], calc, array, number); 1852 1853 if (k >= number || number - k - 1 >= array.count()) 1854 return Value::errorVALUE(); 1855 1856 std::sort(array.begin(), array.end()); 1857 double d = array.at(number - k - 1); 1858 return Value(d); 1859 } 1860 1861 // 1862 // Function: legacaychidist 1863 // 1864 // returns the chi-distribution 1865 // 1866 Value func_legacychidist(valVector args, ValueCalc *calc, FuncExtra *) 1867 { 1868 Value fChi = args[0]; 1869 Value fDF = args[1]; 1870 1871 // fDF < 1 || fDF >= 1.0E5 1872 if (calc->lower(fDF, Value(1)) || (!calc->lower(fDF, Value(1.0E5)))) 1873 return Value::errorVALUE(); 1874 // fChi <= 0.0 1875 if (calc->lower(fChi, Value(0.0)) || (!calc->greater(fChi, Value(0.0)))) 1876 return Value(1.0); 1877 1878 // 1.0 - GetGammaDist (fChi / 2.0, fDF / 2.0, 1.0) 1879 return calc->sub(Value(1.0), calc->GetGammaDist(calc->div(fChi, 2.0), 1880 calc->div(fDF, 2.0), Value(1.0))); 1881 } 1882 1883 // 1884 // Function: legacaychiinv 1885 // 1886 // returns the inverse chi-distribution 1887 // 1888 Value func_legacychiinv(valVector args, ValueCalc *calc, FuncExtra *) 1889 { 1890 Value p = args[0]; 1891 Value DF = args[1]; 1892 1893 Value result; 1894 1895 // constraints 1896 if (calc->lower(DF, 1.0) || calc->greater(DF, 1.0E5) || 1897 calc->greater(p, 1.0) || calc->lower(p, 0.0)) 1898 return Value::errorVALUE(); 1899 1900 bool convergenceError; 1901 1902 result = InverseIterator(func_legacychidist, valVector() << DF, calc).exec(p.asFloat(), DF.asFloat() * 0.5, DF.asFloat(), convergenceError); 1903 1904 if (convergenceError) 1905 return Value::errorVALUE(); 1906 1907 return result; 1908 } 1909 1910 // 1911 // Function: legacy.fdist 1912 // 1913 // returns the f-distribution 1914 // 1915 Value func_legacyfdist(valVector args, ValueCalc *calc, FuncExtra *) 1916 { 1917 Value x = args[0]; 1918 Value fF1 = args[1]; 1919 Value fF2 = args[2]; 1920 1921 // x < 0.0 || fF1 < 1 || fF2 < 1 || fF1 >= 1.0E10 || fF2 >= 1.0E10 1922 if (calc->lower(x, Value(0.0)) || calc->lower(fF1, Value(1)) || calc->lower(fF2, Value(1)) || 1923 (!calc->lower(fF1, Value(1.0E10))) || (!calc->lower(fF2, Value(1.0E10)))) 1924 return Value::errorVALUE(); 1925 1926 // arg = fF2 / (fF2 + fF1 * x) 1927 Value arg = calc->div(fF2, calc->add(fF2, calc->mul(fF1, x))); 1928 // alpha = fF2/2.0 1929 Value alpha = calc->div(fF2, 2.0); 1930 // beta = fF1/2.0 1931 Value beta = calc->div(fF1, 2.0); 1932 return calc->GetBeta(arg, alpha, beta); 1933 } 1934 1935 // 1936 // Function: legacyfinv 1937 // 1938 // returns the inverse legacy f-distribution 1939 // 1940 Value func_legacyfinv(valVector args, ValueCalc *calc, FuncExtra *) 1941 { 1942 Value p = args[0]; 1943 Value f1 = args[1]; 1944 Value f2 = args[2]; 1945 1946 Value result; 1947 1948 //TODO constraints 1949 // if ( calc->lower(DF, 1.0) || calc->greater(DF, 1.0E5) || 1950 // calc->greater(p, 1.0) || calc->lower(p,0.0) ) 1951 // return Value::errorVALUE(); 1952 1953 bool convergenceError; 1954 1955 result = InverseIterator(func_legacyfdist, valVector() << f1 << f2, calc).exec(p.asFloat(), f1.asFloat() * 0.5, f1.asFloat(), convergenceError); 1956 1957 if (convergenceError) 1958 return Value::errorVALUE(); 1959 1960 return result; 1961 } 1962 1963 // 1964 // function: loginv 1965 // 1966 Value func_loginv(valVector args, ValueCalc *calc, FuncExtra *) 1967 { 1968 Value p = args[0]; 1969 1970 // defaults 1971 Value m = Value(0.0); 1972 Value s = Value(1.0); 1973 1974 if (args.count() > 1) 1975 m = args[1]; 1976 if (args.count() > 2) 1977 s = args[2]; 1978 1979 if (calc->lower(p, Value(0)) || calc->greater(p, Value(1))) 1980 return Value::errorVALUE(); 1981 1982 if (!calc->greater(s, Value(0))) 1983 return Value::errorVALUE(); 1984 1985 Value result(0.0); 1986 if (calc->equal(p, Value(1))) //p==1 1987 result = Value::errorVALUE(); 1988 else if (calc->greater(p, Value(0))) { //p>0 1989 Value gaussInv = calc->gaussinv(p); 1990 // exp (gaussInv * s + m) 1991 result = calc->exp(calc->add(calc->mul(s, gaussInv), m)); 1992 } 1993 1994 return result; 1995 } 1996 1997 // 1998 // Function: lognormdist 1999 // 2000 // returns the cumulative lognormal distribution 2001 // 2002 Value func_lognormdist(valVector args, ValueCalc *calc, FuncExtra *) 2003 { 2004 // defaults 2005 Value mue = Value(0); 2006 Value sigma = Value(1); 2007 bool kum = true; 2008 2009 Value x = args[0]; 2010 if (args.count() > 1) 2011 mue = args[1]; 2012 if (args.count() > 2) 2013 sigma = args[2]; 2014 if (args.count() > 3) 2015 kum = calc->conv()->asInteger(args[3]).asInteger(); 2016 2017 if (!kum) { 2018 // TODO implement me !!! 2019 return Value::errorVALUE(); 2020 2021 // check constraints 2022 if (!calc->greater(sigma, 0.0) || (!calc->greater(x, 0.0))) 2023 return Value::errorVALUE(); 2024 } 2025 // non-cumulative 2026 // check constraints 2027 if (calc->lower(x, Value(0.0))) 2028 return Value(0.0); 2029 2030 // (ln(x) - mue) / sigma 2031 Value Y = calc->div(calc->sub(calc->ln(x), mue), sigma); 2032 return calc->add(calc->gauss(Y), 0.5); 2033 } 2034 2035 // 2036 // Function: MEDIAN 2037 // 2038 Value func_median(valVector args, ValueCalc *calc, FuncExtra *) 2039 { 2040 // does NOT support anything other than doubles !!! 2041 List array; 2042 int number = 0; 2043 2044 for (int i = 0; i < args.count(); ++i) 2045 func_array_helper(args[i], calc, array, number); 2046 2047 if (number == 0) 2048 return Value::errorVALUE(); 2049 2050 std::sort(array.begin(), array.end()); 2051 double d; 2052 if (number % 2) // odd 2053 d = array.at((number - 1) / 2); 2054 else // even 2055 d = 0.5 * (array.at(number / 2 - 1) + array.at(number / 2)); 2056 return Value(d); 2057 } 2058 2059 // 2060 // function: mode 2061 // 2062 Value func_mode(valVector args, ValueCalc *calc, FuncExtra *) 2063 { 2064 // does NOT support anything other than doubles !!! 2065 ContentSheet sh; 2066 for (int i = 0; i < args.count(); ++i) 2067 func_mode_helper(args[i], calc, sh); 2068 2069 // retrieve value with max.count 2070 int maxcount = 0; 2071 double max = 0.0; 2072 2073 // check if there is a difference in frequency 2074 ContentSheet::ConstIterator it = sh.constBegin(); 2075 double last = it.value(); // init last with 1st value 2076 bool nodiff = true; // reset flag 2077 2078 for ( ; it != sh.constEnd(); ++it) { 2079 if (it.value() > maxcount) { 2080 max = it.key(); 2081 maxcount = it.value(); 2082 } 2083 if (last != it.value()) 2084 nodiff = false; // set flag 2085 } 2086 2087 // if no diff return error 2088 if (nodiff) 2089 return Value::errorNUM(); 2090 else 2091 return Value(max); 2092 } 2093 2094 // 2095 // function: negbinomdist 2096 // 2097 Value func_negbinomdist(valVector args, ValueCalc *calc, FuncExtra *) 2098 { 2099 double x = calc->conv()->asFloat(args[0]).asFloat(); 2100 double r = calc->conv()->asFloat(args[1]).asFloat(); 2101 double p = calc->conv()->asFloat(args[2]).asFloat(); 2102 2103 if (r < 0.0 || x < 0.0 || p < 0.0 || p > 1.0) 2104 return Value::errorVALUE(); 2105 2106 double q = 1.0 - p; 2107 double res = pow(p, r); 2108 2109 for (double i = 0.0; i < x; ++i) 2110 res *= (i + r) / (i + 1.0) * q; 2111 2112 return Value(res); 2113 // int x = calc->conv()->asInteger (args[0]).asInteger(); 2114 // int r = calc->conv()->asInteger (args[1]).asInteger(); 2115 // Value p = args[2]; 2116 // 2117 // if ((x + r - 1) <= 0) 2118 // return Value::errorVALUE(); 2119 // if (calc->lower (p, Value(0)) || calc->greater (p, Value(1))) 2120 // return Value::errorVALUE(); 2121 // 2122 // Value d1 = calc->combin (x + r - 1, r - 1); 2123 // // d2 = pow (p, r) * pow (1 - p, x) 2124 // Value d2 = calc->mul (calc->pow (p, r), 2125 // calc->pow (calc->sub (Value(1), p), x)); 2126 // 2127 // return calc->mul (d1, d2); 2128 } 2129 2130 // 2131 // Function: normdist 2132 // 2133 // returns the normal cumulative distribution 2134 // 2135 Value func_normdist(valVector args, ValueCalc *calc, FuncExtra *) 2136 { 2137 Value x = args[0]; 2138 Value mue = args[1]; 2139 Value sigma = args[2]; 2140 Value k = args[3]; 2141 2142 if (!calc->greater(sigma, 0.0)) 2143 return Value::errorVALUE(); 2144 2145 // (x - mue) / sigma 2146 Value Y = calc->div(calc->sub(x, mue), sigma); 2147 if (calc->isZero(k)) // density 2148 return calc->div(calc->phi(Y), sigma); 2149 else // distribution 2150 return calc->add(calc->gauss(Y), 0.5); 2151 } 2152 2153 // 2154 // Function: norminv 2155 // 2156 // returns the inverse of the normal cumulative distribution 2157 // 2158 Value func_norminv(valVector args, ValueCalc *calc, FuncExtra *) 2159 { 2160 Value x = args[0]; 2161 Value mue = args[1]; 2162 Value sigma = args[2]; 2163 2164 if (!calc->greater(sigma, 0.0)) 2165 return Value::errorVALUE(); 2166 if (!(calc->greater(x, 0.0) && calc->lower(x, 1.0))) 2167 return Value::errorVALUE(); 2168 2169 // gaussinv (x)*sigma + mue 2170 return calc->add(calc->mul(calc->gaussinv(x), sigma), mue); 2171 } 2172 2173 // 2174 // Function: normsinv 2175 // 2176 // returns the inverse of the standard normal cumulative distribution 2177 // 2178 Value func_normsinv(valVector args, ValueCalc *calc, FuncExtra *) 2179 { 2180 Value x = args[0]; 2181 if (!(calc->greater(x, 0.0) && calc->lower(x, 1.0))) 2182 return Value::errorVALUE(); 2183 2184 return calc->gaussinv(x); 2185 } 2186 2187 // 2188 // Function: percentile 2189 // 2190 // PERCENTILE( data set; alpha ) 2191 // 2192 Value func_percentile(valVector args, ValueCalc *calc, FuncExtra*) 2193 { 2194 double alpha = numToDouble(calc->conv()->toFloat(args[1])); 2195 2196 // create array - does NOT support anything other than doubles !!! 2197 List array; 2198 int number = 0; 2199 2200 func_array_helper(args[0], calc, array, number); 2201 2202 // check constraints - number of values must be > 0 and flag >0 <=4 2203 if (number == 0) 2204 return Value::errorNA(); // or VALUE? 2205 if (alpha < -1e-9 || alpha > 1 + 1e-9) 2206 return Value::errorVALUE(); 2207 2208 // sort values 2209 std::sort(array.begin(), array.end()); 2210 2211 if (number == 1) 2212 return Value(array[0]); // only one value 2213 else { 2214 double r = alpha * (number - 1); 2215 int index = ::floor(r); 2216 double d = r - index; 2217 return Value(array[index] + d * (array[index+1] - array[index])); 2218 } 2219 } 2220 2221 // 2222 // Function: permutationa 2223 // 2224 // returns the number of ordered permutations (with repetition) 2225 Value func_permutationa(valVector args, ValueCalc *calc, FuncExtra *) 2226 { 2227 int n = calc->conv()->toInteger(args[0]); 2228 int k = calc->conv()->toInteger(args[1]); 2229 if (n < 0 || k < 0) 2230 return Value::errorVALUE(); 2231 return calc->pow(Value(n), k); 2232 } 2233 2234 // 2235 // Function: phi 2236 // 2237 // distribution function for a standard normal distribution 2238 // 2239 Value func_phi(valVector args, ValueCalc *calc, FuncExtra *) 2240 { 2241 return calc->phi(args[0]); 2242 } 2243 2244 // 2245 // Function: poisson 2246 // 2247 // returns the Poisson distribution 2248 // 2249 Value func_poisson(valVector args, ValueCalc *calc, FuncExtra *) 2250 { 2251 Value x = args[0]; 2252 Value lambda = args[1]; 2253 Value kum = args[2]; 2254 2255 // lambda < 0.0 || x < 0.0 2256 if (calc->lower(lambda, Value(0.0)) || calc->lower(x, Value(0.0))) 2257 return Value::errorVALUE(); 2258 2259 Value result; 2260 2261 // ex = exp (-lambda) 2262 Value ex = calc->exp(calc->mul(lambda, -1)); 2263 2264 if (calc->isZero(kum)) { // density 2265 if (calc->isZero(lambda)) 2266 result = Value(0); 2267 else 2268 // ex * pow (lambda, x) / fact (x) 2269 result = calc->div(calc->mul(ex, calc->pow(lambda, x)), calc->fact(x)); 2270 } else { // distribution 2271 if (calc->isZero(lambda)) 2272 result = Value(1); 2273 else { 2274 result = Value(1.0); 2275 Value fFak(1.0); 2276 qint64 nEnd = calc->conv()->asInteger(x).asInteger(); 2277 for (qint64 i = 1; i <= nEnd; ++i) { 2278 // fFak *= i 2279 fFak = calc->mul(fFak, (int)i); 2280 // result += pow (lambda, i) / fFak 2281 result = calc->add(result, calc->div(calc->pow(lambda, (int)i), fFak)); 2282 } 2283 result = calc->mul(result, ex); 2284 } 2285 } 2286 2287 return result; 2288 } 2289 2290 // 2291 // Function: rank 2292 // 2293 // rank(rank; ref.;sort order) 2294 Value func_rank(valVector args, ValueCalc *calc, FuncExtra*) 2295 { 2296 double x = calc->conv()->asFloat(args[0]).asFloat(); 2297 2298 // default 2299 bool descending = true; 2300 2301 double count = 1.0; 2302 double val = 0.0; 2303 bool valid = false; // flag 2304 2305 // opt. parameter 2306 if (args.count() > 2) 2307 descending = !calc->conv()->asInteger(args[2]).asInteger(); 2308 2309 // does NOT support anything other than doubles !!! 2310 List array; 2311 int number = 0; 2312 2313 func_array_helper(args[1], calc, array, number); 2314 2315 // sort array 2316 std::sort(array.begin(), array.end()); 2317 2318 for (int i = 0; i < array.count(); ++i) { 2319 if (descending) 2320 val = array[array.count()-count]; 2321 else 2322 val = array[i]; 2323 2324 //debugSheets<<"count ="<<count<<" val = "<<val<<" x = "<<x; 2325 2326 if (x == val) { 2327 valid = true; 2328 break; 2329 } else if ((!descending && x > val) || (descending && x < val)) 2330 ++count; 2331 } 2332 2333 if (valid) 2334 return Value(count); 2335 else 2336 return Value::errorNA(); 2337 } 2338 2339 // 2340 // Function: rsq 2341 // 2342 Value func_rsq(valVector args, ValueCalc *calc, FuncExtra*) 2343 { 2344 const Value matrixA = args[0]; 2345 const Value matrixB = args[1]; 2346 2347 // check constraints - arrays must have the same size 2348 if (matrixA.count() != matrixB.count()) 2349 return Value::errorVALUE(); 2350 2351 double valA = 0.0; // stores current array value 2352 double valB = 0.0; // stores current array value 2353 double count = 0.0; // count fields with numbers 2354 double sumA = 0.0, sumSqrA = 0.0; 2355 double sumB = 0.0, sumSqrB = 0.0; 2356 double sumAB = 0.0; 2357 2358 // calc sums 2359 for (uint v = 0; v < matrixA.count(); ++v) { 2360 Value vA(calc->conv()->asFloat(matrixA.element(v))); 2361 Value vB(calc->conv()->asFloat(matrixB.element(v))); 2362 2363 // TODO add unittest for check 2364 if (!vA.isError() && !vB.isError()) {// only if numbers are in both fields 2365 valA = calc->conv()->asFloat(matrixA.element(v)).asFloat(); 2366 valB = calc->conv()->asFloat(matrixB.element(v)).asFloat(); 2367 ++count; 2368 //debugSheets<<"valA ="<<valA<<" valB ="<<valB; 2369 2370 // value A 2371 sumA += valA; // add sum 2372 sumSqrA += valA * valA; // add square 2373 // value B 2374 sumB += valB; // add sum 2375 sumSqrB += valB * valB; // add square 2376 2377 sumAB += valA * valB; 2378 } 2379 } 2380 2381 // check constraints 2382 if (count < 2) 2383 return Value::errorNA(); 2384 else 2385 return Value((count*sumAB - sumA*sumB) *(count*sumAB - sumA*sumB) / 2386 (count*sumSqrA - sumA*sumA) / (count*sumSqrB - sumB*sumB)); 2387 } 2388 2389 // 2390 // Function: quartile 2391 // 2392 // QUARTILE( data set; flag ) 2393 // 2394 // flag: 2395 // 0 equals MIN() 2396 // 1 25th percentile 2397 // 2 50th percentile equals MEDIAN() 2398 // 3 75th percentile 2399 // 4 equals MAX() 2400 // 2401 Value func_quartile(valVector args, ValueCalc *calc, FuncExtra*) 2402 { 2403 int flag = calc->conv()->asInteger(args[1]).asInteger(); 2404 2405 // create array - does NOT support anything other than doubles !!! 2406 List array; 2407 int number = 0; 2408 2409 func_array_helper(args[0], calc, array, number); 2410 2411 // check constraints - number of values must be > 0 and flag >0 <=4 2412 if (number == 0) 2413 return Value::errorNA(); // or VALUE? 2414 if (flag < 0 || flag > 4) 2415 return Value::errorVALUE(); 2416 2417 // sort values 2418 std::sort(array.begin(), array.end()); 2419 2420 if (number == 1) 2421 return Value(array[0]); // only one value 2422 else { 2423 // 2424 // flag 0 -> MIN() 2425 // 2426 if (flag == 0) 2427 return Value(array[0]); 2428 2429 // 2430 // flag 1 -> 25th percentile 2431 // 2432 else if (flag == 1) { 2433 int nIndex = ::floor(0.25 * (number - 1)); 2434 double diff = 0.25 * (number - 1) - ::floor(0.25 * (number - 1)); 2435 2436 if (diff == 0.0) 2437 return Value(array[nIndex]); 2438 else 2439 return Value(array[nIndex] + diff*(array[nIndex+1] - array[nIndex])); 2440 } 2441 2442 // 2443 // flag 2 -> 50th percentile equals MEDIAN() 2444 // 2445 else if (flag == 2) { 2446 if (number % 2 == 0) 2447 return Value((array[number/2-1] + array[number/2]) / 2.0); 2448 else 2449 return Value(array[(number-1)/2]); 2450 } 2451 2452 // 2453 // flag 3 -> 75thpercentile 2454 // 2455 else if (flag == 3) { 2456 int nIndex = ::floor(0.75 * (number - 1)); 2457 double diff = 0.75 * (number - 1) - ::floor(0.75 * (number - 1)); 2458 2459 if (diff == 0.0) 2460 return Value(array[nIndex]); 2461 else 2462 return Value(array[nIndex] + diff*(array[nIndex+1] - array[nIndex])); 2463 } 2464 2465 // 2466 // flag 4 -> equals MAX() 2467 // 2468 else 2469 return Value(array[number-1]); 2470 } 2471 } 2472 2473 2474 // 2475 // function: skew_est 2476 // 2477 Value func_skew_est(valVector args, ValueCalc *calc, FuncExtra *) 2478 { 2479 int number = calc->count(args); 2480 Value avg = calc->avg(args); 2481 if (number < 3) 2482 return Value::errorVALUE(); 2483 2484 Value res = calc->stddev(args, avg); 2485 if (res.isZero()) 2486 return Value::errorVALUE(); 2487 2488 Value params(Value::Array); 2489 params.setElement(0, 0, avg); 2490 params.setElement(1, 0, res); 2491 Value tskew; 2492 calc->arrayWalk(args, tskew, awSkew, params); 2493 2494 // ((tskew * number) / (number-1)) / (number-2) 2495 return calc->div(calc->div(calc->mul(tskew, number), number - 1), number - 2); 2496 } 2497 2498 // 2499 // function: skew_pop 2500 // 2501 Value func_skew_pop(valVector args, ValueCalc *calc, FuncExtra *) 2502 { 2503 int number = calc->count(args); 2504 Value avg = calc->avg(args); 2505 if (number < 1) 2506 return Value::errorVALUE(); 2507 2508 Value res = calc->stddevP(args, avg); 2509 if (res.isZero()) 2510 return Value::errorVALUE(); 2511 2512 Value params(Value::Array); 2513 params.setElement(0, 0, avg); 2514 params.setElement(1, 0, res); 2515 Value tskew; 2516 calc->arrayWalk(args, tskew, awSkew, params); 2517 2518 // tskew / number 2519 return calc->div(tskew, (double)number); 2520 } 2521 2522 // 2523 // Function: SLOPE 2524 // 2525 Value func_slope(valVector args, ValueCalc* calc, FuncExtra*) 2526 { 2527 int numberY = calc->count(args[0]); 2528 int numberX = calc->count(args[1]); 2529 2530 if (numberY < 1 || numberX < 1 || numberY != numberX) 2531 return Value::errorVALUE(); 2532 2533 Value denominator; 2534 Value avgY = calc->avg(args[0]); 2535 Value avgX = calc->avg(args[1]); 2536 Value nominator = func_covar_helper(args[0], args[1], calc, avgY, avgX); 2537 calc->arrayWalk(args[1], denominator, calc->awFunc("devsq"), avgX); 2538 return calc->div(nominator, denominator); 2539 } 2540 2541 // 2542 // function: small 2543 // 2544 Value func_small(valVector args, ValueCalc *calc, FuncExtra *) 2545 { 2546 // does NOT support anything other than doubles !!! 2547 int k = calc->conv()->asInteger(args[1]).asInteger(); 2548 if (k < 1) 2549 return Value::errorVALUE(); 2550 2551 List array; 2552 int number = 1; 2553 2554 func_array_helper(args[0], calc, array, number); 2555 2556 if (k > number || k - 1 >= array.count()) 2557 return Value::errorVALUE(); 2558 2559 std::sort(array.begin(), array.end()); 2560 double d = array.at(k - 1); 2561 return Value(d); 2562 } 2563 2564 // 2565 // function: standardize 2566 // 2567 Value func_standardize(valVector args, ValueCalc *calc, FuncExtra *) 2568 { 2569 Value x = args[0]; 2570 Value m = args[1]; 2571 Value s = args[2]; 2572 2573 if (!calc->greater(s, Value(0))) // s must be >0 2574 return Value::errorVALUE(); 2575 2576 // (x - m) / s 2577 return calc->div(calc->sub(x, m), s); 2578 } 2579 2580 // 2581 // Function: stddev 2582 // 2583 Value func_stddev(valVector args, ValueCalc *calc, FuncExtra *) 2584 { 2585 return calc->stddev(args, false); 2586 } 2587 2588 // 2589 // Function: stddeva 2590 // 2591 Value func_stddeva(valVector args, ValueCalc *calc, FuncExtra *) 2592 { 2593 return calc->stddev(args); 2594 } 2595 2596 // 2597 // Function: stddevp 2598 // 2599 Value func_stddevp(valVector args, ValueCalc *calc, FuncExtra *) 2600 { 2601 return calc->stddevP(args, false); 2602 } 2603 2604 // 2605 // Function: stddevpa 2606 // 2607 Value func_stddevpa(valVector args, ValueCalc *calc, FuncExtra *) 2608 { 2609 return calc->stddevP(args); 2610 } 2611 2612 // 2613 // Function: normsdist 2614 // 2615 Value func_stdnormdist(valVector args, ValueCalc *calc, FuncExtra *) 2616 { 2617 //returns the cumulative lognormal distribution, mue=0, sigma=1 2618 return calc->add(calc->gauss(args[0]), 0.5); 2619 } 2620 2621 // 2622 // Function: STEYX 2623 // 2624 Value func_steyx(valVector args, ValueCalc* calc, FuncExtra*) 2625 { 2626 int number = calc->count(args[0]); 2627 2628 if (number < 1 || number != calc->count(args[1])) 2629 return Value::errorVALUE(); 2630 2631 Value varX, varY; 2632 Value avgY = calc->avg(args[0]); 2633 Value avgX = calc->avg(args[1]); 2634 Value cov = func_covar_helper(args[0], args[1], calc, avgY, avgX); 2635 calc->arrayWalk(args[0], varY, calc->awFunc("devsq"), avgY); 2636 calc->arrayWalk(args[1], varX, calc->awFunc("devsq"), avgX); 2637 Value numerator = calc->sub(varY, calc->div(calc->sqr(cov), varX)); 2638 Value denominator = calc->sub(Value(number), 2); 2639 return calc->sqrt(calc->div(numerator, denominator)); 2640 } 2641 2642 // 2643 // Function: sumproduct 2644 // 2645 Value func_sumproduct(valVector args, ValueCalc *calc, FuncExtra *) 2646 { 2647 Value result; 2648 calc->twoArrayWalk(args[0], args[1], result, tawSumproduct); 2649 return result; 2650 } 2651 2652 // 2653 // Function: sumx2py2 2654 // 2655 Value func_sumx2py2(valVector args, ValueCalc *calc, FuncExtra *) 2656 { 2657 Value result; 2658 calc->twoArrayWalk(args[0], args[1], result, tawSumx2py2); 2659 return result; 2660 } 2661 2662 // 2663 // Function: sumx2my2 2664 // 2665 Value func_sumx2my2(valVector args, ValueCalc *calc, FuncExtra *) 2666 { 2667 Value result; 2668 calc->twoArrayWalk(args[0], args[1], result, tawSumx2my2); 2669 return result; 2670 } 2671 2672 // 2673 // Function: SUMXMY2 2674 // 2675 Value func_sumxmy2(valVector args, ValueCalc *calc, FuncExtra *) 2676 { 2677 Value result; 2678 calc->twoArrayWalk(args[0], args[1], result, tawSumxmy2); 2679 return result; 2680 } 2681 2682 // 2683 // Function: tdist 2684 // 2685 // returns the t-distribution 2686 // 2687 Value func_tdist(valVector args, ValueCalc *calc, FuncExtra *) 2688 { 2689 Value T = args[0]; 2690 Value fDF = args[1]; 2691 int flag = calc->conv()->asInteger(args[2]).asInteger(); 2692 2693 if (calc->lower(fDF, Value(1)) || (flag != 1 && flag != 2)) 2694 return Value::errorVALUE(); 2695 2696 // arg = fDF / (fDF + T * T) 2697 Value arg = calc->div(fDF, calc->add(fDF, calc->sqr(T))); 2698 2699 Value R; 2700 R = calc->mul(calc->GetBeta(arg, calc->div(fDF, 2.0), Value(0.5)), 0.5); 2701 2702 if (flag == 1) 2703 return R; 2704 return calc->mul(R, 2); // flag is 2 here 2705 } 2706 2707 // 2708 // Function: tinv 2709 // 2710 // returns the inverse t-distribution 2711 // 2712 Value func_tinv(valVector args, ValueCalc *calc, FuncExtra *) 2713 { 2714 Value p = args[0]; 2715 Value DF = args[1]; 2716 2717 Value result; 2718 2719 // constraints 2720 if (calc->lower(DF, 1.0) || calc->greater(DF, 1.0E5) || 2721 calc->greater(p, 1.0) || calc->lower(p, 0.0)) 2722 return Value::errorVALUE(); 2723 2724 bool convergenceError; 2725 2726 result = InverseIterator(func_tdist, valVector() << DF << Value(2), calc).exec(p.asFloat(), DF.asFloat() * 0.5, DF.asFloat(), convergenceError); 2727 2728 if (convergenceError) 2729 return Value::errorVALUE(); 2730 2731 return result; 2732 } 2733 2734 // 2735 // Function: trend 2736 // 2737 // TREND ( knownY [; [knownX] [; [newX] [; allowOsset = TRUE() ] ] ] ) 2738 // 2739 // TODO - check do we need 2d arrays? 2740 // 2741 Value func_trend(valVector args, ValueCalc *calc, FuncExtra *) 2742 { 2743 // default 2744 bool withOffset = true; 2745 2746 if (args.count() > 3) 2747 withOffset = calc->conv()->asInteger(args[3]).asInteger(); 2748 2749 List knownY, knownX, newX; 2750 int knownXcount = 0, newXcount = 0; 2751 2752 // 2753 // knownX 2754 // 2755 if (args[1].isEmpty()) { 2756 // if knownX is empty it has to be set to the sequence 1,2,3... n (n number of counts knownY) 2757 for (uint i = 1; i < args[0].count() + 1; ++i) 2758 knownX.append(i); 2759 } else { 2760 // check constraints / TODO if 2d array, then we must check dimension row&col? 2761 if (args[0].count() != args[1].count()) 2762 return Value::errorNUM(); 2763 2764 // copy array to list 2765 func_array_helper(args[1], calc, knownX, knownXcount); 2766 } 2767 2768 // 2769 // newX 2770 // 2771 if (args[2].isEmpty()) { 2772 for (uint i = 1; i < args[0].count() + 1; ++i) 2773 newX.append(i); 2774 } else { 2775 // copy array to list 2776 func_array_helper(args[2], calc, newX, newXcount); 2777 } 2778 2779 // create the resulting matrix 2780 Value res(Value::Array); 2781 2782 // arrays for slope und intercept 2783 Value Y(Value::Array); 2784 Value X(Value::Array); 2785 2786 Value sumXX(0.0); // stores sum of xi*xi 2787 Value sumYX(0.0); // stores sum of yi*xi 2788 2789 // copy data in arrays 2790 for (uint i = 0; i < args[0].count(); ++i) { 2791 X.setElement(i, 0, Value((double)knownX[i])); 2792 sumXX = calc->add(sumXX, calc->mul(Value((double)knownX[i]), Value((double)knownX[i]))); 2793 } 2794 2795 for (uint i = 0; i < args[0].count(); ++i) { 2796 Y.setElement(i, 0, Value(args[0].element(i))); 2797 // sum yi*xi 2798 sumYX = calc->add(sumYX, calc->mul(Value(args[0].element(i)), Value((double)knownX[i]))); 2799 } 2800 2801 // create parameter for func_slope and func_intercept calls 2802 valVector param; 2803 2804 param.append(Y); 2805 param.append(X); 2806 2807 // a1 = [xy]/[x*x] 2808 Value a1 = calc->div(sumYX, sumXX); 2809 2810 // v2 = INTERCEPT = b 2811 Value v2 = func_intercept(param, calc, 0); // v2 is const, we only need to calc it once 2812 2813 // fill array up with values 2814 for (uint i = 0; i < args[2].count(); ++i) { 2815 Value trend; 2816 Value v1; 2817 2818 if (withOffset) { 2819 v1 = calc->mul(func_slope(param, calc, 0), Value(newX[i])); 2820 trend = Value(calc->add(v1 , v2)); 2821 } else { 2822 // b=0 2823 // x*a1 2824 trend = calc->mul(a1, Value(newX[i])); 2825 } 2826 2827 // set value in res array 2828 res.setElement(i, 0, trend); 2829 } 2830 2831 return (res); // return array 2832 } 2833 2834 // 2835 // Function: trimmean 2836 // 2837 Value func_trimmean(valVector args, ValueCalc *calc, FuncExtra *) 2838 { 2839 Value dataSet = args[0]; 2840 Value cutOffFrac = args[1]; 2841 2842 // constrains 0 <= cutOffFrac < 1 2843 if (calc->lower(cutOffFrac, Value(0)) || !calc->lower(cutOffFrac, Value(1))) 2844 return Value::errorVALUE(); 2845 2846 // cutOff = floor( n*cutOffFrac/2) 2847 int cutOff = floor(calc->div(calc->mul(cutOffFrac , Value((int)dataSet.count())), 2).asFloat()); 2848 2849 double res = 0.0; 2850 2851 // sort parameter into QList array 2852 List array; 2853 int valCount = 0; // stores the number of values in array 2854 2855 func_array_helper(args[0], calc, array, valCount); 2856 2857 if (valCount == 0) 2858 return Value::errorVALUE(); 2859 2860 std::sort(array.begin(), array.end()); 2861 2862 for (int i = cutOff; i < valCount - cutOff ; ++i) 2863 res += array[i]; 2864 2865 res /= (valCount - 2 * cutOff); 2866 2867 return Value(res); 2868 } 2869 2870 // 2871 // Function TTEST 2872 // 2873 Value func_ttest(valVector args, ValueCalc* calc, FuncExtra*) 2874 { 2875 Value x = args[0]; 2876 Value y = args[1]; 2877 int mode = calc->conv()->asInteger(args[2]).asInteger(); 2878 int type = calc->conv()->asInteger(args[3]).asInteger(); 2879 2880 int numX = calc->count(x); 2881 int numY = calc->count(y); 2882 2883 // check mode parameter 2884 if (mode < 1 || mode > 2) 2885 return Value::errorVALUE(); 2886 // check type parameter 2887 if (type < 1 || type > 3) 2888 return Value::errorVALUE(); 2889 // check amount of numbers in sequences 2890 if (numX < 2 || numY < 2 || (type == 1 && numX != numY)) 2891 return Value::errorVALUE(); 2892 2893 Value t; 2894 Value dof; 2895 if (type == 1) { 2896 // paired 2897 dof = calc->sub(Value(numX), 1); 2898 2899 Value mean; 2900 calc->twoArrayWalk(x, y, mean, tawSumxmy); 2901 mean = calc->div(mean, numX); 2902 2903 Value sigma; 2904 calc->twoArrayWalk(x, y, sigma, tawSumxmy2); 2905 sigma = calc->sqrt(calc->sub(calc->div(sigma, numX), calc->sqr(mean))); 2906 2907 t = calc->div(calc->mul(mean, calc->sqrt(dof)), sigma); 2908 } else if (type == 2) { 2909 // independent, equal variances (revised by Jon Cooper) 2910 dof = calc->sub(calc->add(Value(numX), Value(numY)), 2); 2911 2912 Value avgX = calc->avg(x); 2913 Value avgY = calc->avg(y); 2914 Value varX, varY; // summed dev-squares 2915 calc->arrayWalk(x, varX, calc->awFunc("devsq"), avgX); 2916 calc->arrayWalk(y, varY, calc->awFunc("devsq"), avgY); 2917 2918 Value numerator = calc->sub(avgX, avgY); 2919 2920 Value pooled_variance = calc->add(varX, varY); 2921 pooled_variance = calc->div(pooled_variance, dof); 2922 2923 Value denominator = calc->add(calc->div(pooled_variance,Value(numX)), calc->div(pooled_variance,Value(numY))); 2924 denominator = calc->sqrt(denominator); 2925 2926 t = calc->div(numerator, denominator); 2927 } else { 2928 // independent, unequal variances (revised by Jon Cooper) 2929 2930 Value avgX = calc->avg(x); 2931 Value avgY = calc->avg(y); 2932 Value varX, varY; 2933 calc->arrayWalk(x, varX, calc->awFunc("devsq"), avgX); 2934 calc->arrayWalk(y, varY, calc->awFunc("devsq"), avgY); 2935 2936 varX = calc->div(varX,calc->sub(Value(numX),1)); 2937 varY = calc->div(varY,calc->sub(Value(numY),1)); 2938 2939 Value numerator = calc->sub(avgX, avgY); 2940 Value denominator = calc->add(calc->div(varX,Value(numX)), calc->div(varY,Value(numY))); 2941 denominator = calc->sqrt(denominator); 2942 t = calc->div(numerator, denominator); 2943 2944 numerator = calc->add(calc->div(varX,Value(numX)),calc->div(varY,Value(numY))); 2945 numerator = calc->pow(numerator,2); 2946 2947 Value denominator1 = calc->div(calc->pow(calc->div(varX,Value(numX)),2),calc->sub(Value(numX),1)); 2948 Value denominator2 = calc->div(calc->pow(calc->div(varY,Value(numY)),2),calc->sub(Value(numY),1)); 2949 denominator = calc->add(denominator1,denominator2); 2950 dof = calc->div(numerator,denominator); 2951 2952 } 2953 2954 valVector tmp(3); 2955 tmp.insert(0, t); 2956 tmp.insert(1, dof); 2957 tmp.insert(2, Value(mode)); 2958 return func_tdist(tmp, calc, 0); 2959 } 2960 2961 // 2962 // Function: variance 2963 // 2964 Value func_variance(valVector args, ValueCalc *calc, FuncExtra *) 2965 { 2966 int count = calc->count(args, false); 2967 if (count < 2) 2968 return Value::errorVALUE(); 2969 2970 Value result = func_devsq(args, calc, 0); 2971 return calc->div(result, count - 1); 2972 } 2973 2974 // 2975 // Function: vara 2976 // 2977 Value func_variancea(valVector args, ValueCalc *calc, FuncExtra *) 2978 { 2979 int count = calc->count(args); 2980 if (count < 2) 2981 return Value::errorVALUE(); 2982 2983 Value result = func_devsqa(args, calc, 0); 2984 return calc->div(result, count - 1); 2985 } 2986 2987 // 2988 // Function: varp 2989 // 2990 Value func_variancep(valVector args, ValueCalc *calc, FuncExtra *) 2991 { 2992 int count = calc->count(args, false); 2993 if (count == 0) 2994 return Value::errorVALUE(); 2995 2996 Value result = func_devsq(args, calc, 0); 2997 return calc->div(result, count); 2998 } 2999 3000 // 3001 // Function: varpa 3002 // 3003 Value func_variancepa(valVector args, ValueCalc *calc, FuncExtra *) 3004 { 3005 int count = calc->count(args); 3006 if (count == 0) 3007 return Value::errorVALUE(); 3008 3009 Value result = func_devsqa(args, calc, 0); 3010 return calc->div(result, count); 3011 } 3012 3013 // 3014 // Function: weibull 3015 // 3016 // returns the Weibull distribution 3017 // 3018 Value func_weibull(valVector args, ValueCalc *calc, FuncExtra *) 3019 { 3020 Value x = args[0]; 3021 Value alpha = args[1]; 3022 Value beta = args[2]; 3023 Value kum = args[3]; 3024 3025 Value result; 3026 3027 if ((!calc->greater(alpha, 0.0)) || (!calc->greater(beta, 0.0)) || 3028 calc->lower(x, 0.0)) 3029 return Value::errorVALUE(); 3030 3031 // ex = exp (-pow (x / beta, alpha)) 3032 Value ex; 3033 ex = calc->exp(calc->mul(calc->pow(calc->div(x, beta), alpha), -1)); 3034 if (calc->isZero(kum)) { // density 3035 // result = alpha / pow(beta,alpha) * pow(x,alpha-1.0) * ex 3036 result = calc->div(alpha, calc->pow(beta, alpha)); 3037 result = calc->mul(result, calc->mul(calc->pow(x, 3038 calc->sub(alpha, 1)), ex)); 3039 } else // distribution 3040 result = calc->sub(1.0, ex); 3041 3042 return result; 3043 } 3044 3045 // 3046 // Function ZTEST 3047 // 3048 Value func_ztest(valVector args, ValueCalc* calc, FuncExtra*) 3049 { 3050 int number = calc->count(args[0]); 3051 3052 if (number < 2) 3053 return Value::errorVALUE(); 3054 3055 // standard deviation is optional 3056 Value sigma = (args.count() > 2) ? args[2] : calc->stddev(args[0], false); 3057 // z = Ex - mu / sigma * sqrt(N) 3058 Value z = calc->div(calc->mul(calc->sub(calc->avg(args[0]), args[1]), calc->sqrt(Value(number))), sigma); 3059 // result = 1 - ( NORMDIST(-abs(z),0,1,1) + ( 1 - NORMDIST(abs(z),0,1,1) ) ) 3060 return calc->mul(Value(2.0), calc->gauss(calc->abs(z))); 3061 } 3062 3063 #include "statistical.moc"