ref #65, new formula functions: CHIINV and BETADIST
This commit is contained in:
parent
c3424a9a0f
commit
1da129a3df
283
calc.go
283
calc.go
|
|
@ -333,6 +333,7 @@ type formulaFuncs struct {
|
|||
// BESSELJ
|
||||
// BESSELK
|
||||
// BESSELY
|
||||
// BETADIST
|
||||
// BETAINV
|
||||
// BETA.INV
|
||||
// BIN2DEC
|
||||
|
|
@ -348,6 +349,7 @@ type formulaFuncs struct {
|
|||
// CEILING.PRECISE
|
||||
// CHAR
|
||||
// CHIDIST
|
||||
// CHIINV
|
||||
// CHOOSE
|
||||
// CLEAN
|
||||
// CODE
|
||||
|
|
@ -5200,6 +5202,257 @@ func (fn *formulaFuncs) AVERAGEIF(argsList *list.List) formulaArg {
|
|||
return newNumberFormulaArg(sum / count)
|
||||
}
|
||||
|
||||
// getBetaHelperContFrac continued fractions for the beta function.
|
||||
func getBetaHelperContFrac(fX, fA, fB float64) float64 {
|
||||
var a1, b1, a2, b2, fnorm, cfnew, cf, rm float64
|
||||
a1, b1, b2 = 1, 1, 1-(fA+fB)/(fA+1)*fX
|
||||
if b2 == 0 {
|
||||
a2, fnorm, cf = 0, 1, 1
|
||||
} else {
|
||||
a2, fnorm = 1, 1/b2
|
||||
cf = a2 * fnorm
|
||||
}
|
||||
cfnew, rm = 1, 1
|
||||
fMaxIter, fMachEps := 50000.0, 2.22045e-016
|
||||
bfinished := false
|
||||
for rm < fMaxIter && !bfinished {
|
||||
apl2m := fA + 2*rm
|
||||
d2m := rm * (fB - rm) * fX / ((apl2m - 1) * apl2m)
|
||||
d2m1 := -(fA + rm) * (fA + fB + rm) * fX / (apl2m * (apl2m + 1))
|
||||
a1 = (a2 + d2m*a1) * fnorm
|
||||
b1 = (b2 + d2m*b1) * fnorm
|
||||
a2 = a1 + d2m1*a2*fnorm
|
||||
b2 = b1 + d2m1*b2*fnorm
|
||||
if b2 != 0 {
|
||||
fnorm = 1 / b2
|
||||
cfnew = a2 * fnorm
|
||||
bfinished = (math.Abs(cf-cfnew) < math.Abs(cf)*fMachEps)
|
||||
}
|
||||
cf = cfnew
|
||||
rm += 1
|
||||
}
|
||||
return cf
|
||||
}
|
||||
|
||||
// getLanczosSum uses a variant of the Lanczos sum with a rational function.
|
||||
func getLanczosSum(fZ float64) float64 {
|
||||
num := []float64{
|
||||
23531376880.41075968857200767445163675473,
|
||||
42919803642.64909876895789904700198885093,
|
||||
35711959237.35566804944018545154716670596,
|
||||
17921034426.03720969991975575445893111267,
|
||||
6039542586.35202800506429164430729792107,
|
||||
1439720407.311721673663223072794912393972,
|
||||
248874557.8620541565114603864132294232163,
|
||||
31426415.58540019438061423162831820536287,
|
||||
2876370.628935372441225409051620849613599,
|
||||
186056.2653952234950402949897160456992822,
|
||||
8071.672002365816210638002902272250613822,
|
||||
210.8242777515793458725097339207133627117,
|
||||
2.506628274631000270164908177133837338626,
|
||||
}
|
||||
denom := []float64{
|
||||
0,
|
||||
39916800,
|
||||
120543840,
|
||||
150917976,
|
||||
105258076,
|
||||
45995730,
|
||||
13339535,
|
||||
2637558,
|
||||
357423,
|
||||
32670,
|
||||
1925,
|
||||
66,
|
||||
1,
|
||||
}
|
||||
var sumNum, sumDenom, zInv float64
|
||||
if fZ <= 1 {
|
||||
sumNum = num[12]
|
||||
sumDenom = denom[12]
|
||||
for i := 11; i >= 0; i-- {
|
||||
sumNum *= fZ
|
||||
sumNum += num[i]
|
||||
sumDenom *= fZ
|
||||
sumDenom += denom[i]
|
||||
}
|
||||
} else {
|
||||
zInv = 1 / fZ
|
||||
sumNum = num[0]
|
||||
sumDenom = denom[0]
|
||||
for i := 1; i <= 12; i++ {
|
||||
sumNum *= zInv
|
||||
sumNum += num[i]
|
||||
sumDenom *= zInv
|
||||
sumDenom += denom[i]
|
||||
}
|
||||
}
|
||||
return sumNum / sumDenom
|
||||
}
|
||||
|
||||
// getBeta return beta distribution.
|
||||
func getBeta(fAlpha, fBeta float64) float64 {
|
||||
var fA, fB float64
|
||||
if fAlpha > fBeta {
|
||||
fA = fAlpha
|
||||
fB = fBeta
|
||||
} else {
|
||||
fA = fBeta
|
||||
fB = fAlpha
|
||||
}
|
||||
const maxGammaArgument = 171.624376956302
|
||||
if fA+fB < maxGammaArgument {
|
||||
return math.Gamma(fA) / math.Gamma(fA+fB) * math.Gamma(fB)
|
||||
}
|
||||
fg := 6.024680040776729583740234375
|
||||
fgm := fg - 0.5
|
||||
fLanczos := getLanczosSum(fA)
|
||||
fLanczos /= getLanczosSum(fA + fB)
|
||||
fLanczos *= getLanczosSum(fB)
|
||||
fABgm := fA + fB + fgm
|
||||
fLanczos *= math.Sqrt((fABgm / (fA + fgm)) / (fB + fgm))
|
||||
fTempA := fB / (fA + fgm)
|
||||
fTempB := fA / (fB + fgm)
|
||||
fResult := math.Exp(-fA*math.Log1p(fTempA) - fB*math.Log1p(fTempB) - fgm)
|
||||
fResult *= fLanczos
|
||||
return fResult
|
||||
}
|
||||
|
||||
// getBetaDistPDF is an implementation for the Beta probability density
|
||||
// function.
|
||||
func getBetaDistPDF(fX, fA, fB float64) float64 {
|
||||
if fX <= 0 || fX >= 1 {
|
||||
return 0
|
||||
}
|
||||
fLogDblMax, fLogDblMin := math.Log(1.79769e+308), math.Log(2.22507e-308)
|
||||
fLogY := math.Log(0.5 - fX + 0.5)
|
||||
if fX < 0.1 {
|
||||
fLogY = math.Log1p(-fX)
|
||||
}
|
||||
fLogX := math.Log(fX)
|
||||
fAm1LogX := (fA - 1) * fLogX
|
||||
fBm1LogY := (fB - 1) * fLogY
|
||||
fLogBeta := getLogBeta(fA, fB)
|
||||
if fAm1LogX < fLogDblMax && fAm1LogX > fLogDblMin && fBm1LogY < fLogDblMax &&
|
||||
fBm1LogY > fLogDblMin && fLogBeta < fLogDblMax && fLogBeta > fLogDblMin &&
|
||||
fAm1LogX+fBm1LogY < fLogDblMax && fAm1LogX+fBm1LogY > fLogDblMin {
|
||||
return math.Pow(fX, fA-1) * math.Pow(0.5-fX+0.5, fB-1) / getBeta(fA, fB)
|
||||
}
|
||||
return math.Exp(fAm1LogX + fBm1LogY - fLogBeta)
|
||||
}
|
||||
|
||||
// getLogBeta return beta with logarithm.
|
||||
func getLogBeta(fAlpha, fBeta float64) float64 {
|
||||
var fA, fB float64
|
||||
if fAlpha > fBeta {
|
||||
fA, fB = fAlpha, fBeta
|
||||
} else {
|
||||
fA, fB = fBeta, fAlpha
|
||||
}
|
||||
fg := 6.024680040776729583740234375
|
||||
fgm := fg - 0.5
|
||||
fLanczos := getLanczosSum(fA)
|
||||
fLanczos /= getLanczosSum(fA + fB)
|
||||
fLanczos *= getLanczosSum(fB)
|
||||
fLogLanczos := math.Log(fLanczos)
|
||||
fABgm := fA + fB + fgm
|
||||
fLogLanczos += 0.5 * (math.Log(fABgm) - math.Log(fA+fgm) - math.Log(fB+fgm))
|
||||
fTempA := fB / (fA + fgm)
|
||||
fTempB := fA / (fB + fgm)
|
||||
fResult := -fA*math.Log1p(fTempA) - fB*math.Log1p(fTempB) - fgm
|
||||
fResult += fLogLanczos
|
||||
return fResult
|
||||
}
|
||||
|
||||
// getBetaDist is an implementation for the beta distribution function.
|
||||
func getBetaDist(fXin, fAlpha, fBeta float64) float64 {
|
||||
if fXin <= 0 {
|
||||
return 0
|
||||
}
|
||||
if fXin >= 1 {
|
||||
return 1
|
||||
}
|
||||
if fBeta == 1 {
|
||||
return math.Pow(fXin, fAlpha)
|
||||
}
|
||||
if fAlpha == 1 {
|
||||
return -math.Expm1(fBeta * math.Log1p(-fXin))
|
||||
}
|
||||
var fResult float64
|
||||
fY, flnY := (0.5-fXin)+0.5, math.Log1p(-fXin)
|
||||
fX, flnX := fXin, math.Log(fXin)
|
||||
fA, fB := fAlpha, fBeta
|
||||
bReflect := fXin > fAlpha/(fAlpha+fBeta)
|
||||
if bReflect {
|
||||
fA = fBeta
|
||||
fB = fAlpha
|
||||
fX = fY
|
||||
fY = fXin
|
||||
flnX = flnY
|
||||
flnY = math.Log(fXin)
|
||||
}
|
||||
fResult = getBetaHelperContFrac(fX, fA, fB) / fA
|
||||
fP, fQ := fA/(fA+fB), fB/(fA+fB)
|
||||
var fTemp float64
|
||||
if fA > 1 && fB > 1 && fP < 0.97 && fQ < 0.97 {
|
||||
fTemp = getBetaDistPDF(fX, fA, fB) * fX * fY
|
||||
} else {
|
||||
fTemp = math.Exp(fA*flnX + fB*flnY - getLogBeta(fA, fB))
|
||||
}
|
||||
fResult *= fTemp
|
||||
if bReflect {
|
||||
fResult = 0.5 - fResult + 0.5
|
||||
}
|
||||
return fResult
|
||||
}
|
||||
|
||||
// BETADIST function calculates the cumulative beta probability density
|
||||
// function for a supplied set of parameters. The syntax of the function is:
|
||||
//
|
||||
// BETADIST(x,alpha,beta,[A],[B])
|
||||
//
|
||||
func (fn *formulaFuncs) BETADIST(argsList *list.List) formulaArg {
|
||||
if argsList.Len() < 3 {
|
||||
return newErrorFormulaArg(formulaErrorVALUE, "BETADIST requires at least 3 arguments")
|
||||
}
|
||||
if argsList.Len() > 5 {
|
||||
return newErrorFormulaArg(formulaErrorVALUE, "BETADIST requires at most 5 arguments")
|
||||
}
|
||||
x := argsList.Front().Value.(formulaArg).ToNumber()
|
||||
if x.Type != ArgNumber {
|
||||
return x
|
||||
}
|
||||
alpha := argsList.Front().Next().Value.(formulaArg).ToNumber()
|
||||
if alpha.Type != ArgNumber {
|
||||
return alpha
|
||||
}
|
||||
beta := argsList.Front().Next().Next().Value.(formulaArg).ToNumber()
|
||||
if beta.Type != ArgNumber {
|
||||
return beta
|
||||
}
|
||||
if alpha.Number <= 0 || beta.Number <= 0 {
|
||||
return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM)
|
||||
}
|
||||
a, b := newNumberFormulaArg(0), newNumberFormulaArg(1)
|
||||
if argsList.Len() > 3 {
|
||||
if a = argsList.Front().Next().Next().Next().Value.(formulaArg).ToNumber(); a.Type != ArgNumber {
|
||||
return a
|
||||
}
|
||||
}
|
||||
if argsList.Len() == 5 {
|
||||
if b = argsList.Back().Value.(formulaArg).ToNumber(); b.Type != ArgNumber {
|
||||
return b
|
||||
}
|
||||
}
|
||||
if x.Number < a.Number || x.Number > b.Number {
|
||||
return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM)
|
||||
}
|
||||
if a.Number == b.Number {
|
||||
return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM)
|
||||
}
|
||||
return newNumberFormulaArg(getBetaDist((x.Number-a.Number)/(b.Number-a.Number), alpha.Number, beta.Number))
|
||||
}
|
||||
|
||||
// d1mach returns double precision real machine constants.
|
||||
func d1mach(i int) float64 {
|
||||
arr := []float64{
|
||||
|
|
@ -5603,6 +5856,32 @@ func (fn *formulaFuncs) CHIDIST(argsList *list.List) formulaArg {
|
|||
return newNumberFormulaArg(1 - (incompleteGamma(degress.Number/2, x.Number/2) / math.Gamma(degress.Number/2)))
|
||||
}
|
||||
|
||||
// CHIINV function calculates the inverse of the right-tailed probability of
|
||||
// the Chi-Square Distribution. The syntax of the function is:
|
||||
//
|
||||
// CHIINV(probability,degrees_freedom)
|
||||
//
|
||||
func (fn *formulaFuncs) CHIINV(argsList *list.List) formulaArg {
|
||||
if argsList.Len() != 2 {
|
||||
return newErrorFormulaArg(formulaErrorVALUE, "CHIINV requires 2 numeric arguments")
|
||||
}
|
||||
probability := argsList.Front().Value.(formulaArg).ToNumber()
|
||||
if probability.Type != ArgNumber {
|
||||
return probability
|
||||
}
|
||||
if probability.Number <= 0 || probability.Number > 1 {
|
||||
return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM)
|
||||
}
|
||||
degress := argsList.Back().Value.(formulaArg).ToNumber()
|
||||
if degress.Type != ArgNumber {
|
||||
return degress
|
||||
}
|
||||
if degress.Number < 1 {
|
||||
return newErrorFormulaArg(formulaErrorNUM, formulaErrorNUM)
|
||||
}
|
||||
return newNumberFormulaArg(gammainv(1-probability.Number, 0.5*degress.Number, 2.0))
|
||||
}
|
||||
|
||||
// confidence is an implementation of the formula functions CONFIDENCE and
|
||||
// CONFIDENCE.NORM.
|
||||
func (fn *formulaFuncs) confidence(name string, argsList *list.List) formulaArg {
|
||||
|
|
@ -6511,7 +6790,6 @@ func (fn *formulaFuncs) LOGNORMdotDIST(argsList *list.List) formulaArg {
|
|||
}
|
||||
return newNumberFormulaArg((1 / (math.Sqrt(2*math.Pi) * stdDev.Number * x.Number)) *
|
||||
math.Exp(0-(math.Pow((math.Log(x.Number)-mean.Number), 2)/(2*math.Pow(stdDev.Number, 2)))))
|
||||
|
||||
}
|
||||
|
||||
// LOGNORMDIST function calculates the Cumulative Log-Normal Distribution
|
||||
|
|
@ -10812,7 +11090,8 @@ func (fn *formulaFuncs) prepareXlookupArgs(argsList *list.List) formulaArg {
|
|||
|
||||
// xlookup is an implementation of the formula function XLOOKUP.
|
||||
func (fn *formulaFuncs) xlookup(lookupRows, lookupCols, returnArrayRows, returnArrayCols, matchIdx int,
|
||||
condition1, condition2, condition3, condition4 bool, returnArray formulaArg) formulaArg {
|
||||
condition1, condition2, condition3, condition4 bool, returnArray formulaArg,
|
||||
) formulaArg {
|
||||
result := [][]formulaArg{}
|
||||
for rowIdx, row := range returnArray.Matrix {
|
||||
for colIdx, cell := range row {
|
||||
|
|
|
|||
48
calc_test.go
48
calc_test.go
|
|
@ -784,6 +784,20 @@ func TestCalcCellValue(t *testing.T) {
|
|||
"=AVERAGEA(A1)": "1",
|
||||
"=AVERAGEA(A1:A2)": "1.5",
|
||||
"=AVERAGEA(D2:F9)": "12671.375",
|
||||
// BETADIST
|
||||
"=BETADIST(0.4,4,5)": "0.4059136",
|
||||
"=BETADIST(0.4,4,5,0,1)": "0.4059136",
|
||||
"=BETADIST(0.4,4,5,0.4,1)": "0",
|
||||
"=BETADIST(1,2,2,1,3)": "0",
|
||||
"=BETADIST(0.4,4,5,0.2,0.4)": "1",
|
||||
"=BETADIST(0.4,4,1)": "0.0256",
|
||||
"=BETADIST(0.4,1,5)": "0.92224",
|
||||
"=BETADIST(3,4,6,2,4)": "0.74609375",
|
||||
"=BETADIST(0.4,2,100)": "1",
|
||||
"=BETADIST(0.75,3,4)": "0.96240234375",
|
||||
"=BETADIST(0.2,0.7,4)": "0.71794309318323",
|
||||
"=BETADIST(0.01,3,4)": "1.955359E-05",
|
||||
"=BETADIST(0.75,130,140)": "1",
|
||||
// BETAINV
|
||||
"=BETAINV(0.2,4,5,0,1)": "0.303225844664082",
|
||||
// BETA.INV
|
||||
|
|
@ -791,6 +805,11 @@ func TestCalcCellValue(t *testing.T) {
|
|||
// CHIDIST
|
||||
"=CHIDIST(0.5,3)": "0.918891411654676",
|
||||
"=CHIDIST(8,3)": "0.0460117056892315",
|
||||
// CHIINV
|
||||
"=CHIINV(0.5,1)": "0.454936423119572",
|
||||
"=CHIINV(0.75,1)": "0.101531044267622",
|
||||
"=CHIINV(0.1,2)": "4.605170185988088",
|
||||
"=CHIINV(0.8,2)": "0.446287102628419",
|
||||
// CONFIDENCE
|
||||
"=CONFIDENCE(0.05,0.07,100)": "0.0137197479028414",
|
||||
// CONFIDENCE.NORM
|
||||
|
|
@ -2322,6 +2341,19 @@ func TestCalcCellValue(t *testing.T) {
|
|||
"=AVERAGE(H1)": "AVERAGE divide by zero",
|
||||
// AVERAGEA
|
||||
"=AVERAGEA(H1)": "AVERAGEA divide by zero",
|
||||
// BETADIST
|
||||
"=BETADIST()": "BETADIST requires at least 3 arguments",
|
||||
"=BETADIST(0.4,4,5,0,1,0)": "BETADIST requires at most 5 arguments",
|
||||
"=BETADIST(\"\",4,5,0,1)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=BETADIST(0.4,\"\",5,0,1)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=BETADIST(0.4,4,\"\",0,1)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=BETADIST(0.4,4,5,\"\",1)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=BETADIST(0.4,4,5,0,\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=BETADIST(2,4,5,3,1)": "#NUM!",
|
||||
"=BETADIST(2,4,5,0,1)": "#NUM!",
|
||||
"=BETADIST(0.4,0,5,0,1)": "#NUM!",
|
||||
"=BETADIST(0.4,4,0,0,1)": "#NUM!",
|
||||
"=BETADIST(0.4,4,5,0.4,0.4)": "#NUM!",
|
||||
// BETAINV
|
||||
"=BETAINV()": "BETAINV requires at least 3 arguments",
|
||||
"=BETAINV(0.2,4,5,0,1,0)": "BETAINV requires at most 5 arguments",
|
||||
|
|
@ -2357,6 +2389,13 @@ func TestCalcCellValue(t *testing.T) {
|
|||
"=CHIDIST()": "CHIDIST requires 2 numeric arguments",
|
||||
"=CHIDIST(\"\",3)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=CHIDIST(0.5,\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
// CHIINV
|
||||
"=CHIINV()": "CHIINV requires 2 numeric arguments",
|
||||
"=CHIINV(\"\",1)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=CHIINV(0.5,\"\")": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
"=CHIINV(0,1)": "#NUM!",
|
||||
"=CHIINV(2,1)": "#NUM!",
|
||||
"=CHIINV(0.5,0.5)": "#NUM!",
|
||||
// CONFIDENCE
|
||||
"=CONFIDENCE()": "CONFIDENCE requires 3 numeric arguments",
|
||||
"=CONFIDENCE(\"\",0.07,100)": "strconv.ParseFloat: parsing \"\": invalid syntax",
|
||||
|
|
@ -4388,6 +4427,15 @@ func TestGetYearDays(t *testing.T) {
|
|||
}
|
||||
}
|
||||
|
||||
func TestCalcGetBetaHelperContFrac(t *testing.T) {
|
||||
assert.Equal(t, 1.0, getBetaHelperContFrac(1, 0, 1))
|
||||
}
|
||||
|
||||
func TestCalcGetBetaDistPDF(t *testing.T) {
|
||||
assert.Equal(t, 0.0, getBetaDistPDF(0.5, 2000, 3))
|
||||
assert.Equal(t, 0.0, getBetaDistPDF(0, 1, 0))
|
||||
}
|
||||
|
||||
func TestCalcD1mach(t *testing.T) {
|
||||
assert.Equal(t, 0.0, d1mach(6))
|
||||
}
|
||||
|
|
|
|||
Loading…
Reference in New Issue