From afa0a718d5d8d5866f1188a1e702fc6a3bfb4df0 Mon Sep 17 00:00:00 2001 From: "dependabot[bot]" <49699333+dependabot[bot]@users.noreply.github.com> Date: Fri, 21 Aug 2026 17:13:09 +0000 Subject: [PATCH] Bump github.com/montanaflynn/stats from 0.12.2 to 0.12.4 Bumps [github.com/montanaflynn/stats](https://github.com/montanaflynn/stats) from 0.12.2 to 0.12.4. - [Release notes](https://github.com/montanaflynn/stats/releases) - [Changelog](https://github.com/montanaflynn/stats/blob/master/CHANGELOG.md) - [Commits](https://github.com/montanaflynn/stats/compare/v0.12.2...v0.12.4) --- updated-dependencies: - dependency-name: github.com/montanaflynn/stats dependency-version: 0.12.4 dependency-type: direct:production update-type: version-update:semver-patch ... Signed-off-by: dependabot[bot] --- go.mod | 2 +- go.sum | 4 +- .../montanaflynn/stats/CHANGELOG.md | 69 ++++++------------- .../montanaflynn/stats/DOCUMENTATION.md | 64 +++++++++-------- .../stats/geometric_distribution.go | 9 +-- .../github.com/montanaflynn/stats/interp.go | 29 ++++++-- vendor/github.com/montanaflynn/stats/norm.go | 64 ++++++++++++----- .../montanaflynn/stats/percentile.go | 6 +- .../montanaflynn/stats/percentile_of_score.go | 2 +- .../montanaflynn/stats/percentile_weighted.go | 2 +- vendor/github.com/montanaflynn/stats/ttest.go | 2 +- vendor/modules.txt | 2 +- 12 files changed, 145 insertions(+), 110 deletions(-) diff --git a/go.mod b/go.mod index 3925508..6e707e6 100644 --- a/go.mod +++ b/go.mod @@ -6,7 +6,7 @@ require ( github.com/alecthomas/kingpin/v2 v2.4.0 github.com/illumos/go-kstat v0.0.0-20210513183136-173c9b0a9973 github.com/kubeservice-stack/common v1.9.1 - github.com/montanaflynn/stats v0.12.2 + github.com/montanaflynn/stats v0.12.4 github.com/prometheus/client_golang v1.24.0 github.com/prometheus/common v0.70.0 github.com/prometheus/exporter-toolkit v0.17.1 diff --git a/go.sum b/go.sum index ed49b98..ebdb80b 100644 --- a/go.sum +++ b/go.sum @@ -41,8 +41,8 @@ github.com/mdlayher/socket v0.6.0 h1:ScZPaAGyO1icQnbFrhPM8mnXyMu9qukC1K4ZoM2IQKU github.com/mdlayher/socket v0.6.0/go.mod h1:q7vozUAnxSqnjHc12Fik5yUKIzfZ8ITCfMkhOtE9z18= github.com/mdlayher/vsock v1.3.0 h1:bqQfZ1OznI03y6YiXp2sze05RVdzLn/zsfjnjd4+ivI= github.com/mdlayher/vsock v1.3.0/go.mod h1:WsuksavOvwCnV5UqGHUkvAvCy+Dqy81y4goKQTzxxNY= -github.com/montanaflynn/stats v0.12.2 h1:qHR+IveGjTbO+lnrz1nKR+xpIcOtovJ5Xu0cst99h80= -github.com/montanaflynn/stats v0.12.2/go.mod h1:etXPPgVO6n31NxCd9KQUMvCM+ve0ruNzt6R8Bnaayow= +github.com/montanaflynn/stats v0.12.4 h1:amtNRsti20yIhcrkfUJGwoYqBR82jKQFE8SNNYVgGn0= +github.com/montanaflynn/stats v0.12.4/go.mod h1:etXPPgVO6n31NxCd9KQUMvCM+ve0ruNzt6R8Bnaayow= github.com/munnerz/goautoneg v0.0.0-20191010083416-a7dc8b61c822 h1:C3w9PqII01/Oq1c1nUAm88MOHcQC9l5mIlSMApZMrHA= github.com/munnerz/goautoneg v0.0.0-20191010083416-a7dc8b61c822/go.mod h1:+n7T8mK8HuQTcFwEeznm/DIxMOiR9yIdICNftLE1DvQ= github.com/mwitkow/go-conntrack v0.0.0-20190716064945-2f068394615f h1:KUppIJq7/+SVif2QVs3tOP0zanoHgBEVAwHxUSIzRqU= diff --git a/vendor/github.com/montanaflynn/stats/CHANGELOG.md b/vendor/github.com/montanaflynn/stats/CHANGELOG.md index ea3553a..eebf097 100644 --- a/vendor/github.com/montanaflynn/stats/CHANGELOG.md +++ b/vendor/github.com/montanaflynn/stats/CHANGELOG.md @@ -7,16 +7,29 @@ - -## [v0.12.2] - 2026-07-17 + + + +## [v0.12.4] - 2026-08-17 ### Fix -- Regression stability and invalid domains ([#124](https://github.com/montanaflynn/stats/issues/124)) +- Make Interp robust to extreme values and exact knot hits ([#132](https://github.com/montanaflynn/stats/issues/132)) +- Reject NaN percent in PercentileWeighted ([#131](https://github.com/montanaflynn/stats/issues/131)) +- Reject NaN percent in Percentile and PercentileNearestRank ([#130](https://github.com/montanaflynn/stats/issues/130)) +- Interp panics when x or xp contains a NaN ([#129](https://github.com/montanaflynn/stats/issues/129)) - -## [v0.12.1] - 2026-07-16 + +## [v0.12.3] - 2026-08-10 ### Fix -- Stop Entropy from mutating its input slice ([#123](https://github.com/montanaflynn/stats/issues/123)) +- ProbGeom off-by-one that drops the first interval term ([#127](https://github.com/montanaflynn/stats/issues/127)) +- NormIsf ignores loc, NormPpf uncorrected above the median ([#126](https://github.com/montanaflynn/stats/issues/126)) +- Normal tail collapses to 0/-Inf past ~8 sigma ([#125](https://github.com/montanaflynn/stats/issues/125)) + + + +## [v0.12.2] - 2026-07-17 +### Fix +- Regression stability and invalid domains ([#124](https://github.com/montanaflynn/stats/issues/124)) @@ -28,22 +41,6 @@ ## [v0.12.0] - 2026-07-16 - -## [v0.12.0] - 2026-07-16 - - -## [v0.11.0] - 2026-07-13 -### Add -- Add Interp for piecewise-linear interpolation ([#121](https://github.com/montanaflynn/stats/issues/121)) -- Add Histogram with equal-width bins ([#120](https://github.com/montanaflynn/stats/issues/120)) -- Add KendallTau rank correlation coefficient ([#119](https://github.com/montanaflynn/stats/issues/119)) -- Add SEM, RMS, Product, and PercentileOfScore ([#118](https://github.com/montanaflynn/stats/issues/118)) -- Add MovingMedian, MovingMin, MovingMax, MovingSum, and EWMA ([#117](https://github.com/montanaflynn/stats/issues/117)) -- Add TrimmedMean and Winsorize robust statistics ([#116](https://github.com/montanaflynn/stats/issues/116)) -- Add Kurtosis, PopulationKurtosis, and SampleKurtosis ([#115](https://github.com/montanaflynn/stats/issues/115)) -- Add Clip and Rescale elementwise transforms ([#114](https://github.com/montanaflynn/stats/issues/114)) - - ## [v0.11.0] - 2026-07-13 ### Add @@ -77,26 +74,6 @@ - Correct AutoCorrelation lag handling ([#83](https://github.com/montanaflynn/stats/issues/83)) ([#95](https://github.com/montanaflynn/stats/issues/95)) - -## [v0.10.0] - 2026-07-10 -### Add -- Add MovingAverage and MovingStdDev ([#112](https://github.com/montanaflynn/stats/issues/112)) -- Add ZScore and Rank functions ([#111](https://github.com/montanaflynn/stats/issues/111)) -- Add WeightedMean and CoefficientOfVariation ([#110](https://github.com/montanaflynn/stats/issues/110)) -- Add ArgMax, ArgMin and Range functions ([#109](https://github.com/montanaflynn/stats/issues/109)) -- Add CumulativeProduct, CumulativeMax and CumulativeMin ([#108](https://github.com/montanaflynn/stats/issues/108)) -- Add Diff and PercentChange functions ([#107](https://github.com/montanaflynn/stats/issues/107)) -- Add weighted percentile function ([#102](https://github.com/montanaflynn/stats/issues/102)) -- Add NormSample function for normal distribution sampling ([#100](https://github.com/montanaflynn/stats/issues/100)) -- Add Z-test and T-test functions ([#99](https://github.com/montanaflynn/stats/issues/99)) -- Add Spearman rank correlation function ([#98](https://github.com/montanaflynn/stats/issues/98)) - -### Fix -- Stabilize GeometricMean and add input validation -- Use math.Round to avoid ARM64 FMA fusion miscompile ([#97](https://github.com/montanaflynn/stats/issues/97)) -- Correct AutoCorrelation lag handling ([#83](https://github.com/montanaflynn/stats/issues/83)) ([#95](https://github.com/montanaflynn/stats/issues/95)) - - ## [v0.9.0] - 2026-03-24 ### Add @@ -638,16 +615,14 @@ - Merge pull request [#4](https://github.com/montanaflynn/stats/issues/4) from saromanov/sample -[Unreleased]: https://github.com/montanaflynn/stats/compare/v0.12.2...HEAD +[Unreleased]: https://github.com/montanaflynn/stats/compare/v0.12.4...HEAD +[v0.12.4]: https://github.com/montanaflynn/stats/compare/v0.12.3...v0.12.4 +[v0.12.3]: https://github.com/montanaflynn/stats/compare/v0.12.2...v0.12.3 [v0.12.2]: https://github.com/montanaflynn/stats/compare/v0.12.1...v0.12.2 [v0.12.1]: https://github.com/montanaflynn/stats/compare/v0.12.0...v0.12.1 -[v0.12.1]: https://github.com/montanaflynn/stats/compare/v0.12.0...v0.12.1 [v0.12.0]: https://github.com/montanaflynn/stats/compare/v0.11.0...v0.12.0 -[v0.12.0]: https://github.com/montanaflynn/stats/compare/v0.11.0...v0.12.0 -[v0.11.0]: https://github.com/montanaflynn/stats/compare/v0.10.0...v0.11.0 [v0.11.0]: https://github.com/montanaflynn/stats/compare/v0.10.0...v0.11.0 [v0.10.0]: https://github.com/montanaflynn/stats/compare/v0.9.0...v0.10.0 -[v0.10.0]: https://github.com/montanaflynn/stats/compare/v0.9.0...v0.10.0 [v0.9.0]: https://github.com/montanaflynn/stats/compare/v0.8.2...v0.9.0 [v0.8.2]: https://github.com/montanaflynn/stats/compare/v0.8.1...v0.8.2 [v0.8.1]: https://github.com/montanaflynn/stats/compare/v0.8.0...v0.8.1 diff --git a/vendor/github.com/montanaflynn/stats/DOCUMENTATION.md b/vendor/github.com/montanaflynn/stats/DOCUMENTATION.md index 0dc1c0d..ee153d6 100644 --- a/vendor/github.com/montanaflynn/stats/DOCUMENTATION.md +++ b/vendor/github.com/montanaflynn/stats/DOCUMENTATION.md @@ -83,7 +83,7 @@ MIT License Copyright (c) 2014-2026 Montana Flynn (func [ExpGeom](/geometric_distribution.go?s=652:700#L27) +## func [ExpGeom](/geometric_distribution.go?s=816:864#L28) ``` go func ExpGeom(p float64) (exp float64, err error) ``` @@ -525,7 +525,7 @@ InterQuartileRange finds the range between Q1 and Q3 -## func [Interp](/interp.go?s=547:600#L12) +## func [Interp](/interp.go?s=635:688#L16) ``` go func Interp(x, xp, fp Float64Data) ([]float64, error) ``` @@ -536,6 +536,7 @@ fp[len(xp)-1], so no extrapolation is performed. Unlike numpy's interp, which silently returns nonsense for unsorted coordinates, xp must be strictly increasing or ErrBounds is returned. An empty x or xp returns ErrEmptyInput and xp and fp of different lengths return ErrSize. +A NaN in xp returns ErrBounds and a NaN in x gives a NaN in the output. @@ -729,7 +730,7 @@ returned. An empty input returns ErrEmptyInput. -## func [Ncr](/norm.go?s=7623:7645#L245) +## func [Ncr](/norm.go?s=8827:8849#L277) ``` go func Ncr(n, r int) int ``` @@ -747,7 +748,7 @@ For more information please visit: func [NormCdf](/norm.go?s=2065:2124#L58) +## func [NormCdf](/norm.go?s=2034:2093#L59) ``` go func NormCdf(x float64, loc float64, scale float64) float64 ``` @@ -755,7 +756,7 @@ NormCdf is the cumulative distribution function. -## func [NormEntropy](/norm.go?s=6012:6064#L186) +## func [NormEntropy](/norm.go?s=7117:7169#L219) ``` go func NormEntropy(loc float64, scale float64) float64 ``` @@ -763,7 +764,7 @@ NormEntropy is the differential entropy of the RV. -## func [NormFit](/norm.go?s=6297:6336#L193) +## func [NormFit](/norm.go?s=7402:7441#L226) ``` go func NormFit(data []float64) [2]float64 ``` @@ -773,7 +774,7 @@ Returns array of Mean followed by Standard Deviation. -## func [NormInterval](/norm.go?s=7215:7286#L227) +## func [NormInterval](/norm.go?s=8320:8391#L260) ``` go func NormInterval(alpha float64, loc float64, scale float64) [2]float64 ``` @@ -781,15 +782,15 @@ NormInterval finds endpoints of the range that contains alpha percent of the dis -## func [NormIsf](/norm.go?s=4569:4632#L143) +## func [NormIsf](/norm.go?s=5589:5648#L177) ``` go -func NormIsf(p float64, loc float64, scale float64) (x float64) +func NormIsf(p float64, loc float64, scale float64) float64 ``` NormIsf is the inverse survival function (inverse of sf). -## func [NormLogCdf](/norm.go?s=2255:2317#L63) +## func [NormLogCdf](/norm.go?s=2218:2280#L64) ``` go func NormLogCdf(x float64, loc float64, scale float64) float64 ``` @@ -805,7 +806,7 @@ NormLogPdf is the log of the probability density function. -## func [NormLogSf](/norm.go?s=2662:2723#L73) +## func [NormLogSf](/norm.go?s=2664:2725#L78) ``` go func NormLogSf(x float64, loc float64, scale float64) float64 ``` @@ -813,7 +814,7 @@ NormLogSf is the log of the survival function. -## func [NormMean](/norm.go?s=6799:6848#L212) +## func [NormMean](/norm.go?s=7904:7953#L245) ``` go func NormMean(loc float64, scale float64) float64 ``` @@ -821,7 +822,7 @@ NormMean is the mean/expected value of the distribution. -## func [NormMedian](/norm.go?s=6670:6721#L207) +## func [NormMedian](/norm.go?s=7775:7826#L240) ``` go func NormMedian(loc float64, scale float64) float64 ``` @@ -829,7 +830,7 @@ NormMedian is the median of the distribution. -## func [NormMoment](/norm.go?s=4933:4991#L152) +## func [NormMoment](/norm.go?s=6038:6096#L185) ``` go func NormMoment(n int, loc float64, scale float64) float64 ``` @@ -846,7 +847,7 @@ NormPdf is the probability density function. -## func [NormPpf](/norm.go?s=3093:3156#L81) +## func [NormPpf](/norm.go?s=3828:3891#L107) ``` go func NormPpf(p float64, loc float64, scale float64) (x float64) ``` @@ -875,7 +876,7 @@ with the given mean (loc) and standard deviation (scale). -## func [NormSf](/norm.go?s=2489:2547#L68) +## func [NormSf](/norm.go?s=2498:2556#L73) ``` go func NormSf(x float64, loc float64, scale float64) float64 ``` @@ -883,7 +884,7 @@ NormSf is the survival function (also defined as 1 - cdf, but sf is sometimes mo -## func [NormStats](/norm.go?s=5516:5584#L168) +## func [NormStats](/norm.go?s=6621:6689#L201) ``` go func NormStats(loc float64, scale float64, moments string) []float64 ``` @@ -894,7 +895,7 @@ Returns array of m v s k in that order. -## func [NormStd](/norm.go?s=7053:7101#L222) +## func [NormStd](/norm.go?s=8158:8206#L255) ``` go func NormStd(loc float64, scale float64) float64 ``` @@ -902,7 +903,7 @@ NormStd is the standard deviation of the distribution. -## func [NormVar](/norm.go?s=6914:6962#L217) +## func [NormVar](/norm.go?s=8019:8067#L250) ``` go func NormVar(loc float64, scale float64) float64 ``` @@ -953,7 +954,7 @@ Algorithm (for percent p and sorted data of length n): -## func [PercentileNearestRank](/percentile.go?s=1382:1476#L55) +## func [PercentileNearestRank](/percentile.go?s=1405:1499#L55) ``` go func PercentileNearestRank(input Float64Data, percent float64) (percentile float64, err error) ``` @@ -970,7 +971,7 @@ relative to a slice of floats, defined as the percentage of values strictly below the score plus half the percentage of values equal to the score. The result is between 0 and 100. This matches the behavior of Python's -scipy.stats.percentileofscore with kind="rank". +scipy.stats.percentileofscore with kind="mean". @@ -1268,7 +1269,7 @@ the same result as Mean. -## func [VarGeom](/geometric_distribution.go?s=885:933#L37) +## func [VarGeom](/geometric_distribution.go?s=1049:1097#L38) ``` go func VarGeom(p float64) (exp float64, err error) ``` @@ -2196,25 +2197,30 @@ Series is a container for a series of data -### func [ExponentialRegression](/regression.go?s=1061:1129#L49) +### func [ExponentialRegression](/regression.go?s=1269:1337#L54) ``` go func ExponentialRegression(s Series) (regressions Series, err error) ``` -ExponentialRegression returns an exponential regression on data series +ExponentialRegression returns an exponential regression on data series. +A non-positive Y value returns ErrYCoord, and a series without at least two +distinct X values returns ErrBounds. -### func [LinearRegression](/regression.go?s=262:325#L14) +### func [LinearRegression](/regression.go?s=333:396#L15) ``` go func LinearRegression(s Series) (regressions Series, err error) ``` -LinearRegression finds the least squares linear regression on data series +LinearRegression finds the least squares linear regression on data series. +A series without at least two distinct X values returns ErrBounds. -### func [LogarithmicRegression](/regression.go?s=1875:1943#L84) +### func [LogarithmicRegression](/regression.go?s=2368:2436#L98) ``` go func LogarithmicRegression(s Series) (regressions Series, err error) ``` -LogarithmicRegression returns an logarithmic regression on data series +LogarithmicRegression returns a logarithmic regression on data series. +A non-positive X value or a series without at least two distinct X values +returns ErrBounds. diff --git a/vendor/github.com/montanaflynn/stats/geometric_distribution.go b/vendor/github.com/montanaflynn/stats/geometric_distribution.go index db785dd..f1dfebe 100644 --- a/vendor/github.com/montanaflynn/stats/geometric_distribution.go +++ b/vendor/github.com/montanaflynn/stats/geometric_distribution.go @@ -12,14 +12,15 @@ func ProbGeom(a int, b int, p float64) (prob float64, err error) { return math.NaN(), ErrBounds } - prob = 0 q := 1 - p // probability of failure - for k := a + 1; k <= b; k++ { - prob = prob + p*math.Pow(q, float64(k-1)) + if a == b { + return p * math.Pow(q, float64(a-1)), nil } - return prob, nil + // closed form of the sum p*q^(k-1) over k = a..b; expm1/log1p keep + // 1-q^n accurate where direct subtraction would cancel + return math.Pow(q, float64(a-1)) * -math.Expm1(float64(b-a+1)*math.Log1p(-p)), nil } // ProbGeom generates the expectation or average number of trials diff --git a/vendor/github.com/montanaflynn/stats/interp.go b/vendor/github.com/montanaflynn/stats/interp.go index 86a4c3a..df75063 100644 --- a/vendor/github.com/montanaflynn/stats/interp.go +++ b/vendor/github.com/montanaflynn/stats/interp.go @@ -1,6 +1,9 @@ package stats -import "sort" +import ( + "math" + "sort" +) // Interp calculates the one-dimensional piecewise-linear interpolant to a // function with given discrete data points (xp, fp), evaluated at each x. @@ -9,6 +12,7 @@ import "sort" // which silently returns nonsense for unsorted coordinates, xp must be // strictly increasing or ErrBounds is returned. An empty x or xp returns // ErrEmptyInput and xp and fp of different lengths return ErrSize. +// A NaN in xp returns ErrBounds and a NaN in x gives a NaN in the output. func Interp(x, xp, fp Float64Data) ([]float64, error) { if x.Len() == 0 || xp.Len() == 0 { @@ -19,8 +23,9 @@ func Interp(x, xp, fp Float64Data) ([]float64, error) { return nil, ErrSize } - for i := 1; i < xp.Len(); i++ { - if xp[i] <= xp[i-1] { + // NaN loses every comparison, so the ordering check can't catch it + for i := 0; i < xp.Len(); i++ { + if math.IsNaN(xp[i]) || (i > 0 && xp[i] <= xp[i-1]) { return nil, ErrBounds } } @@ -29,6 +34,8 @@ func Interp(x, xp, fp Float64Data) ([]float64, error) { for i, xv := range x { switch { + case math.IsNaN(xv): + output[i] = math.NaN() case xv <= xp[0]: output[i] = fp[0] case xv >= xp[xp.Len()-1]: @@ -37,8 +44,22 @@ func Interp(x, xp, fp Float64Data) ([]float64, error) { // The first index with xp[j] >= xv, which the clamping // above guarantees is within [1, len(xp)-1] j := sort.SearchFloat64s(xp, xv) + if xv == xp[j] { + // An exact knot hit returns fp[j] exactly, with no + // interpolation arithmetic that could lose precision + output[i] = fp[j] + continue + } t := (xv - xp[j-1]) / (xp[j] - xp[j-1]) - output[i] = fp[j-1] + t*(fp[j]-fp[j-1]) + if math.IsInf(xp[j]-xp[j-1], 1) { + // The knot spacing overflows float64, so halve each + // term before dividing; halving is exact for the huge + // values that make an overflowing difference possible + t = (xv/2 - xp[j-1]/2) / (xp[j]/2 - xp[j-1]/2) + } + // The symmetric form stays finite for any finite fp where + // fp[j]-fp[j-1] would overflow, since 0 < t < 1 here + output[i] = fp[j-1]*(1-t) + fp[j]*t } } diff --git a/vendor/github.com/montanaflynn/stats/norm.go b/vendor/github.com/montanaflynn/stats/norm.go index 620c5ac..0d0fd08 100644 --- a/vendor/github.com/montanaflynn/stats/norm.go +++ b/vendor/github.com/montanaflynn/stats/norm.go @@ -51,27 +51,53 @@ func NormPdf(x float64, loc float64, scale float64) float64 { // NormLogPdf is the log of the probability density function. func NormLogPdf(x float64, loc float64, scale float64) float64 { - return math.Log((math.Pow(math.E, -(math.Pow(x-loc, 2))/(2*math.Pow(scale, 2)))) / (scale * math.Sqrt(2*math.Pi))) + z := (x - loc) / scale + return -0.5*z*z - math.Log(scale) - 0.5*math.Log(2*math.Pi) } // NormCdf is the cumulative distribution function. func NormCdf(x float64, loc float64, scale float64) float64 { - return 0.5 * (1 + math.Erf((x-loc)/(scale*math.Sqrt(2)))) + return 0.5 * math.Erfc(-(x-loc)/(scale*math.Sqrt2)) } // NormLogCdf is the log of the cumulative distribution function. func NormLogCdf(x float64, loc float64, scale float64) float64 { - return math.Log(0.5 * (1 + math.Erf((x-loc)/(scale*math.Sqrt(2))))) + z := (x - loc) / scale + if z > 0 { + return math.Log1p(-0.5 * math.Erfc(z/math.Sqrt2)) + } + return normLogTail(-z) } // NormSf is the survival function (also defined as 1 - cdf, but sf is sometimes more accurate). func NormSf(x float64, loc float64, scale float64) float64 { - return 1 - 0.5*(1+math.Erf((x-loc)/(scale*math.Sqrt(2)))) + return 0.5 * math.Erfc((x-loc)/(scale*math.Sqrt2)) } // NormLogSf is the log of the survival function. func NormLogSf(x float64, loc float64, scale float64) float64 { - return math.Log(1 - 0.5*(1+math.Erf((x-loc)/(scale*math.Sqrt(2))))) + z := (x - loc) / scale + if z < 0 { + return math.Log1p(-0.5 * math.Erfc(-z/math.Sqrt2)) + } + return normLogTail(z) +} + +// normSmallestNormal is the smallest positive normal float64; below it math.Erfc +// keeps only a handful of significant bits. +const normSmallestNormal = 2.2250738585072014e-308 + +// normLogTail returns log(sf(z)) for z >= 0. +func normLogTail(z float64) float64 { + if q := 0.5 * math.Erfc(z/math.Sqrt2); q >= normSmallestNormal { + return math.Log(q) + } + // math.Erfc has decayed into the subnormals, so switch to the Mills ratio + // expansion sf(z) = pdf(z)/z * (1 - 1/z^2 + 3/z^4 - 15/z^6 + 105/z^8 - ...), + // whose first dropped term is below 1e-12 this far out. + r := 1 / (z * z) + return -0.5*z*z - math.Log(z) - 0.5*math.Log(2*math.Pi) + + math.Log1p(r*(-1+r*(3+r*(-15+r*105)))) } // NormPpf is the point percentile function. @@ -132,7 +158,15 @@ func NormPpf(p float64, loc float64, scale float64) (x float64) { (((((b1*r+b2)*r+b3)*r+b4)*r+b5)*r + 1) } - e := 0.5*math.Erfc(-x/math.Sqrt2) - p + // Halley correction on cdf(x)-p. Above the median cdf(x) and p have both + // already rounded to 1, so the difference is taken between the survival + // functions instead; 1-p is exact for p >= 0.5. + var e float64 + if p > 0.5 { + e = (1 - p) - 0.5*math.Erfc(x/math.Sqrt2) + } else { + e = 0.5*math.Erfc(-x/math.Sqrt2) - p + } u := e * math.Sqrt(2*math.Pi) * math.Exp(x*x/2) x = x - u/(1+x*u/2) @@ -140,11 +174,10 @@ func NormPpf(p float64, loc float64, scale float64) (x float64) { } // NormIsf is the inverse survival function (inverse of sf). -func NormIsf(p float64, loc float64, scale float64) (x float64) { - if -NormPpf(p, loc, scale) == 0 { - return 0 - } - return -NormPpf(p, loc, scale) +func NormIsf(p float64, loc float64, scale float64) float64 { + // isf(p) == ppf(1-p), reached by reflecting the standard normal so that + // loc stays out of the negation and 1-p is never formed. + return loc - scale*NormPpf(p, 0, 1) } // NormMoment approximates the non-central (raw) moment of order n. @@ -225,11 +258,10 @@ func NormStd(loc float64, scale float64) float64 { // NormInterval finds endpoints of the range that contains alpha percent of the distribution. func NormInterval(alpha float64, loc float64, scale float64) [2]float64 { - q1 := (1.0 - alpha) / 2 - q2 := (1.0 + alpha) / 2 - a := NormPpf(q1, loc, scale) - b := NormPpf(q2, loc, scale) - return [2]float64{a, b} + // Derive both endpoints from the lower tail: (1+alpha)/2 rounds to 1 once + // alpha is within an ulp of it, which sends the upper endpoint to +Inf. + z := NormPpf((1.0-alpha)/2, 0, 1) + return [2]float64{loc + scale*z, loc - scale*z} } // factorial is the naive factorial algorithm. diff --git a/vendor/github.com/montanaflynn/stats/percentile.go b/vendor/github.com/montanaflynn/stats/percentile.go index 5bb4d7b..a7f27d3 100644 --- a/vendor/github.com/montanaflynn/stats/percentile.go +++ b/vendor/github.com/montanaflynn/stats/percentile.go @@ -27,7 +27,7 @@ func Percentile(input Float64Data, percent float64) (percentile float64, err err return input[0], nil } - if percent <= 0 || percent > 100 { + if math.IsNaN(percent) || percent <= 0 || percent > 100 { return math.NaN(), BoundsErr } @@ -62,8 +62,8 @@ func PercentileNearestRank(input Float64Data, percent float64) (percentile float return math.NaN(), EmptyInputErr } - // Return error for less than 0 or greater than 100 percentages - if percent < 0 || percent > 100 { + // Return error for NaN, less than 0, or greater than 100 percentages + if math.IsNaN(percent) || percent < 0 || percent > 100 { return math.NaN(), BoundsErr } diff --git a/vendor/github.com/montanaflynn/stats/percentile_of_score.go b/vendor/github.com/montanaflynn/stats/percentile_of_score.go index a958c95..1e46e0c 100644 --- a/vendor/github.com/montanaflynn/stats/percentile_of_score.go +++ b/vendor/github.com/montanaflynn/stats/percentile_of_score.go @@ -7,7 +7,7 @@ import "math" // values strictly below the score plus half the percentage of // values equal to the score. The result is between 0 and 100. // This matches the behavior of Python's -// scipy.stats.percentileofscore with kind="rank". +// scipy.stats.percentileofscore with kind="mean". func PercentileOfScore(input Float64Data, score float64) (float64, error) { if input.Len() == 0 { return math.NaN(), ErrEmptyInput diff --git a/vendor/github.com/montanaflynn/stats/percentile_weighted.go b/vendor/github.com/montanaflynn/stats/percentile_weighted.go index beb09a0..992f93c 100644 --- a/vendor/github.com/montanaflynn/stats/percentile_weighted.go +++ b/vendor/github.com/montanaflynn/stats/percentile_weighted.go @@ -26,7 +26,7 @@ func PercentileWeighted(data, weights Float64Data, percent float64) (percentile return math.NaN(), ErrSize } - if percent <= 0 || percent > 100 { + if math.IsNaN(percent) || percent <= 0 || percent > 100 { return math.NaN(), ErrBounds } diff --git a/vendor/github.com/montanaflynn/stats/ttest.go b/vendor/github.com/montanaflynn/stats/ttest.go index f36ea32..728014a 100644 --- a/vendor/github.com/montanaflynn/stats/ttest.go +++ b/vendor/github.com/montanaflynn/stats/ttest.go @@ -76,7 +76,7 @@ func regIncBeta(a, b, x float64) float64 { } lbeta := lgammaBeta(a, b) - front := math.Exp(math.Log(x)*a + math.Log(1-x)*b - lbeta) / a + front := math.Exp(math.Log(x)*a+math.Log(1-x)*b-lbeta) / a // Use Lentz's continued fraction algorithm f := 1.0 diff --git a/vendor/modules.txt b/vendor/modules.txt index 8e72978..3107847 100644 --- a/vendor/modules.txt +++ b/vendor/modules.txt @@ -55,7 +55,7 @@ github.com/mdlayher/socket # github.com/mdlayher/vsock v1.3.0 ## explicit; go 1.25.0 github.com/mdlayher/vsock -# github.com/montanaflynn/stats v0.12.2 +# github.com/montanaflynn/stats v0.12.4 ## explicit; go 1.13 github.com/montanaflynn/stats # github.com/munnerz/goautoneg v0.0.0-20191010083416-a7dc8b61c822