From 46336bc788ce344533524a29bc9858d183f91aeb Mon Sep 17 00:00:00 2001 From: xuri Date: Mon, 28 Mar 2022 08:13:47 +0800 Subject: [PATCH] ref #65, new formula functions: CHISQ.DIST.RT CHISQ.DIST and GAMMALN.PRECISE --- calc.go | 216 +++++++++++++++++++++++++++++++++++++++++++++++++++ calc_test.go | 31 ++++++++ 2 files changed, 247 insertions(+) diff --git a/calc.go b/calc.go index a87fa2f5..84f85682 100644 --- a/calc.go +++ b/calc.go @@ -357,6 +357,8 @@ type formulaFuncs struct { // CHIDIST // CHIINV // CHITEST +// CHISQ.DIST +// CHISQ.DIST.RT // CHISQ.TEST // CHOOSE // CLEAN @@ -449,6 +451,7 @@ type formulaFuncs struct { // GAMMA.INV // GAMMAINV // GAMMALN +// GAMMALN.PRECISE // GAUSS // GCD // GEOMEAN @@ -6416,6 +6419,200 @@ func (fn *formulaFuncs) CHITEST(argsList *list.List) formulaArg { return fn.CHIDIST(args) } +// getGammaSeries calculates a power-series of the gamma function. +func getGammaSeries(fA, fX float64) float64 { + var ( + fHalfMachEps = 2.22045e-016 / 2 + fDenomfactor = fA + fSummand = 1 / fA + fSum = fSummand + nCount = 1 + ) + for fSummand/fSum > fHalfMachEps && nCount <= 10000 { + fDenomfactor = fDenomfactor + 1 + fSummand = fSummand * fX / fDenomfactor + fSum = fSum + fSummand + nCount = nCount + 1 + } + return fSum +} + +// getGammaContFraction returns continued fraction with odd items of the gamma +// function. +func getGammaContFraction(fA, fX float64) float64 { + var ( + fBigInv = 2.22045e-016 + fHalfMachEps = fBigInv / 2 + fBig = 1 / fBigInv + fCount = 0.0 + fY = 1 - fA + fDenom = fX + 2 - fA + fPkm1 = fX + 1 + fPkm2 = 1.0 + fQkm1 = fDenom * fX + fQkm2 = fX + fApprox = fPkm1 / fQkm1 + bFinished = false + ) + for !bFinished && fCount < 10000 { + fCount = fCount + 1 + fY = fY + 1 + fDenom = fDenom + 2 + var ( + fNum = fY * fCount + f1 = fPkm1 * fDenom + f2 = fPkm2 * fNum + fPk = math.Nextafter(f1, f1) - math.Nextafter(f2, f2) + f3 = fQkm1 * fDenom + f4 = fQkm2 * fNum + fQk = math.Nextafter(f3, f3) - math.Nextafter(f4, f4) + ) + if fQk != 0 { + var fR = fPk / fQk + bFinished = math.Abs((fApprox-fR)/fR) <= fHalfMachEps + fApprox = fR + } + fPkm2, fPkm1, fQkm2, fQkm1 = fPkm1, fPk, fQkm1, fQk + if math.Abs(fPk) > fBig { + // reduce a fraction does not change the value + fPkm2 = fPkm2 * fBigInv + fPkm1 = fPkm1 * fBigInv + fQkm2 = fQkm2 * fBigInv + fQkm1 = fQkm1 * fBigInv + } + } + return fApprox +} + +// getLogGammaHelper is a part of implementation of the function getLogGamma. +func getLogGammaHelper(fZ float64) float64 { + var _fg = 6.024680040776729583740234375 + var zgHelp = fZ + _fg - 0.5 + return math.Log(getLanczosSum(fZ)) + (fZ-0.5)*math.Log(zgHelp) - zgHelp +} + +// getGammaHelper is a part of implementation of the function getLogGamma. +func getGammaHelper(fZ float64) float64 { + var ( + gamma = getLanczosSum(fZ) + fg = 6.024680040776729583740234375 + zgHelp = fZ + fg - 0.5 + // avoid intermediate overflow + halfpower = math.Pow(zgHelp, fZ/2-0.25) + ) + gamma *= halfpower + gamma /= math.Exp(zgHelp) + gamma *= halfpower + if fZ <= 20 && fZ == math.Floor(fZ) { + gamma = math.Round(gamma) + } + return gamma +} + +// getLogGamma calculates the natural logarithm of the gamma function. +func getLogGamma(fZ float64) float64 { + var fMaxGammaArgument = 171.624376956302 + if fZ >= fMaxGammaArgument { + return getLogGammaHelper(fZ) + } + if fZ >= 1.0 { + return math.Log(getGammaHelper(fZ)) + } + if fZ >= 0.5 { + return math.Log(getGammaHelper(fZ+1) / fZ) + } + return getLogGammaHelper(fZ+2) - math.Log(fZ+1) - math.Log(fZ) +} + +// getLowRegIGamma returns lower regularized incomplete gamma function. +func getLowRegIGamma(fA, fX float64) float64 { + fLnFactor := fA*math.Log(fX) - fX - getLogGamma(fA) + fFactor := math.Exp(fLnFactor) + if fX > fA+1 { + return 1 - fFactor*getGammaContFraction(fA, fX) + } + return fFactor * getGammaSeries(fA, fX) +} + +// getChiSqDistCDF returns left tail for the Chi-Square distribution. +func getChiSqDistCDF(fX, fDF float64) float64 { + if fX <= 0 { + return 0 + } + return getLowRegIGamma(fDF/2, fX/2) +} + +// getChiSqDistPDF calculates the probability density function for the +// Chi-Square distribution. +func getChiSqDistPDF(fX, fDF float64) float64 { + if fDF*fX > 1391000 { + return math.Exp((0.5*fDF-1)*math.Log(fX*0.5) - 0.5*fX - math.Log(2) - getLogGamma(0.5*fDF)) + } + var fCount, fValue float64 + if math.Mod(fDF, 2) < 0.5 { + fValue = 0.5 + fCount = 2 + } else { + fValue = 1 / math.Sqrt(fX*2*math.Pi) + fCount = 1 + } + for fCount < fDF { + fValue *= fX / fCount + fCount += 2 + } + if fX >= 1425 { + fValue = math.Exp(math.Log(fValue) - fX/2) + } else { + fValue *= math.Exp(-fX / 2) + } + return fValue +} + +// CHISQdotDIST function calculates the Probability Density Function or the +// Cumulative Distribution Function for the Chi-Square Distribution. The +// syntax of the function is: +// +// CHISQ.DIST(x,degrees_freedom,cumulative) +// +func (fn *formulaFuncs) CHISQdotDIST(argsList *list.List) formulaArg { + if argsList.Len() != 3 { + return newErrorFormulaArg(formulaErrorVALUE, "CHISQ.DIST requires 3 arguments") + } + var x, degrees, cumulative formulaArg + if x = argsList.Front().Value.(formulaArg).ToNumber(); x.Type != ArgNumber { + return x + } + if degrees = argsList.Front().Next().Value.(formulaArg).ToNumber(); degrees.Type != ArgNumber { + return degrees + } + if cumulative = argsList.Back().Value.(formulaArg).ToBool(); cumulative.Type == ArgError { + return cumulative + } + if x.Number < 0 { + return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) + } + maxDeg := math.Pow10(10) + if degrees.Number < 1 || degrees.Number >= maxDeg { + return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) + } + if cumulative.Number == 1 { + return newNumberFormulaArg(getChiSqDistCDF(x.Number, degrees.Number)) + } + return newNumberFormulaArg(getChiSqDistPDF(x.Number, degrees.Number)) +} + +// CHISQdotDISTdotRT function calculates the right-tailed probability of the +// Chi-Square Distribution. The syntax of the function is: +// +// CHISQ.DIST.RT(x,degrees_freedom) +// +func (fn *formulaFuncs) CHISQdotDISTdotRT(argsList *list.List) formulaArg { + if argsList.Len() != 2 { + return newErrorFormulaArg(formulaErrorVALUE, "CHISQ.DIST.RT requires 2 numeric arguments") + } + return fn.CHIDIST(argsList) +} + // CHISQdotTEST function performs the chi-square test on two supplied data sets // (of observed and expected frequencies), and returns the probability that // the differences between the sets are simply due to sampling error. The @@ -7033,6 +7230,25 @@ func (fn *formulaFuncs) GAMMALN(argsList *list.List) formulaArg { return newErrorFormulaArg(formulaErrorVALUE, "GAMMALN requires 1 numeric argument") } +// GAMMALNdotPRECISE function returns the natural logarithm of the Gamma +// Function, Γ(n). The syntax of the function is: +// +// GAMMALN.PRECISE(x) +// +func (fn *formulaFuncs) GAMMALNdotPRECISE(argsList *list.List) formulaArg { + if argsList.Len() != 1 { + return newErrorFormulaArg(formulaErrorVALUE, "GAMMALN.PRECISE requires 1 numeric argument") + } + x := argsList.Front().Value.(formulaArg).ToNumber() + if x.Type != ArgNumber { + return x + } + if x.Number <= 0 { + return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM) + } + return newNumberFormulaArg(getLogGamma(x.Number)) +} + // GAUSS function returns the probability that a member of a standard normal // population will fall between the mean and a specified number of standard // deviations from the mean. The syntax of the function is: diff --git a/calc_test.go b/calc_test.go index 308db55d..306e24d9 100644 --- a/calc_test.go +++ b/calc_test.go @@ -851,6 +851,20 @@ func TestCalcCellValue(t *testing.T) { "=CHIINV(0.75,1)": "0.101531044267622", "=CHIINV(0.1,2)": "4.60517018598809", "=CHIINV(0.8,2)": "0.446287102628419", + // CHISQ.DIST + "=CHISQ.DIST(0,2,TRUE)": "0", + "=CHISQ.DIST(4,1,TRUE)": "0.954499736103642", + "=CHISQ.DIST(1180,1180,FALSE)": "0.00821093706387967", + "=CHISQ.DIST(2,1,FALSE)": "0.103776874355149", + "=CHISQ.DIST(3,2,FALSE)": "0.111565080074215", + "=CHISQ.DIST(2,3,FALSE)": "0.207553748710297", + "=CHISQ.DIST(1425,1,FALSE)": "3.88315098887099E-312", + "=CHISQ.DIST(3,2,TRUE)": "0.77686983985157", + // CHISQ.DIST.RT + "=CHISQ.DIST.RT(0.5,3)": "0.918891411654676", + "=CHISQ.DIST.RT(8,3)": "0.0460117056892314", + "=CHISQ.DIST.RT(40,4)": "4.32842260712097E-08", + "=CHISQ.DIST.RT(42,4)": "1.66816329414062E-08", // CONFIDENCE "=CONFIDENCE(0.05,0.07,100)": "0.0137197479028414", // CONFIDENCE.NORM @@ -918,6 +932,9 @@ func TestCalcCellValue(t *testing.T) { // GAMMALN "=GAMMALN(4.5)": "2.45373657084244", "=GAMMALN(INT(1))": "0", + // GAMMALN.PRECISE + "=GAMMALN.PRECISE(0.4)": "0.796677817701784", + "=GAMMALN.PRECISE(4.5)": "2.45373657084244", // GAUSS "=GAUSS(-5)": "-0.499999713348428", "=GAUSS(0)": "0", @@ -2523,6 +2540,17 @@ func TestCalcCellValue(t *testing.T) { "=CHIINV(0,1)": "#NUM!", "=CHIINV(2,1)": "#NUM!", "=CHIINV(0.5,0.5)": "#NUM!", + // CHISQ.DIST + "=CHISQ.DIST()": "CHISQ.DIST requires 3 arguments", + "=CHISQ.DIST(\"\",2,TRUE)": "strconv.ParseFloat: parsing \"\": invalid syntax", + "=CHISQ.DIST(3,\"\",TRUE)": "strconv.ParseFloat: parsing \"\": invalid syntax", + "=CHISQ.DIST(3,2,\"\")": "strconv.ParseBool: parsing \"\": invalid syntax", + "=CHISQ.DIST(-1,2,TRUE)": "#NUM!", + "=CHISQ.DIST(3,0,TRUE)": "#NUM!", + // CHISQ.DIST.RT + "=CHISQ.DIST.RT()": "CHISQ.DIST.RT requires 2 numeric arguments", + "=CHISQ.DIST.RT(\"\",3)": "strconv.ParseFloat: parsing \"\": invalid syntax", + "=CHISQ.DIST.RT(0.5,\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax", // CONFIDENCE "=CONFIDENCE()": "CONFIDENCE requires 3 numeric arguments", "=CONFIDENCE(\"\",0.07,100)": "strconv.ParseFloat: parsing \"\": invalid syntax", @@ -2621,6 +2649,9 @@ func TestCalcCellValue(t *testing.T) { "=GAMMALN(F1)": "GAMMALN requires 1 numeric argument", "=GAMMALN(0)": "#N/A", "=GAMMALN(INT(0))": "#N/A", + // GAMMALN.PRECISE + "=GAMMALN.PRECISE()": "GAMMALN.PRECISE requires 1 numeric argument", + "=GAMMALN.PRECISE(\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax", // GAUSS "=GAUSS()": "GAUSS requires 1 numeric argument", "=GAUSS(\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax",