LibreOffice Module sc (master) 1
op_statistical_helpers.hxx
Go to the documentation of this file.
1/* -*- Mode: C++; tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 4 -*- */
2/*
3 * This file is part of the LibreOffice project.
4 *
5 * This Source Code Form is subject to the terms of the Mozilla Public
6 * License, v. 2.0. If a copy of the MPL was not distributed with this
7 * file, You can obtain one at http://mozilla.org/MPL/2.0/.
8 */
9
10#pragma once
11
12const char MinDecl[] = "#define Min 2.22507e-308\n";
13const char fBigInvDecl[] ="#define fBigInv 2.22045e-016\n";
14const char fMachEpsDecl[] ="#define fMachEps 2.22045e-016\n";
15const char fLogDblMaxDecl[] ="#define fLogDblMax log(1.79769e+308)\n";
16const char fHalfMachEpsDecl[] ="#define fHalfMachEps 0.5*2.22045e-016\n";
17const char fMaxGammaArgumentDecl[] =
18"#define fMaxGammaArgument 171.624376956302\n";
19const char GetValueDecl[] =
20"double GetValue( double x,double fp,double fDF );\n";
21const char GetValue[] =
22"double GetValue( double x,double fp,double fDF )\n"
23"{\n"
24" return fp - 2 * GetTDist(x, fDF);\n"
25"}\n";
26const char fmin_countDecl[] = "double fmin_count(double a, double b, __private int *p);\n";
27const char fmin_count[] =
28"double fmin_count(double a, double b, __private int *p) {\n"
29" double result = fmin(a, b);\n"
30" bool t = isnan(result);\n"
31" (*p) += t?0:1;\n"
32" return result;\n"
33"}\n";
34const char fmax_countDecl[] = "double fmax_count(double a, double b, __private int *p);\n";
35const char fmax_count[] =
36"double fmax_count(double a, double b, __private int *p) {\n"
37" double result = fmax(a, b);\n"
38" bool t = isnan(result);\n"
39" (*p) += t?0:1;\n"
40" return result;\n"
41"}\n";
42const char fsum_countDecl[] = "double fsum_count(double a, double b, __private int *p);\n";
43const char fsum_count[] =
44"double fsum_count(double a, double b, __private int *p) {\n"
45" bool t = isnan(a);\n"
46" (*p) += t?0:1;\n"
47" return t?b:a+b;\n"
48"}\n";
49const char GetGammaSeriesDecl[] =
50"double GetGammaSeries( double fA, double fX );\n";
51const char GetGammaSeries[] =
52"double GetGammaSeries( double fA, double fX )\n"
53"{\n"
54" double fDenomfactor = fA;\n"
55" double fSummand = 1.0/fA;\n"
56" double fSum = fSummand;\n"
57" int nCount=1;\n"
58" do\n"
59" {\n"
60" fDenomfactor = fDenomfactor + 1.0;\n"
61" fSummand = fSummand * fX/fDenomfactor;\n"
62" fSum = fSum + fSummand;\n"
63" nCount = nCount+1;\n"
64" } while ( fSummand/fSum > fHalfMachEps && nCount<=10000);\n"
65" if (nCount>10000)\n"
66" {\n"
67" }\n"
68" return fSum;\n"
69"}\n";
70const char GetGammaContFractionDecl[] = "double GetGammaContFraction( double "
71"fA, double fX );\n";
72const char GetGammaContFraction[] =
73"double GetGammaContFraction( double fA, double fX )\n"
74"{\n"
75" double fBig = 1.0/fBigInv;\n"
76" double fCount = 0.0;\n"
77" double fNum = 0.0;\n"
78" double fY = 1.0 - fA;\n"
79" double fDenom = fX + 2.0-fA;\n"
80" double fPk = 0.0;\n"
81" double fPkm1 = fX + 1.0;\n"
82" double fPkm2 = 1.0;\n"
83" double fQk = 1.0;\n"
84" double fQkm1 = fDenom * fX;\n"
85" double fQkm2 = fX;\n"
86" double fApprox = fPkm1/fQkm1;\n"
87" bool bFinished = false;\n"
88" double fR = 0.0;\n"
89" do\n"
90" {\n"
91" fCount = fCount +1.0;\n"
92" fY = fY+ 1.0;\n"
93" fNum = fY * fCount;\n"
94" fDenom = fDenom +2.0;\n"
95" fPk = fPkm1 * fDenom - fPkm2 * fNum;\n"
96" fQk = fQkm1 * fDenom - fQkm2 * fNum;\n"
97" if (fQk != 0.0)\n"
98" {\n"
99" fR = fPk/fQk;\n"
100" bFinished = (fabs( (fApprox - fR)/fR ) <= fHalfMachEps);\n"
101" fApprox = fR;\n"
102" }\n"
103" fPkm2 = fPkm1;\n"
104" fPkm1 = fPk;\n"
105" fQkm2 = fQkm1;\n"
106" fQkm1 = fQk;\n"
107" if (fabs(fPk) > fBig)\n"
108" {\n"
109" fPkm2 = fPkm2 * fBigInv;\n"
110" fPkm1 = fPkm1 * fBigInv;\n"
111" fQkm2 = fQkm2 * fBigInv;\n"
112" fQkm1 = fQkm1 * fBigInv;\n"
113" }\n"
114" } while (!bFinished && fCount<10000);\n"
115" if (!bFinished)\n"
116" {\n"
117" }\n"
118" return fApprox;\n"
119"}\n";
120const char GetLowRegIGammaDecl[] = "double GetLowRegIGamma( double "
121"fA, double fX );\n";
122const char GetLowRegIGamma[] =
123"double GetLowRegIGamma( double fA, double fX )\n"
124"{\n"
125" double fLnFactor = fA * log(fX) - fX - lgamma(fA);\n"
126" double fFactor = exp(fLnFactor);\n"
127" if (fX>fA+1.0) \n"
128" return 1.0 - fFactor * GetGammaContFraction(fA,fX);\n"
129" else\n"
130" return fFactor * GetGammaSeries(fA,fX);\n"
131"}\n";
132const char GetGammaDistDecl[] = "double GetGammaDist( double fX, double "
133"fAlpha, double fLambda );\n";
134const char GetGammaDist[] =
135"double GetGammaDist( double fX, double fAlpha, double fLambda )\n"
136"{\n"
137" if (fX <= 0.0)\n"
138" return 0.0;\n"
139" else\n"
140" return GetLowRegIGamma( fAlpha, fX / fLambda);\n"
141"}\n";
142const char GetGammaDistPDFDecl[] = "double GetGammaDistPDF( double fX"
143", double fAlpha, double fLambda );\n";
144const char GetGammaDistPDF[] =
145"double GetGammaDistPDF( double fX, double fAlpha, double fLambda )\n"
146"{\n"
147" if (fX < 0.0)\n"
148" return 0.0;\n"
149" else if (fX == 0)\n"
150" {\n"
151" if (fAlpha < 1.0)\n"
152" {\n"
153" return HUGE_VAL;\n"
154" }\n"
155" else if (fAlpha == 1)\n"
156" {\n"
157" return (1.0 / fLambda);\n"
158" }\n"
159" else\n"
160" {\n"
161" return 0.0;\n"
162" }\n"
163" }\n"
164" else\n"
165" {\n"
166" double fXr = fX / fLambda;\n"
167" if (fXr > 1.0)\n"
168" {\n"
169" if (log(fXr) * (fAlpha-1.0) < fLogDblMax &&"
170"fAlpha < fMaxGammaArgument)\n"
171" {\n"
172" return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
173"fLambda / tgamma(fAlpha);\n"
174" }\n"
175" else\n"
176" {\n"
177" return exp( (fAlpha-1.0) * log(fXr) - fXr - "
178"log(fLambda) - lgamma(fAlpha));\n"
179" }\n"
180" }\n"
181" else\n"
182" {\n"
183" if (fAlpha<fMaxGammaArgument)\n"
184" {\n"
185" return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
186"fLambda / tgamma(fAlpha);\n"
187" }\n"
188" else\n"
189" {\n"
190" return pow( fXr, fAlpha-1.0) * exp(-fXr) / "
191"fLambda / exp( lgamma(fAlpha));\n"
192" }\n"
193" }\n"
194" }\n"
195"}\n";
196const char GetBetaDistDecl[] =
197"double GetBetaDist(double fXin, double fAlpha, double fBeta);\n";
198const char GetBetaDist[] =
199"double GetBetaDist(double fXin, double fAlpha, double fBeta)\n"
200"{\n"
201" if (fXin <= 0.0)\n"
202" return 0.0;\n"
203" if (fXin >= 1.0)\n"
204" return 1.0;\n"
205" if (fBeta == 1.0)\n"
206" return pow(fXin, fAlpha);\n"
207" if (fAlpha == 1.0)\n"
208" return -expm1(fBeta * log1p(-fXin));\n"
209" double fResult;\n"
210" double fY = (0.5-fXin)+0.5;\n"
211" double flnY = log1p(-fXin);\n"
212" double fX = fXin;\n"
213" double flnX = log(fXin);\n"
214" double fA = fAlpha;\n"
215" double fB = fBeta;\n"
216" bool bReflect = fXin > fAlpha/(fAlpha+fBeta);\n"
217" if (bReflect)\n"
218" {\n"
219" fA = fBeta;\n"
220" fB = fAlpha;\n"
221" fX = fY;\n"
222" fY = fXin;\n"
223" flnX = flnY;\n"
224" flnY = log(fXin);\n"
225" }\n"
226" fResult = lcl_GetBetaHelperContFrac(fX,fA,fB)/fA;\n"
227" double fP = fA/(fA+fB);\n"
228" double fQ = fB/(fA+fB);\n"
229" if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
230" fResult *= GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
231" else\n"
232" fResult *= pow(exp(1.0),(fA*flnX + fB*flnY - GetLogBeta(fA,fB)));\n"
233" if (bReflect)\n"
234" fResult = 0.5 - fResult + 0.5;\n"
235" if (fResult > 1.0)\n"
236" fResult = 1.0;\n"
237" if (fResult < 0.0)\n"
238" fResult = 0.0;\n"
239" return fResult;\n"
240"}\n";
241
242const char GetFDistDecl[] =
243 "double GetFDist(double x, double fF1, double fF2);\n";
244const char GetFDist[] =
245"double GetFDist(double x, double fF1, double fF2)\n"
246"{\n"
247" double arg = fF2/(fF2+fF1*x);\n"
248" double alpha = fF2/2.0;\n"
249" double beta = fF1/2.0;\n"
250" return (GetBetaDist(arg, alpha, beta));\n"
251"}\n";
252const char GetGammaInvValueDecl[] = "double"
253" GetGammaInvValue(double fAlpha,double fBeta,double fX1 );\n";
254const char GetGammaInvValue[] =
255"double GetGammaInvValue(double fAlpha,double fBeta,double fX1 )\n"
256"{\n"
257" if (fX1 <= 0.0)\n"
258" return 0.0;\n"
259" else\n"
260" {\n"
261" double fX=fX1/fBeta;\n"
262" double fLnFactor = fAlpha * log(fX) - fX - lgamma(fAlpha);\n"
263" double fFactor = exp(fLnFactor);\n"
264" if (fX>fAlpha+1.0)\n"
265" return 1.0 - fFactor * GetGammaContFraction(fAlpha,fX);\n"
266" else\n"
267" return fFactor * GetGammaSeries(fAlpha,fX);\n"
268" }\n"
269"}\n";
270const char GetFInvValueDecl[] = "double GetFInvValue(double fF1,double fF2"
271",double fX1 );";
272const char GetFInvValue[] =
273"double GetFInvValue(double fF1,double fF2,double fX1 )\n"
274"{\n"
275" double arg = fF2/(fF2+fF1*fX1);\n"
276" double alpha = fF2/2.0;\n"
277" double beta = fF1/2.0;\n"
278" double fXin,fAlpha,fBeta;\n"
279" fXin=arg;fAlpha=alpha;fBeta=beta;\n"
280" if (fXin <= 0.0)\n"
281" return 0.0;\n"
282" if (fXin >= 1.0)\n"
283" return 1.0;\n"
284" if (fBeta == 1.0)\n"
285" return pow(fXin, fAlpha);\n"
286" if (fAlpha == 1.0)\n"
287" return -expm1(fBeta * log1p(-fXin));\n"
288" double fResult;\n"
289" double fY = (0.5-fXin)+0.5;\n"
290" double flnY = log1p(-fXin);\n"
291" double fX = fXin;\n"
292" double flnX = log(fXin);\n"
293" double fA = fAlpha;\n"
294" double fB = fBeta;\n"
295" bool bReflect = fXin > fAlpha/(fAlpha+fBeta);\n"
296" if (bReflect)\n"
297" {\n"
298" fA = fBeta;\n"
299" fB = fAlpha;\n"
300" fX = fY;\n"
301" fY = fXin;\n"
302" flnX = flnY;\n"
303" flnY = log(fXin);\n"
304" }\n"
305" fResult = lcl_GetBetaHelperContFrac(fX,fA,fB);\n"
306" fResult = fResult/fA;\n"
307" double fP = fA/(fA+fB);\n"
308" double fQ = fB/(fA+fB);\n"
309" double fTemp;\n"
310" if (fA > 1.0 && fB > 1.0 && fP < 0.97 && fQ < 0.97)\n"
311" fTemp = GetBetaDistPDF(fX,fA,fB)*fX*fY;\n"
312" else\n"
313" fTemp = exp(fA*flnX + fB*flnY - GetLogBeta(fA,fB));\n"
314" fResult *= fTemp;\n"
315" if (bReflect)\n"
316" fResult = 0.5 - fResult + 0.5;\n"
317" if (fResult > 1.0)\n"
318" fResult = 1.0;\n"
319" if (fResult < 0.0)\n"
320" fResult = 0.0;\n"
321" return fResult;\n"
322"}\n";
323const char GetBinomDistPMFDecl[] =
324 "double GetBinomDistPMF(double x, double n, double p);";
325const char GetBinomDistPMF[] =
326"double GetBinomDistPMF(double x, double n, double p)\n"
327"{\n"
328" double q = (0.5 - p) + 0.5;\n"
329" double fFactor = pow(q, n);\n"
330" if (fFactor <= Min)\n"
331" {\n"
332" fFactor = pow(p, n);\n"
333" if (fFactor <= Min)\n"
334" return GetBetaDistPDF(p, x + 1.0, n - x + 1.0)/(n + 1.0);\n"
335" else\n"
336" {\n"
337" uint max = (uint)(n - x);\n"
338" for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
339" fFactor *= (n - i)/(double)(i + 1)*q/p;\n"
340" return fFactor;\n"
341" }\n"
342" }\n"
343" else\n"
344" {\n"
345" uint max = (uint)x;\n"
346" for (uint i = 0; i < max && fFactor > 0.0; ++i)\n"
347" fFactor *= (n - i)/(double)(i + 1)*p/q;\n"
348" return fFactor;\n"
349" }\n"
350"}\n";
351
352const char lcl_GetBinomDistRangeDecl[] =
353 "double lcl_GetBinomDistRange(double n, \n"
354"double xs, double xe, double fFactor, double p, double q);";
355const char lcl_GetBinomDistRange[] =
356"double lcl_GetBinomDistRange(double n, double xs, double xe,\n"
357" double fFactor, double p, double q)\n"
358"{\n"
359" uint i;\n"
360" double fSum;\n"
361" uint nXs = (uint)xs;\n"
362" for (i = 1; i <= nXs && fFactor > 0.0; ++i)\n"
363" fFactor *= (n - i + 1)/(double)(i)*p/q;\n"
364" fSum = fFactor;\n"
365" uint nXe =(uint)xe;\n"
366" for (i = nXs + 1; i <= nXe && fFactor > 0.0; ++i)\n"
367" {\n"
368" fFactor *= (n - i + 1)/(double)(i)*p/q;\n"
369" fSum += fFactor;\n"
370" }\n"
371" return (fSum > 1.0) ? 1.0 : fSum;\n"
372"}\n";
373
374const char GetLogGammaDecl[] = "double GetLogGamma(double fZ);\n";
375const char GetLogGamma[] =
376"double GetLogGamma(double fZ)\n"
377"{\n"
378" if (fZ >= fMaxGammaArgument)\n"
379" return lcl_GetLogGammaHelper(fZ);\n"
380" if (fZ >= 1.0)\n"
381" return log(lcl_GetGammaHelper(fZ));\n"
382" if (fZ >= 0.5)\n"
383" return log( lcl_GetGammaHelper(fZ+1) / fZ);\n"
384" return lcl_GetLogGammaHelper(fZ+2) - log(fZ+1) - log(fZ);\n"
385"}\n";
386
387const char GetChiDistDecl[] = "double GetChiDist(double fX, double fDF);\n";
388const char GetChiDist[] =
389"double GetChiDist(double fX, double fDF)\n"
390"{\n"
391" if (fX <= 0.0)\n"
392" return 1.0;\n"
393" else\n"
394" return GetUpRegIGamma( fDF/2.0, fX/2.0);\n"
395"}\n";
396
397const char GetChiSqDistCDFDecl[] =
398"double GetChiSqDistCDF(double fX, double fDF);\n";
399const char GetChiSqDistCDF[] =
400"double GetChiSqDistCDF(double fX, double fDF)\n"
401"{\n"
402" if (fX <= 0.0)\n"
403" return 0.0;"
404" else\n"
405" return GetLowRegIGamma( fDF/2.0, fX/2.0);\n"
406"}\n";
407
408const char GetChiSqDistPDFDecl[] =
409"double GetChiSqDistPDF(double fX, double fDF);\n";
410const char GetChiSqDistPDF[] =
411"double GetChiSqDistPDF(double fX, double fDF)\n"
412"{\n"
413" double fValue;\n"
414" if (fX <= 0.0)\n"
415" return 0.0;\n"
416" if (fDF*fX > 1391000.0)\n"
417" {\n"
418" fValue = exp((0.5*fDF - 1) * log(fX*0.5) - 0.5 * fX - log(2.0)"
419" - lgamma(0.5*fDF));\n"
420" }\n"
421" else\n"
422" {\n"
423" double fCount;\n"
424" if (fmod(fDF,2.0)<0.5)\n"
425" {\n"
426" fValue = 0.5;\n"
427" fCount = 2.0;\n"
428" }\n"
429" else\n"
430" {\n"
431" fValue = 1.0/sqrt(fX*2*M_PI);\n"
432" fCount = 1.0;\n"
433" }\n"
434" while ( fCount < fDF)\n"
435" {\n"
436" fValue *= (fX / fCount);\n"
437" fCount += 2.0;\n"
438" }\n"
439" if (fX>=1425.0)\n"
440" fValue = exp(log(fValue)-fX/2);\n"
441" else\n"
442" fValue *= exp(-fX/2);\n"
443" }\n"
444" return fValue;\n"
445"}\n";
446
447const char lcl_IterateInverseBetaInvDecl[] =
448"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
449" double fBeta, double fAx, double fBx, bool *rConvError );\n";
450const char lcl_IterateInverseBetaInv[] =
451"static double lcl_IterateInverseBetaInv(double fp, double fAlpha, \n"
452" double fBeta, double fAx, double fBx, bool *rConvError )\n"
453"{\n"
454" *rConvError = false;\n"
455" double fYEps = 1.0E-307;\n"
456" double fXEps = fMachEps;\n"
457" if(!(fAx < fBx))\n"
458" {\n"
459" //print error\n"
460" }\n"
461" double fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
462" double fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
463" double fTemp;\n"
464" unsigned short nCount;\n"
465" for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
466" nCount++)\n"
467" {\n"
468" if (fabs(fAy) <= fabs(fBy))\n"
469" {\n"
470" fTemp = fAx;\n"
471" fAx += 2.0 * (fAx - fBx);\n"
472" if (fAx < 0.0)\n"
473" fAx = 0.0;\n"
474" fBx = fTemp;\n"
475" fBy = fAy;\n"
476" fAy = fp - GetBetaDist(fAx, fAlpha, fBeta);\n"
477" }\n"
478" else\n"
479" {\n"
480" fTemp = fBx;\n"
481" fBx += 2.0 * (fBx - fAx);\n"
482" fAx = fTemp;\n"
483" fAy = fBy;\n"
484" fBy = fp - GetBetaDist(fBx, fAlpha, fBeta);\n"
485" }\n"
486" }\n"
487" if (fAy == 0.0)\n"
488" return fAx;\n"
489" if (fBy == 0.0)\n"
490" return fBx;\n"
491" if (!lcl_HasChangeOfSign( fAy, fBy))\n"
492" {\n"
493" *rConvError = true;\n"
494" return 0.0;\n"
495" }\n"
496" double fPx = fAx;\n"
497" double fPy = fAy;\n"
498" double fQx = fBx;\n"
499" double fQy = fBy;\n"
500" double fRx = fAx;\n"
501" double fRy = fAy;\n"
502" double fSx = 0.5 * (fAx + fBx);\n"
503" bool bHasToInterpolate = true;\n"
504" nCount = 0;\n"
505" while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
506" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
507" {\n"
508" if (bHasToInterpolate)\n"
509" {\n"
510" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
511" {\n"
512" fSx = fPx*fRy*fQy/(fRy-fPy)/(fQy-fPy)\n"
513" + fRx*fQy*fPy/(fQy-fRy)/(fPy-fRy)\n"
514" + fQx*fPy*fRy/(fPy-fQy)/(fRy-fQy);\n"
515" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
516" }\n"
517" else\n"
518" bHasToInterpolate = false;\n"
519" }\n"
520" if(!bHasToInterpolate)\n"
521" {\n"
522" fSx = 0.5 * (fAx + fBx);\n"
523" fPx = fAx; fPy = fAy;\n"
524" fQx = fBx; fQy = fBy;\n"
525" bHasToInterpolate = true;\n"
526" }\n"
527" fPx = fQx; fQx = fRx; fRx = fSx;\n"
528" fPy = fQy; fQy = fRy; fRy = fp - GetBetaDist(fSx, fAlpha, fBeta);\n"
529" if (lcl_HasChangeOfSign( fAy, fRy))\n"
530" {\n"
531" fBx = fRx; fBy = fRy;\n"
532" }\n"
533" else\n"
534" {\n"
535" fAx = fRx; fAy = fRy;\n"
536" }\n"
537" bHasToInterpolate = bHasToInterpolate && (fabs(fRy) *"
538" 2.0 <= fabs(fQy));\n"
539" ++nCount;\n"
540" }\n"
541" return fRx;\n"
542"}\n";
543
544const char lcl_IterateInverseChiInvDecl[] =
545"static double lcl_IterateInverseChiInv"
546"(double fp, double fdf, double fAx, double fBx, bool *rConvError);\n";
547const char lcl_IterateInverseChiInv[] =
548"static double lcl_IterateInverseChiInv"
549"(double fp, double fdf, double fAx, double fBx, bool *rConvError)\n"
550"{\n"
551" *rConvError = false;\n"
552" double fYEps = 1.0E-307;\n"
553" double fXEps = fMachEps;\n"
554" if(!(fAx < fBx))\n"
555" {\n"
556" //print error\n"
557" }"
558" double fAy = fp - GetChiDist(fAx, fdf);\n"
559" double fBy = fp - GetChiDist(fBx, fdf);\n"
560" double fTemp;\n"
561" unsigned short nCount;\n"
562" for (nCount = 0; nCount < 1000 && "
563"!lcl_HasChangeOfSign(fAy,fBy); nCount++)\n"
564" {\n"
565" if (fabs(fAy) <= fabs(fBy))\n"
566" {\n"
567" fTemp = fAx;\n"
568" fAx += 2.0 * (fAx - fBx);\n"
569" if (fAx < 0.0)\n"
570" fAx = 0.0;\n"
571" fBx = fTemp;\n"
572" fBy = fAy;\n"
573" fAy = fp - GetChiDist(fAx, fdf);\n"
574" }\n"
575" else\n"
576" {\n"
577" fTemp = fBx;\n"
578" fBx += 2.0 * (fBx - fAx);\n"
579" fAx = fTemp;\n"
580" fAy = fBy;\n"
581" fBy = fp - GetChiDist(fBx, fdf);\n"
582" }\n"
583" }\n"
584" if (fAy == 0.0)\n"
585" return fAx;\n"
586" if (fBy == 0.0)\n"
587" return fBx;\n"
588" if (!lcl_HasChangeOfSign( fAy, fBy))\n"
589" {\n"
590" *rConvError = true;\n"
591" return 0.0;\n"
592" }\n"
593" double fPx = fAx;\n"
594" double fPy = fAy;\n"
595" double fQx = fBx;\n"
596" double fQy = fBy;\n"
597" double fRx = fAx;\n"
598" double fRy = fAy;\n"
599" double fSx = 0.5 * (fAx + fBx);\n"
600" bool bHasToInterpolate = true;\n"
601" nCount = 0;\n"
602" while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
603" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
604" {\n"
605" if (bHasToInterpolate)\n"
606" {\n"
607" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
608" {\n"
609" fSx = fPx * fRy * fQy/(fRy-fPy)/(fQy-fPy)\n"
610" + fRx * fQy * fPy/(fQy-fRy)/(fPy-fRy)\n"
611" + fQx * fPy * fRy/(fPy-fQy)/(fRy-fQy);\n"
612" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
613" }\n"
614" else\n"
615" bHasToInterpolate = false;\n"
616" }\n"
617" if(!bHasToInterpolate)\n"
618" {\n"
619" fSx = 0.5 * (fAx + fBx);\n"
620" fPx = fAx; fPy = fAy;\n"
621" fQx = fBx; fQy = fBy;\n"
622" bHasToInterpolate = true;\n"
623" }\n"
624" fPx = fQx; fQx = fRx; fRx = fSx;\n"
625" fPy = fQy; fQy = fRy; fRy = fp - GetChiDist(fSx, fdf);\n"
626" if (lcl_HasChangeOfSign( fAy, fRy))\n"
627" {\n"
628" fBx = fRx; fBy = fRy;\n"
629" }\n"
630" else\n"
631" {\n"
632" fAx = fRx; fAy = fRy;\n"
633" }\n"
634" bHasToInterpolate = bHasToInterpolate && (fabs(fRy)"
635" * 2.0 <= fabs(fQy));\n"
636" ++nCount;\n"
637" }\n"
638" return fRx;\n"
639"}\n";
640
641const char lcl_IterateInverseChiSQInvDecl[] =
642"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
643" double fAx, double fBx, bool *rConvError );\n";
644const char lcl_IterateInverseChiSQInv[] =
645"static double lcl_IterateInverseChiSQInv( double fp, double fdf, \n"
646" double fAx, double fBx, bool *rConvError )\n"
647"{\n"
648" *rConvError = false;\n"
649" double fYEps = 1.0E-307;\n"
650" double fXEps = fMachEps;\n"
651
652" if(!(fAx < fBx))\n"
653" {\n"
654" //print error\n"
655" }\n"
656" double fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
657" double fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
658" double fTemp;\n"
659" unsigned short nCount;\n"
660" for (nCount = 0; nCount < 1000 && !lcl_HasChangeOfSign(fAy,fBy);"
661" nCount++)\n"
662" {\n"
663" if (fabs(fAy) <= fabs(fBy))\n"
664" {\n"
665" fTemp = fAx;\n"
666" fAx += 2.0 * (fAx - fBx);\n"
667" if (fAx < 0.0)\n"
668" fAx = 0.0;\n"
669" fBx = fTemp;\n"
670" fBy = fAy;\n"
671" fAy = fp - GetChiSqDistCDF(fAx, fdf);\n"
672" }\n"
673" else\n"
674" {\n"
675" fTemp = fBx;\n"
676" fBx += 2.0 * (fBx - fAx);\n"
677" fAx = fTemp;\n"
678" fAy = fBy;\n"
679" fBy = fp - GetChiSqDistCDF(fBx, fdf);\n"
680" }\n"
681" }\n"
682" if (fAy == 0.0)\n"
683" return fAx;\n"
684" if (fBy == 0.0)\n"
685" return fBx;\n"
686" if (!lcl_HasChangeOfSign( fAy, fBy))\n"
687" {\n"
688" *rConvError = true;\n"
689" return 0.0;\n"
690" }\n"
691" double fPx = fAx;\n"
692" double fPy = fAy;\n"
693" double fQx = fBx;\n"
694" double fQy = fBy;\n"
695" double fRx = fAx;\n"
696" double fRy = fAy;\n"
697" double fSx = 0.5 * (fAx + fBx);\n"
698" bool bHasToInterpolate = true;\n"
699" nCount = 0;\n"
700" while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
701" (fBx-fAx) > fmax( fabs(fAx), fabs(fBx)) * fXEps )\n"
702" {\n"
703" if (bHasToInterpolate)\n"
704" {\n"
705" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
706" {\n"
707" fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)\n"
708" + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)\n"
709" + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n"
710" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
711" }\n"
712" else\n"
713" bHasToInterpolate = false;\n"
714" }\n"
715" if(!bHasToInterpolate)\n"
716" {\n"
717" fSx = 0.5 * (fAx + fBx);\n"
718" fPx = fAx; fPy = fAy;\n"
719" fQx = fBx; fQy = fBy;\n"
720" bHasToInterpolate = true;\n"
721" }\n"
722" fPx = fQx; fQx = fRx; fRx = fSx;\n"
723" fPy = fQy; fQy = fRy; fRy = fp - GetChiSqDistCDF(fSx, fdf);\n"
724" if (lcl_HasChangeOfSign( fAy, fRy))\n"
725" {\n"
726" fBx = fRx; fBy = fRy;\n"
727" }\n"
728" else\n"
729" {\n"
730" fAx = fRx; fAy = fRy;\n"
731" }\n"
732" bHasToInterpolate = bHasToInterpolate && (fabs(fRy) * 2.0"
733" <= fabs(fQy));\n"
734" ++nCount;\n"
735" }\n"
736" return fRx;\n"
737"}\n";
738
739const char gaussinvDecl[] = "double gaussinv(double x);\n";
740const char gaussinv[] =
741"double gaussinv(double x)\n"
742"{\n"
743" double q,t,z;\n"
744" q=x-0.5;\n"
745" if(fabs(q)<=.425)\n"
746" {\n"
747" t=0.180625-q*q;\n"
748" z=\n"
749" q*\n"
750" (\n"
751" (\n"
752" (\n"
753" (\n"
754" (\n"
755" (\n"
756" (\n"
757" t*2509.0809287301226727+33430.575583588128105\n"
758" )\n"
759" *t+67265.770927008700853\n"
760" )\n"
761" *t+45921.953931549871457\n"
762" )\n"
763" *t+13731.693765509461125\n"
764" )\n"
765" *t+1971.5909503065514427\n"
766" )\n"
767" *t+133.14166789178437745\n"
768" )\n"
769" *t+3.387132872796366608\n"
770" )\n"
771" /\n"
772" (\n"
773" (\n"
774" (\n"
775" (\n"
776" (\n"
777" (\n"
778" (\n"
779" t*5226.495278852854561+28729.085735721942674\n"
780" )\n"
781" *t+39307.89580009271061\n"
782" )\n"
783" *t+21213.794301586595867\n"
784" )\n"
785" *t+5394.1960214247511077\n"
786" )\n"
787" *t+687.1870074920579083\n"
788" )\n"
789" *t+42.313330701600911252\n"
790" )\n"
791" *t+1.0);\n"
792" }\n"
793" else\n"
794" {\n"
795" if(q>0) t=1-x;\n"
796" else t=x;\n"
797" t=sqrt(-log(t));\n"
798" if(t<=5.0)\n"
799" {\n"
800" t+=-1.6;\n"
801" z=\n"
802" (\n"
803" (\n"
804" (\n"
805" (\n"
806" (\n"
807" (\n"
808" (\n"
809" t*7.7454501427834140764e-4+0.0227238449892691845833\n"
810" )\n"
811" *t+0.24178072517745061177\n"
812" )\n"
813" *t+1.27045825245236838258\n"
814" )\n"
815" *t+3.64784832476320460504\n"
816" )\n"
817" *t+5.7694972214606914055\n"
818" )\n"
819" *t+4.6303378461565452959\n"
820" )\n"
821" *t+1.42343711074968357734\n"
822" )\n"
823" /\n"
824" (\n"
825" (\n"
826" (\n"
827" (\n"
828" (\n"
829" (\n"
830" (\n"
831" t*1.05075007164441684324e-9+5.475938084995344946e-4\n"
832" )\n"
833" *t+0.0151986665636164571966\n"
834" )\n"
835" *t+0.14810397642748007459\n"
836" )\n"
837" *t+0.68976733498510000455\n"
838" )\n"
839" *t+1.6763848301838038494\n"
840" )\n"
841" *t+2.05319162663775882187\n"
842" )\n"
843" *t+1.0);\n"
844" }\n"
845" else\n"
846" {\n"
847" t+=-5.0;\n"
848" z=\n"
849" (\n"
850" (\n"
851" (\n"
852" (\n"
853" (\n"
854" (\n"
855" (\n"
856" t*2.01033439929228813265e-7+2.71155556874348757815e-5\n"
857" )\n"
858" *t+0.0012426609473880784386\n"
859" )\n"
860" *t+0.026532189526576123093\n"
861" )\n"
862" *t+0.29656057182850489123\n"
863" )\n"
864" *t+1.7848265399172913358\n"
865" )\n"
866" *t+5.4637849111641143699\n"
867" )\n"
868" *t+6.6579046435011037772\n"
869" )\n"
870" /\n"
871" (\n"
872" (\n"
873" (\n"
874" (\n"
875" (\n"
876" (\n"
877" (\n"
878" t*2.04426310338993978564e-15+1.4215117583164458887e-7\n"
879" )\n"
880" *t+1.8463183175100546818e-5\n"
881" )\n"
882" *t+7.868691311456132591e-4\n"
883" )\n"
884" *t+0.0148753612908506148525\n"
885" )\n"
886" *t+0.13692988092273580531\n"
887" )\n"
888" *t+0.59983220655588793769\n"
889" )\n"
890" *t+1.0);\n"
891" }\n"
892" if(q<0.0) z=-z;\n"
893" }\n"
894" return z;\n"
895"}\n";
896
897const char lcl_GetLogGammaHelperDecl[] =
898"static double lcl_GetLogGammaHelper(double fZ);\n";
899const char lcl_GetLogGammaHelper[] =
900"static double lcl_GetLogGammaHelper(double fZ)\n"
901"{\n"
902" double fg = 6.024680040776729583740234375;\n"
903" double fZgHelp = fZ + fg - 0.5;\n"
904" return log( lcl_getLanczosSum(fZ)) + (fZ-0.5) * log(fZgHelp) - fZgHelp;\n"
905"}\n";
906const char lcl_GetGammaHelperDecl[] =
907"static double lcl_GetGammaHelper(double fZ);\n";
908const char lcl_GetGammaHelper[] =
909"static double lcl_GetGammaHelper(double fZ)\n"
910"{\n"
911" double fGamma = lcl_getLanczosSum(fZ);\n"
912" double fg = 6.024680040776729583740234375;\n"
913" double fZgHelp = fZ + fg - 0.5;\n"
914" double fHalfpower = pow( fZgHelp, fZ/2 - 0.25);\n"
915" fGamma *= fHalfpower;\n"
916" fGamma = fGamma/exp(fZgHelp);\n"
917" fGamma *= fHalfpower;\n"
918" fGamma = 120.4;\n"
919" if (fZ <= 20.0 && fZ == (int)fZ)\n"
920" {\n"
921" fGamma = (int)(fGamma+0.5);\n"
922" }\n"
923" return fGamma;\n"
924"}\n";
925const char lcl_getLanczosSumDecl[] =
926"static double lcl_getLanczosSum(double fZ);\n";
927const char lcl_getLanczosSum[] =
928"static double lcl_getLanczosSum(double fZ) \n"
929"{ \n"
930" double fNum[13] ={ \n"
931" 23531376880.41075968857200767445163675473, \n"
932" 42919803642.64909876895789904700198885093, \n"
933" 35711959237.35566804944018545154716670596, \n"
934" 17921034426.03720969991975575445893111267, \n"
935" 6039542586.35202800506429164430729792107, \n"
936" 1439720407.311721673663223072794912393972, \n"
937" 248874557.8620541565114603864132294232163, \n"
938" 31426415.58540019438061423162831820536287, \n"
939" 2876370.628935372441225409051620849613599, \n"
940" 186056.2653952234950402949897160456992822, \n"
941" 8071.672002365816210638002902272250613822, \n"
942" 210.8242777515793458725097339207133627117, \n"
943" 2.506628274631000270164908177133837338626 \n"
944" }; \n"
945" double fDenom[13] = { \n"
946" 0,\n"
947" 39916800,\n"
948" 120543840,\n"
949" 150917976,\n"
950" 105258076,\n"
951" 45995730,\n"
952" 13339535,\n"
953" 2637558,\n"
954" 357423,\n"
955" 32670,\n"
956" 1925,\n"
957" 66,\n"
958" 1\n"
959" };\n"
960" double fSumNum;\n"
961" double fSumDenom;\n"
962" int nI;\n"
963" if (fZ<=1.0)\n"
964" {\n"
965" fSumNum = fNum[12];\n"
966" fSumDenom = fDenom[12];\n"
967" nI = 11;\n"
968" fSumNum = fSumNum*fZ+fNum[nI];\n"
969" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
970" nI = 10;\n"
971" fSumNum = fSumNum*fZ+fNum[nI];\n"
972" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
973" nI = 9;\n"
974" fSumNum = fSumNum*fZ+fNum[nI];\n"
975" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
976" nI = 8;\n"
977" fSumNum = fSumNum*fZ+fNum[nI];\n"
978" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
979" nI = 7;\n"
980" fSumNum = fSumNum*fZ+fNum[nI];\n"
981" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
982" nI = 6;\n"
983" fSumNum = fSumNum*fZ+fNum[nI];\n"
984" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
985" nI = 5;\n"
986" fSumNum = fSumNum*fZ+fNum[nI];\n"
987" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
988" nI = 4;\n"
989" fSumNum = fSumNum*fZ+fNum[nI];\n"
990" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
991" nI = 3;\n"
992" fSumNum = fSumNum*fZ+fNum[nI];\n"
993" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
994" nI = 2;\n"
995" fSumNum = fSumNum*fZ+fNum[nI];\n"
996" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
997" nI = 1;\n"
998" fSumNum = fSumNum*fZ+fNum[nI];\n"
999" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
1000" nI = 0;\n"
1001" fSumNum = fSumNum*fZ+fNum[nI];\n"
1002" fSumDenom = fSumDenom*fZ+fDenom[nI];\n"
1003" }\n"
1004" if (fZ>1.0)\n"
1005" {\n"
1006" double fZInv = 1.0/fZ;\n"
1007" fSumNum = fNum[0];\n"
1008" fSumDenom = fDenom[0];\n"
1009" nI = 1;\n"
1010" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1011" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1012" nI = 2;\n"
1013" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1014" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1015" nI = 3;\n"
1016" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1017" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1018" nI = 4;\n"
1019" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1020" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1021" nI = 5;\n"
1022" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1023" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1024" nI = 6;\n"
1025" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1026" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1027" nI = 7;\n"
1028" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1029" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1030" nI = 8;\n"
1031" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1032" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1033" nI = 9;\n"
1034" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1035" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1036" nI = 10;\n"
1037" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1038" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1039" nI = 11;\n"
1040" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1041" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1042" nI = 12;\n"
1043" fSumNum = fSumNum*fZInv+fNum[nI];\n"
1044" fSumDenom = fSumDenom*fZInv+fDenom[nI];\n"
1045" }\n"
1046" return fSumNum/fSumDenom;\n"
1047"}\n";
1048
1049const char GetUpRegIGammaDecl[] =
1050" double GetUpRegIGamma( double fA, double fX ) ;\n";
1051const char GetUpRegIGamma[] =
1052"double GetUpRegIGamma( double fA, double fX )\n"
1053"{\n"
1054" double fLnFactor= fA*log(fX)-fX-lgamma(fA);\n"
1055" double fFactor = exp(fLnFactor); \n"
1056" if (fX>fA+1.0) \n"
1057" return fFactor * GetGammaContFraction(fA,fX);\n"
1058" else \n"
1059" return 1.0 -fFactor * GetGammaSeries(fA,fX);\n"
1060"}\n";
1061
1062const char lcl_HasChangeOfSignDecl[] =
1063"static inline bool lcl_HasChangeOfSign( double u, double w );\n";
1064const char lcl_HasChangeOfSign[] =
1065"static inline bool lcl_HasChangeOfSign( double u, double w )\n"
1066"{\n"
1067" return (u < 0.0 && w > 0.0) || (u > 0.0 && w < 0.0);\n"
1068"}\n";
1069
1070const char GetTDistDecl[] =" double GetTDist(double T, double fDF);\n";
1071const char GetTDist[] =
1072"double GetTDist(double T, double fDF)\n"
1073"{\n"
1074" return 0.5 * GetBetaDist(fDF/(fDF+T*T),fDF/2.0, 0.5);\n"
1075"}\n";
1076
1077const char GetBetaDecl[] =" double GetBeta(double fAlpha, double fBeta);\n";
1078const char GetBeta[] =
1079"double GetBeta(double fAlpha, double fBeta)\n"
1080"{\n"
1081" double fA;\n"
1082" double fB;\n"
1083" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
1084" double fAB = fA + fB;\n"
1085
1086" if (fAB < fMaxGammaArgument)\n"
1087" return tgamma(fA)/tgamma(fAB)*tgamma(fB);\n"
1088" double fgm = 5.524680040776729583740234375;\n"
1089" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n"
1090" /lcl_getLanczosSum(fAB);\n"
1091" fLanczos *= sqrt(((fAB + fgm)/(fA + fgm))/(fB + fgm));\n"
1092" return fLanczos * pow(exp(1.0),(-fA*log1p(fB/(fA + fgm)))"
1093" - fB*log1p(fA/(fB + fgm)) - fgm);\n"
1094"}\n";
1095
1096const char GetLogBetaDecl[] =
1097" double GetLogBeta(double fAlpha, double fBeta);\n";
1098const char GetLogBeta[] =
1099"double GetLogBeta(double fAlpha, double fBeta)\n"
1100"{\n"
1101" double fA;\n"
1102" double fB;\n"
1103" fAlpha>fBeta?(fA = fAlpha,fB = fBeta):(fA = fBeta,fB = fAlpha);\n"
1104" double fgm = 5.524680040776729583740234375;\n"
1105
1106" double fLanczos = lcl_getLanczosSum(fA)*lcl_getLanczosSum(fB)\n"
1107" /lcl_getLanczosSum(fA + fB);\n"
1108" double fResult= -fA *log1p(fB/(fA + fgm))"
1109"-fB *log1p(fA/(fB + fgm))-fgm;\n"
1110" fResult += log(fLanczos)+0.5*(log(fA + fB + fgm) - log(fA + fgm)\n"
1111" - log(fB + fgm));\n"
1112" return fResult;\n"
1113"}\n";
1114
1115const char GetBetaDistPDFDecl[] =
1116"double GetBetaDistPDF(double fX, double fA, double fB);\n";
1117const char GetBetaDistPDF[] =
1118"double GetBetaDistPDF(double fX, double fA, double fB)\n"
1119"{\n"
1120" if (fA == 1.0) \n"
1121" {\n"
1122" if (fB == 1.0)\n"
1123" return 1.0;\n"
1124" if (fB == 2.0)\n"
1125" return -2.0*fX + 2.0;\n"
1126" if (fX == 1.0 && fB < 1.0)\n"
1127" {\n"
1128" return HUGE_VAL;\n"
1129" }\n"
1130" if (fX <= 0.01)\n"
1131" return fB + fB * expm1((fB-1.0) * log1p(-fX));\n"
1132" else \n"
1133" return fB * pow(0.5-fX+0.5,fB-1.0);\n"
1134" }\n"
1135" if (fB == 1.0) \n"
1136" {\n"
1137" if (fA == 2.0)\n"
1138" return fA * fX;\n"
1139" if (fX == 0.0 && fA < 1.0)\n"
1140" {\n"
1141" return HUGE_VAL;\n"
1142" }\n"
1143" return fA * pow(fX,fA-1);\n"
1144" }\n"
1145" if (fX <= 0.0)\n"
1146" {\n"
1147" if (fA < 1.0 && fX == 0.0)\n"
1148" {\n"
1149" return HUGE_VAL;\n"
1150" }\n"
1151" else\n"
1152" return 0.0;\n"
1153" }\n"
1154" if (fX >= 1.0)\n"
1155" {\n"
1156" if (fB < 1.0 && fX == 1.0)\n"
1157" {\n"
1158" return HUGE_VAL;\n"
1159" }\n"
1160" else \n"
1161" return 0.0;\n"
1162" }\n"
1163" double fLogDblMax = log( 1.79769e+308 );\n"
1164" double fLogDblMin = log( 2.22507e-308 );\n"
1165" double fLogY = (fX < 0.1) ? log1p(-fX) : log(0.5-fX+0.5);\n"
1166" double fLogX = log(fX);\n"
1167" double fAm1LogX = (fA-1.0) * fLogX;\n"
1168" double fBm1LogY = (fB-1.0) * fLogY;\n"
1169" double fLogBeta = GetLogBeta(fA,fB);\n"
1170" if ( fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin\n"
1171" && fBm1LogY < fLogDblMax && fBm1LogY > fLogDblMin\n"
1172" && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin\n"
1173" && fAm1LogX + fBm1LogY < fLogDblMax && fAm1LogX + fBm1LogY > \n"
1174" fLogDblMin)\n"
1175" return pow(fX,fA-1.0)*pow(0.5-fX+0.5,fB-1.0)/GetBeta(fA,fB);\n"
1176" else \n"
1177" return exp( fAm1LogX + fBm1LogY - fLogBeta);\n"
1178"}\n";
1179
1180const char lcl_GetBetaHelperContFracDecl[] =
1181"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB);\n";
1182const char lcl_GetBetaHelperContFrac[] =
1183"double lcl_GetBetaHelperContFrac(double fX, double fA, double fB)\n"
1184"{ \n"
1185
1186" double a1, b1, a2, b2, fnorm, apl2m, d2m, d2m1, cfnew, cf;\n"
1187" a1 = 1.0; b1 = 1.0;\n"
1188" b2 = 1.0 - (fA+fB)/(fA+1.0)*fX;\n"
1189" b2==0.0?(a2 = 0.0,fnorm = 1.0,cf = 1.0):\n"
1190" (a2 = 1.0,fnorm = 1.0/b2,cf = a2*fnorm);\n"
1191" cfnew = 1.0;\n"
1192" double rm = 1.0;\n"
1193" double fMaxIter = 50000.0;\n"
1194" bool bfinished = false;\n"
1195" do\n"
1196" {\n"
1197" apl2m = fA + 2.0*rm;\n"
1198" d2m = (rm*(fB-rm))*fX/(apl2m*(apl2m-1.0));\n"
1199" d2m1 = -((fA+rm)*(fA+rm+fB))*fX/(apl2m*(apl2m+1.0));\n"
1200" a1 = (a2+d2m*a1)*fnorm;\n"
1201" b1 = (b2+d2m*b1)*fnorm;\n"
1202" a2 = a1 + d2m1*a2*fnorm;\n"
1203" b2 = b1 + d2m1*b2*fnorm;\n"
1204" if (b2 != 0.0) \n"
1205" {\n"
1206" fnorm = 1.0/b2;\n"
1207" cfnew = a2*fnorm;\n"
1208" bfinished = (fabs(cf-cfnew) < fabs(cf)*fMachEps);\n"
1209" }\n"
1210" cf = cfnew;\n"
1211" rm += 1.0;\n"
1212" }\n"
1213" while (rm < fMaxIter && !bfinished);\n"
1214" return cf;\n"
1215"}\n";
1216
1217const char lcl_IterateInverseDecl[] =
1218"double lcl_IterateInverse("
1219"double fAx, double fBx, bool* rConvError,double fp,double fDF );\n";
1220const char lcl_IterateInverse[] =
1221"double lcl_IterateInverse( "
1222"double fAx, double fBx, bool* rConvError,double fp,double fDF )\n"
1223"{\n"
1224" *rConvError = false;\n"
1225" double fYEps = 1.0E-307;\n"
1226" double fXEps =DBL_EPSILON;\n"
1227" if(fAx>fBx)\n"
1228" return DBL_MAX;\n"
1229" double fAy = GetValue(fAx,fp,fDF);\n"
1230" double fBy = GetValue(fBx,fp,fDF);\n"
1231" double fTemp;\n"
1232" unsigned short nCount;\n"
1233" double inter;\n"
1234" bool sign;\n"
1235" for (nCount =0;nCount<1000&&!lcl_HasChangeOfSign(fAy,fBy);nCount++)\n"
1236" {\n"
1237" inter = 2.0 * (fAx - fBx);\n"
1238" if (fabs(fAy) <= fabs(fBy)) \n"
1239" {\n"
1240" sign = true;\n"
1241" fTemp = fAx;\n"
1242" fAx += inter;\n"
1243" if (fAx < 0.0)\n"
1244" fAx = 0.0;\n"
1245" fBx = fTemp;\n"
1246" fBy = fAy;\n"
1247" fTemp = fAx;\n"
1248" }\n"
1249" else\n"
1250" {\n"
1251" sign = false;\n"
1252" fTemp = fBx;\n"
1253" fBx -= inter;\n"
1254" fAx = fTemp;\n"
1255" fAy = fBy;\n"
1256" fTemp = fBx;\n"
1257" }\n"
1258" fTemp = GetValue(fTemp,fp,fDF);\n"
1259" sign?(fAy = fTemp):(fBy = fTemp);\n"
1260" }\n"
1261" if (fAy == 0.0)\n"
1262" return fAx;\n"
1263" if (fBy == 0.0)\n"
1264" return fBx;\n"
1265" if (!lcl_HasChangeOfSign( fAy, fBy))\n"
1266" {\n"
1267" *rConvError = true;\n"
1268" return 0.0;\n"
1269" }\n"
1270" double fPx = fAx;\n"
1271" double fPy = fAy;\n"
1272" double fQx = fBx;\n"
1273" double fQy = fBy;\n"
1274" double fRx = fAx;\n"
1275" double fRy = fAy;\n"
1276" double fSx = 0.5 * (fAx + fBx); \n"
1277" bool bHasToInterpolate = true;\n"
1278" nCount = 0;\n"
1279" while ( nCount < 500 && fabs(fRy) > fYEps &&\n"
1280" (fBx-fAx) > max( fabs(fAx), fabs(fBx)) * fXEps )\n"
1281" {\n"
1282" if (bHasToInterpolate)\n"
1283" {\n"
1284" if (fPy!=fQy && fQy!=fRy && fRy!=fPy)\n"
1285" {\n"
1286" fSx = fPx * fRy * fQy / (fRy-fPy) / (fQy-fPy)\n"
1287" + fRx * fQy * fPy / (fQy-fRy) / (fPy-fRy)\n"
1288" + fQx * fPy * fRy / (fPy-fQy) / (fRy-fQy);\n"
1289" bHasToInterpolate = (fAx < fSx) && (fSx < fBx);\n"
1290" }\n"
1291" else\n"
1292" bHasToInterpolate = false;\n"
1293" }\n"
1294" if(!bHasToInterpolate)\n"
1295" {\n"
1296" fSx = 0.5 * (fAx + fBx);\n"
1297" \n"
1298" fPx = fAx; fPy = fAy;\n"
1299" fQx = fBx; fQy = fBy;\n"
1300" bHasToInterpolate = true;\n"
1301" }\n"
1302" fPx = fQx; fQx = fRx; fRx = fSx;\n"
1303" fPy = fQy; fQy = fRy; fRy = GetValue(fSx,fp,fDF);\n"
1304" lcl_HasChangeOfSign( fAy, fRy)?(fBx = fRx,fBy = fRy):\n"
1305" (fAx = fRx,fAy = fRy);\n"
1306" bHasToInterpolate =\n"
1307" bHasToInterpolate && (fabs(fRy) * 2.0 <= fabs(fQy));\n"
1308" ++nCount;\n"
1309" }\n"
1310" return fRx;\n"
1311"}\n";
1312const char phiDecl[] =
1313"double phi(double x);\n";
1314const char phi[] =
1315"double phi(double x)\n"
1316"{\n"
1317" return 0.39894228040143268 * exp(-(x * x) / 2.0);\n"
1318"}\n";
1319const char taylorDecl[] =
1320"double taylor(double* pPolynom, uint nMax, double x);\n";
1321const char taylor[] =
1322"double taylor(double* pPolynom, uint nMax, double x)\n"
1323"{\n"
1324" double nVal = pPolynom[nMax];\n"
1325" for (short i = nMax-1; i >= 0; i--)\n"
1326" {\n"
1327" nVal = pPolynom[i] + (nVal * x);\n"
1328" }\n"
1329" return nVal;\n"
1330"}";
1331const char gaussDecl[] = "double gauss(double x);\n";
1332const char gauss[] =
1333"double gauss(double x)\n"
1334"{\n"
1335" double xAbs = fabs(x);\n"
1336" uint xShort = (uint)(floor(xAbs));\n"
1337" double nVal = 0.0;\n"
1338" if (xShort == 0)\n"
1339" {\n"
1340" double t0[] =\n"
1341" { 0.39894228040143268, -0.06649038006690545, 0.00997355701003582,\n"
1342" -0.00118732821548045, 0.00011543468761616, -0.00000944465625950,\n"
1343" 0.00000066596935163, -0.00000004122667415, 0.00000000227352982,\n"
1344" 0.00000000011301172, 0.00000000000511243, -0.00000000000021218 };\n"
1345" nVal = taylor(t0, 11, (xAbs * xAbs)) * xAbs;\n"
1346" }\n"
1347" else if ((xShort >= 1) && (xShort <= 2))\n"
1348" {\n"
1349" double t2[] =\n"
1350" { 0.47724986805182079, 0.05399096651318805, -0.05399096651318805,\n"
1351" 0.02699548325659403, -0.00449924720943234, -0.00224962360471617,\n"
1352" 0.00134977416282970, -0.00011783742691370, -0.00011515930357476,\n"
1353" 0.00003704737285544, 0.00000282690796889, -0.00000354513195524,\n"
1354" 0.00000037669563126, 0.00000019202407921, -0.00000005226908590,\n"
1355" -0.00000000491799345, 0.00000000366377919, -0.00000000015981997,\n"
1356" -0.00000000017381238, 0.00000000002624031, 0.00000000000560919,\n"
1357" -0.00000000000172127, -0.00000000000008634, 0.00000000000007894 };\n"
1358" nVal = taylor(t2, 23, (xAbs - 2.0));\n"
1359" }\n"
1360" else if ((xShort >= 3) && (xShort <= 4))\n"
1361" {\n"
1362" double t4[] =\n"
1363" { 0.49996832875816688, 0.00013383022576489, -0.00026766045152977,\n"
1364" 0.00033457556441221, -0.00028996548915725, 0.00018178605666397,\n"
1365" -0.00008252863922168, 0.00002551802519049, -0.00000391665839292,\n"
1366" -0.00000074018205222, 0.00000064422023359, -0.00000017370155340,\n"
1367" 0.00000000909595465, 0.00000000944943118, -0.00000000329957075,\n"
1368" 0.00000000029492075, 0.00000000011874477, -0.00000000004420396,\n"
1369" 0.00000000000361422, 0.00000000000143638, -0.00000000000045848 };\n"
1370" nVal = taylor(t4, 20, (xAbs - 4.0));\n"
1371" }\n"
1372" else\n"
1373" {\n"
1374" double asympt[] = { -1.0, 1.0, -3.0, 15.0, -105.0 };\n"
1375" nVal = 0.5 + phi(xAbs) * taylor(asympt, 4, 1.0/(xAbs * xAbs))/xAbs;\n"
1376" }\n"
1377" if (x < 0.0)\n"
1378" return -nVal;\n"
1379" else\n"
1380" return nVal;\n"
1381"}\n";
1382
1383/* vim:set shiftwidth=4 softtabstop=4 expandtab: */
lcl_GetGammaHelperDecl
const char lcl_GetGammaHelperDecl[]
Definition: op_statistical_helpers.hxx:906
GetGammaContFraction
const char GetGammaContFraction[]
Definition: op_statistical_helpers.hxx:72
GetFInvValue
const char GetFInvValue[]
Definition: op_statistical_helpers.hxx:272
lcl_GetLogGammaHelperDecl
const char lcl_GetLogGammaHelperDecl[]
Definition: op_statistical_helpers.hxx:897
GetGammaInvValueDecl
const char GetGammaInvValueDecl[]
Definition: op_statistical_helpers.hxx:252
GetLogGamma
const char GetLogGamma[]
Definition: op_statistical_helpers.hxx:375
fMachEpsDecl
const char fMachEpsDecl[]
Definition: op_statistical_helpers.hxx:14
GetLogBetaDecl
const char GetLogBetaDecl[]
Definition: op_statistical_helpers.hxx:1096
lcl_IterateInverseBetaInv
const char lcl_IterateInverseBetaInv[]
Definition: op_statistical_helpers.hxx:450
GetBetaDecl
const char GetBetaDecl[]
Definition: op_statistical_helpers.hxx:1077
GetChiSqDistPDF
const char GetChiSqDistPDF[]
Definition: op_statistical_helpers.hxx:410
lcl_IterateInverseDecl
const char lcl_IterateInverseDecl[]
Definition: op_statistical_helpers.hxx:1217
taylor
const char taylor[]
Definition: op_statistical_helpers.hxx:1321
lcl_IterateInverseBetaInvDecl
const char lcl_IterateInverseBetaInvDecl[]
Definition: op_statistical_helpers.hxx:447
lcl_GetGammaHelper
const char lcl_GetGammaHelper[]
Definition: op_statistical_helpers.hxx:908
taylorDecl
const char taylorDecl[]
Definition: op_statistical_helpers.hxx:1319
GetLogBeta
const char GetLogBeta[]
Definition: op_statistical_helpers.hxx:1098
GetBetaDistPDFDecl
const char GetBetaDistPDFDecl[]
Definition: op_statistical_helpers.hxx:1115
GetUpRegIGamma
const char GetUpRegIGamma[]
Definition: op_statistical_helpers.hxx:1051
fMaxGammaArgumentDecl
const char fMaxGammaArgumentDecl[]
Definition: op_statistical_helpers.hxx:17
gaussinvDecl
const char gaussinvDecl[]
Definition: op_statistical_helpers.hxx:739
fsum_countDecl
const char fsum_countDecl[]
Definition: op_statistical_helpers.hxx:42
GetFInvValueDecl
const char GetFInvValueDecl[]
Definition: op_statistical_helpers.hxx:270
lcl_GetBetaHelperContFracDecl
const char lcl_GetBetaHelperContFracDecl[]
Definition: op_statistical_helpers.hxx:1180
GetBetaDistDecl
const char GetBetaDistDecl[]
Definition: op_statistical_helpers.hxx:196
GetGammaInvValue
const char GetGammaInvValue[]
Definition: op_statistical_helpers.hxx:254
gaussinv
const char gaussinv[]
Definition: op_statistical_helpers.hxx:740
GetLogGammaDecl
const char GetLogGammaDecl[]
Definition: op_statistical_helpers.hxx:374
GetBinomDistPMF
const char GetBinomDistPMF[]
Definition: op_statistical_helpers.hxx:325
GetLowRegIGamma
const char GetLowRegIGamma[]
Definition: op_statistical_helpers.hxx:122
gaussDecl
const char gaussDecl[]
Definition: op_statistical_helpers.hxx:1331
lcl_IterateInverseChiInv
const char lcl_IterateInverseChiInv[]
Definition: op_statistical_helpers.hxx:547
MinDecl
const char MinDecl[]
Definition: op_statistical_helpers.hxx:12
GetTDist
const char GetTDist[]
Definition: op_statistical_helpers.hxx:1071
GetChiDistDecl
const char GetChiDistDecl[]
Definition: op_statistical_helpers.hxx:387
lcl_IterateInverseChiSQInvDecl
const char lcl_IterateInverseChiSQInvDecl[]
Definition: op_statistical_helpers.hxx:641
fmax_count
const char fmax_count[]
Definition: op_statistical_helpers.hxx:35
GetValueDecl
const char GetValueDecl[]
Definition: op_statistical_helpers.hxx:19
GetBeta
const char GetBeta[]
Definition: op_statistical_helpers.hxx:1078
GetChiSqDistPDFDecl
const char GetChiSqDistPDFDecl[]
Definition: op_statistical_helpers.hxx:408
lcl_IterateInverse
const char lcl_IterateInverse[]
Definition: op_statistical_helpers.hxx:1220
lcl_GetBinomDistRange
const char lcl_GetBinomDistRange[]
Definition: op_statistical_helpers.hxx:355
phi
const char phi[]
Definition: op_statistical_helpers.hxx:1314
lcl_IterateInverseChiInvDecl
const char lcl_IterateInverseChiInvDecl[]
Definition: op_statistical_helpers.hxx:544
gauss
const char gauss[]
Definition: op_statistical_helpers.hxx:1332
GetLowRegIGammaDecl
const char GetLowRegIGammaDecl[]
Definition: op_statistical_helpers.hxx:120
lcl_getLanczosSumDecl
const char lcl_getLanczosSumDecl[]
Definition: op_statistical_helpers.hxx:925
fsum_count
const char fsum_count[]
Definition: op_statistical_helpers.hxx:43
GetGammaSeriesDecl
const char GetGammaSeriesDecl[]
Definition: op_statistical_helpers.hxx:49
fmin_countDecl
const char fmin_countDecl[]
Definition: op_statistical_helpers.hxx:26
GetGammaDistPDFDecl
const char GetGammaDistPDFDecl[]
Definition: op_statistical_helpers.hxx:142
lcl_HasChangeOfSignDecl
const char lcl_HasChangeOfSignDecl[]
Definition: op_statistical_helpers.hxx:1062
lcl_GetLogGammaHelper
const char lcl_GetLogGammaHelper[]
Definition: op_statistical_helpers.hxx:899
GetFDistDecl
const char GetFDistDecl[]
Definition: op_statistical_helpers.hxx:242
lcl_GetBinomDistRangeDecl
const char lcl_GetBinomDistRangeDecl[]
Definition: op_statistical_helpers.hxx:352
GetValue
const char GetValue[]
Definition: op_statistical_helpers.hxx:21
fHalfMachEpsDecl
const char fHalfMachEpsDecl[]
Definition: op_statistical_helpers.hxx:16
lcl_HasChangeOfSign
const char lcl_HasChangeOfSign[]
Definition: op_statistical_helpers.hxx:1064
GetChiSqDistCDFDecl
const char GetChiSqDistCDFDecl[]
Definition: op_statistical_helpers.hxx:397
fLogDblMaxDecl
const char fLogDblMaxDecl[]
Definition: op_statistical_helpers.hxx:15
GetTDistDecl
const char GetTDistDecl[]
Definition: op_statistical_helpers.hxx:1070
GetGammaContFractionDecl
const char GetGammaContFractionDecl[]
Definition: op_statistical_helpers.hxx:70
lcl_IterateInverseChiSQInv
const char lcl_IterateInverseChiSQInv[]
Definition: op_statistical_helpers.hxx:644
GetBetaDist
const char GetBetaDist[]
Definition: op_statistical_helpers.hxx:198
fmin_count
const char fmin_count[]
Definition: op_statistical_helpers.hxx:27
lcl_getLanczosSum
const char lcl_getLanczosSum[]
Definition: op_statistical_helpers.hxx:927
GetGammaDistPDF
const char GetGammaDistPDF[]
Definition: op_statistical_helpers.hxx:144
GetChiSqDistCDF
const char GetChiSqDistCDF[]
Definition: op_statistical_helpers.hxx:399
GetGammaSeries
const char GetGammaSeries[]
Definition: op_statistical_helpers.hxx:51
fBigInvDecl
const char fBigInvDecl[]
Definition: op_statistical_helpers.hxx:13
GetGammaDistDecl
const char GetGammaDistDecl[]
Definition: op_statistical_helpers.hxx:132
GetChiDist
const char GetChiDist[]
Definition: op_statistical_helpers.hxx:388
GetUpRegIGammaDecl
const char GetUpRegIGammaDecl[]
Definition: op_statistical_helpers.hxx:1049
GetBinomDistPMFDecl
const char GetBinomDistPMFDecl[]
Definition: op_statistical_helpers.hxx:323
fmax_countDecl
const char fmax_countDecl[]
Definition: op_statistical_helpers.hxx:34
phiDecl
const char phiDecl[]
Definition: op_statistical_helpers.hxx:1312
GetGammaDist
const char GetGammaDist[]
Definition: op_statistical_helpers.hxx:134
GetBetaDistPDF
const char GetBetaDistPDF[]
Definition: op_statistical_helpers.hxx:1117
GetFDist
const char GetFDist[]
Definition: op_statistical_helpers.hxx:244
lcl_GetBetaHelperContFrac
const char lcl_GetBetaHelperContFrac[]
Definition: op_statistical_helpers.hxx:1182
Generated on Sun Jul 30 2023 04:28:14 for LibreOffice Module sc (master) by doxygen 1.9.3