From d4a8fc5b884c0b4f1a56f843d34cf59ba2f435bc Mon Sep 17 00:00:00 2001 From: Masahiro Nagano Date: Mon, 12 Oct 2026 17:09:48 +0900 Subject: [PATCH] feat: add population and sample standard deviation calculations --- README.md | 40 ++++++++++++++++++++-- sorted.go | 56 +++++++++++++++++++++++++++++++ sorted_test.go | 91 ++++++++++++++++++++++++++++++++++++++++++++++++++ 3 files changed, 185 insertions(+), 2 deletions(-) diff --git a/README.md b/README.md index 745cfec..5d5881a 100644 --- a/README.md +++ b/README.md @@ -1,8 +1,8 @@ # sampdo -`sampdo` is a lightweight statistical aggregation library for Go. It collects a large number of `float64` samples and efficiently computes minimum, maximum, mean, median, and arbitrary percentiles from a sorted view. For non-negative values it uses radix sort, which is faster than a comparison sort. +`sampdo` is a lightweight statistical aggregation library for Go. It collects a large number of `float64` samples and efficiently computes minimum, maximum, mean, median, population/sample standard deviation, and arbitrary percentiles from a sorted view. For non-negative values it uses radix sort, which is faster than a comparison sort. -`sampdo` は Go 用の軽量な数値集計ライブラリです。大量の `float64` サンプルに対して、ソート済みの状態から最小値・最大値・平均値・中央値・任意のパーセンタイルを効率的に取得できます。非負数に対しては基数ソートを使い、比較ソートより高速に集計できます。 +`sampdo` は Go 用の軽量な数値集計ライブラリです。大量の `float64` サンプルに対して、ソート済みの状態から最小値・最大値・平均値・中央値・母標準偏差・標本標準偏差・任意のパーセンタイルを効率的に取得できます。非負数に対しては基数ソートを使い、比較ソートより高速に集計できます。 ## Features / 特徴 @@ -107,9 +107,45 @@ P99 : 9.91 | `Min()` | Smallest value; returns `0` when empty. | 最小値。空の場合は `0`。 | | `Max()` | Largest value; returns `0` when empty. | 最大値。空の場合は `0`。 | | `Mean()` | Arithmetic mean; sums in input order. Returns `0` when empty. | 算術平均。入力順に合計し、空の場合は `0`。 | +| `PopulationStdDev()` | Population standard deviation, divisor n. | 母標準偏差。分母は n。 | +| `SampleStdDev()` | Sample standard deviation, divisor n-1; requires two points. | 標本標準偏差。分母は n-1。2件以上が必要。 | | `Median()` | 50th percentile; averages the two middle values for an even count. | 50 パーセンタイル。要素数が偶数の場合は中央 2 値の平均。 | | `Percentile(p)` | Returns the `p`th percentile (0–100) using linear interpolation; returns `0` for out-of-range `p`. | `p` パーセンタイル（0〜100）を線形補間で返す。範囲外は `0`。 | +### Standard deviation / 標準偏差 + +```go +population, err := sorted.PopulationStdDev() +if err != nil { + log.Fatal(err) +} +sample, err := sorted.SampleStdDev() +if err != nil { + log.Fatal(err) +} +fmt.Printf("Population: %.6f, Sample: %.6f\n", population, sample) +// For samples 1 through 10: Population: 2.872281, Sample: 3.027650 +``` + +`PopulationStdDev()` computes `sqrt(sum((x - mean)^2) / n)`; +`SampleStdDev()` computes `sqrt(sum((x - mean)^2) / (n - 1))`. +Both return `(float64, error)` in the same units as the input, using O(n) time +and O(1) extra space per call. Calculations center and scale the values to reduce +rounding errors and avoid intermediate overflow/underflow. + +`PopulationStdDev()` は母標準偏差（分母 `n`）、`SampleStdDev()` は標本標準偏差 +（分母 `n-1`、不偏分散の平方根）を返します。戻り値は `(float64, error)` で、 +標準偏差の単位は入力と同じです。呼び出しごとに O(n) 時間・O(1) 追加メモリで計算します。 +値の中心を移動してスケールを調整し、丸め誤差や中間計算のオーバーフロー・アンダーフローを抑えます。 + +- Empty views return `0, ErrNoPoints` for both methods. + 空の `Sorted` は、どちらも `0, ErrNoPoints` を返します。 +- A single finite point has population standard deviation `0`; sample standard + deviation returns `0, ErrInsufficientPoints` for any single point. + 有限値1件の母標準偏差は `0`。標本標準偏差は値によらず1件では `0, ErrInsufficientPoints` です。 +- With sufficient points, any infinity yields `NaN, nil`. + 必要な件数を満たしていて無限大を含む場合は `NaN, nil` を返します。 + ### How percentiles are calculated / パーセンタイルの計算方法 `Percentile(p)` computes the position `k` with the formula below and linearly interpolates between the two surrounding points. diff --git a/sorted.go b/sorted.go index 4258724..b35954f 100644 --- a/sorted.go +++ b/sorted.go @@ -8,6 +8,9 @@ import ( var ErrNoPoints = errors.New("no points available") var ErrInvalidPercentile = errors.New("invalid percentile") +// ErrInsufficientPoints indicates that sample standard deviation needs at least two points. +var ErrInsufficientPoints = errors.New("at least two points required") + type Sorted struct { points []float64 sum float64 @@ -46,6 +49,59 @@ func (s *Sorted) Mean() (float64, error) { return sum / float64(len(s.points)), nil } +// PopulationStdDev returns the population standard deviation (divisor n). +// Empty input returns 0 and ErrNoPoints. Input containing infinity returns NaN. +func (s *Sorted) PopulationStdDev() (float64, error) { + return s.stdDev(0) +} + +// SampleStdDev returns the sample standard deviation (divisor n-1). +// Empty input returns 0 and ErrNoPoints; one point returns 0 and +// ErrInsufficientPoints. Otherwise, input containing infinity returns NaN. +func (s *Sorted) SampleStdDev() (float64, error) { + return s.stdDev(1) +} + +func (s *Sorted) stdDev(correction int) (float64, error) { + n := len(s.points) + if n == 0 { + return 0, ErrNoPoints + } + if n <= correction { + return 0, ErrInsufficientPoints + } + low, high := s.points[0], s.points[n-1] + if math.IsInf(low, 0) || math.IsInf(high, 0) { + return math.NaN(), nil + } + if low == high { + return 0, nil + } + // Center before scaling to preserve small differences on a large offset. + // Halving the endpoints first avoids overflow for extreme finite inputs. + center := low/2 + high/2 + scale := math.Max(math.Abs(low-center), math.Abs(high-center)) + // Compensate the sum of normalized offsets; the cached input-order sum + // can overflow or lose the very differences needed by this calculation. + sum, compensation := 0.0, 0.0 + for _, v := range s.points { + y := (v-center)/scale - compensation + next := sum + y + compensation = (next - sum) - y + sum = next + } + mean := sum / float64(n) + squares, compensation := 0.0, 0.0 + for _, v := range s.points { + d := (v-center)/scale - mean + y := d*d - compensation + next := squares + y + compensation = (next - squares) - y + squares = next + } + return scale * math.Sqrt(squares/float64(n-correction)), nil +} + func (s *Sorted) Median() (float64, error) { if len(s.points) == 0 { return 0, ErrNoPoints diff --git a/sorted_test.go b/sorted_test.go index 934275b..9927d06 100644 --- a/sorted_test.go +++ b/sorted_test.go @@ -38,6 +38,97 @@ func TestSortedBasic(t *testing.T) { } +func TestSortedStdDev(t *testing.T) { + tests := []struct { + name string + points []float64 + population, sample float64 + }{ + {"textbook", []float64{2, 4, 4, 4, 5, 5, 7, 9}, 2, math.Sqrt(32.0 / 7)}, + {"negative", []float64{-5, -3, -1}, math.Sqrt(8.0 / 3), 2}, + {"mixed", []float64{-1, 0, 1}, math.Sqrt(2.0 / 3), 1}, + {"constant", []float64{42, 42, 42}, 0, 0}, + {"zeros", []float64{math.Copysign(0, -1), 0}, 0, 0}, + {"large offset", []float64{1e16, 1e16 + 2, 1e16 + 4}, math.Sqrt(8.0 / 3), 2}, + {"large magnitude", []float64{-1e200, 1e200}, 1e200, math.Sqrt2 * 1e200}, + {"small magnitude", []float64{-1e-200, 1e-200}, 1e-200, math.Sqrt2 * 1e-200}, + {"overflowing sum", []float64{math.MaxFloat64, math.MaxFloat64}, 0, 0}, + {"extreme range", []float64{-math.MaxFloat64, math.MaxFloat64}, math.MaxFloat64, math.Inf(1)}, + {"subnormal", []float64{0, 2 * math.SmallestNonzeroFloat64}, math.SmallestNonzeroFloat64, math.SmallestNonzeroFloat64}, + {"positive infinity", []float64{1, math.Inf(1)}, math.NaN(), math.NaN()}, + {"negative infinity", []float64{math.Inf(-1), 1}, math.NaN(), math.NaN()}, + {"both infinities", []float64{math.Inf(-1), math.Inf(1)}, math.NaN(), math.NaN()}, + {"equal infinities", []float64{math.Inf(1), math.Inf(1)}, math.NaN(), math.NaN()}, + } + for _, tt := range tests { + t.Run(tt.name, func(t *testing.T) { + data := New() + require.NoError(t, data.Append(tt.points...)) + sorted, err := data.Sorted() + require.NoError(t, err) + for _, view := range []*Sorted{sorted, {points: slices.Clone(sorted.points)}} { + population, err := view.PopulationStdDev() + require.NoError(t, err) + sample, err := view.SampleStdDev() + require.NoError(t, err) + for i, got := range []float64{population, sample} { + want := []float64{tt.population, tt.sample}[i] + switch { + case math.IsNaN(want): + require.True(t, math.IsNaN(got)) + case want == 0 || math.IsInf(want, 0) || math.Abs(want) < 1e-300: + require.Equal(t, want, got) + default: + require.InEpsilon(t, want, got, 1e-14) + } + } + } + }) + } +} + +func TestSortedStdDevSmallInputs(t *testing.T) { + s := &Sorted{} + got, err := s.PopulationStdDev() + require.ErrorIs(t, err, ErrNoPoints) + require.Zero(t, got) + got, err = s.SampleStdDev() + require.ErrorIs(t, err, ErrNoPoints) + require.Zero(t, got) + for _, point := range []float64{42, math.MaxFloat64, math.Inf(1)} { + s = &Sorted{points: []float64{point}} + got, err = s.PopulationStdDev() + require.NoError(t, err) + if math.IsInf(point, 0) { + require.True(t, math.IsNaN(got)) + } else { + require.Zero(t, got) + } + got, err = s.SampleStdDev() + require.ErrorIs(t, err, ErrInsufficientPoints) + require.Zero(t, got) + } +} + +func TestSortedStdDevFrozenView(t *testing.T) { + data := New() + require.NoError(t, data.Append(3, 1)) + before, err := data.Sorted() + require.NoError(t, err) + require.NoError(t, data.Append(5)) + after, err := data.Sorted() + require.NoError(t, err) + got, err := before.PopulationStdDev() + require.NoError(t, err) + require.Equal(t, 1.0, got) + got, err = before.SampleStdDev() + require.NoError(t, err) + require.Equal(t, math.Sqrt2, got) + got, err = after.SampleStdDev() + require.NoError(t, err) + require.Equal(t, 2.0, got) +} + func TestSortedMedian(t *testing.T) { tests := []struct { name string