Skip to content

Commit a43256b

Browse files
JimmyMAnderssonJimmyMAndersson
committed
fix(quantile)!: correct off-by-one index error in all estimation methods
BREAKING CHANGE: QuantileEstimationMethod.closestOrOddIndexed renamed to .closest. Fixes: #177
1 parent 54c192a commit a43256b

2 files changed

Lines changed: 137 additions & 39 deletions

File tree

‎Sources/StatKit/Descriptive Statistics/Quantile.swift‎

Lines changed: 28 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -10,10 +10,14 @@ public extension Collection {
1010
/// Since quantiles have no meaning on empty collections or probabilities outside of range [0, 1],
1111
/// this method returns a NaN for calls under any such conditions.
1212
///
13-
/// - important: This method has undefined behavior on collections containing elements that are incomparable.
14-
/// For example, an array containing NaN's will produce unpredictable results.
13+
/// Collections containing NaN values will also produce a NaN result.
1514
///
16-
/// - complexity: O(n log n), where n is the length of the collection.
15+
/// Infinite values are handled correctly as long as the chosen estimation method does not require
16+
/// arithmetic between conflicting infinities (e.g. `+∞ − ∞`).
17+
/// If such arithmetic is unavoidable for the requested probability and method, this method returns NaN.
18+
///
19+
/// - important: This method has undefined behavior on collections containing elements that are incomparable
20+
/// for reasons other than NaN (e.g. custom types with a broken `Comparable` implementation).
1721
func quantile<T: Comparable & ConvertibleToReal>(
1822
probability: Double,
1923
of variable: KeyPath<Element, T>,
@@ -23,50 +27,52 @@ public extension Collection {
2327
guard
2428
probability.isFinite,
2529
0 <= probability,
26-
probability <= 1
30+
probability <= 1,
31+
!self.isEmpty,
32+
!self.contains(where: { $0[keyPath: variable].realValue.isNaN })
2733
else { return .signalingNaN }
2834

2935
let ordered = self.sorted { lhs, rhs in
3036
lhs[keyPath: variable] < rhs[keyPath: variable]
3137
}
3238

33-
guard !ordered.isEmpty else { return .signalingNaN }
34-
3539
if probability == 1 {
3640
guard let element = ordered.last else {
37-
fatalError("Could not fetch last element of Sequence.")
41+
fatalError("Could not fetch greatest element of collection.")
3842
}
3943
return element[keyPath: variable].realValue
4044
}
4145

4246
if probability == 0 {
4347
guard let element = ordered.first else {
44-
fatalError("Could not fetch first element of Sequence.")
48+
fatalError("Could not fetch smallest element of collection.")
4549
}
4650
return element[keyPath: variable].realValue
4751
}
4852

4953
switch method {
5054
case .inverseEmpiricalCDF:
51-
let h = Double(ordered.count) * probability + 0.5
52-
let index = Int((h - 0.5).rounded(.up)) - 1
55+
let h = Double(ordered.count - 1) * probability
56+
let index = Int(h.rounded(.up))
5357
return ordered[index][keyPath: variable].realValue
5458

5559
case .averagedInverseEmpiricalCDF:
56-
let h = Double(ordered.count) * probability + 0.5
57-
let firstIndex = Int((h - 0.5).rounded(.up)) - 1
58-
let secondIndex = Int((h + 0.5).rounded(.down)) - 1
59-
return ordered[firstIndex...secondIndex].mean(variable: variable)
60+
let h = Double(ordered.count - 1) * probability
61+
let firstIndex = Int(h.rounded(.down))
62+
let secondIndex = Int(h.rounded(.up))
63+
let firstElement = ordered[firstIndex][keyPath: variable]
64+
let secondElement = ordered[secondIndex][keyPath: variable]
65+
return (firstElement.realValue + secondElement.realValue) / 2
6066

61-
case .closestOrOddIndexed:
62-
let h = Double(ordered.count) * probability
63-
let index = Int(h.rounded(.toNearestOrEven)) - 1
67+
case .closest:
68+
let h = Double(ordered.count - 1) * probability
69+
let index = Int(h.rounded(.toNearestOrEven))
6470
return ordered[index][keyPath: variable].realValue
6571

6672
case .lerpInverseEmpiricalCDF:
67-
let h = Double(ordered.count) * probability
68-
let firstIndex = Int(h.rounded(.down)) - 1
69-
let secondIndex = Int(h.rounded(.up)) - 1
73+
let h = Double(ordered.count - 1) * probability
74+
let firstIndex = Int(h.rounded(.down))
75+
let secondIndex = Int(h.rounded(.up))
7076
let firstElement = ordered[firstIndex][keyPath: variable]
7177
let secondElement = ordered[secondIndex][keyPath: variable]
7278
let difference = secondElement - firstElement
@@ -84,8 +90,8 @@ public enum QuantileEstimationMethod: CaseIterable, Sendable {
8490
case averagedInverseEmpiricalCDF
8591
/// Computes the quantile by rounding to the closest observation.
8692
///
87-
/// In case of a tie, the odd index element will be chosen.
88-
case closestOrOddIndexed
93+
/// In case of a tie, the even index element will be chosen.
94+
case closest
8995
/// Computes the quantile usign the inverse empirical CDF, and linearly interpolates at discontinuities.
9096
case lerpInverseEmpiricalCDF
9197
}

‎Tests/StatKitTests/Descriptive Statistics Tests/QuantileTests.swift‎

Lines changed: 109 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -6,22 +6,66 @@ struct QuantileTests {
66
@Test(
77
"Valid data returns correct quantile",
88
arguments: [
9-
([1, 4, 2, 6, 4, 3], 1, .inverseEmpiricalCDF, 6),
10-
([1, 4, 2, 6, 4, 3], 0, .inverseEmpiricalCDF, 1),
11-
([1, 4, 2, 6, 4, 3], 1, .averagedInverseEmpiricalCDF, 6),
12-
([1, 4, 2, 6, 4, 3], 0, .averagedInverseEmpiricalCDF, 1),
13-
([1, 4, 2, 6, 4, 3], 1, .closestOrOddIndexed, 6),
14-
([1, 4, 2, 6, 4, 3], 0, .closestOrOddIndexed, 1),
15-
([1, 4, 2, 6, 4, 3], 1, .lerpInverseEmpiricalCDF, 6),
16-
([1, 4, 2, 6, 4, 3], 0, .lerpInverseEmpiricalCDF, 1),
17-
([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 0.5, .inverseEmpiricalCDF, 5),
18-
([1, 2, 3, 4, 5, 6, 7, 8, 9], 0.5, .inverseEmpiricalCDF, 5),
19-
([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 0.5, .averagedInverseEmpiricalCDF, 5.5),
20-
([1, 2, 3, 4, 5, 6, 7, 8, 9], 0.5, .averagedInverseEmpiricalCDF, 5),
21-
([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 0.5, .closestOrOddIndexed, 5),
22-
([1, 2, 3, 4, 5, 6, 7, 8, 9], 0.5, .closestOrOddIndexed, 4),
23-
([1, 2, 3, 4, 5, 6, 7, 8, 9, 10], 0.5, .lerpInverseEmpiricalCDF, 5),
24-
([1, 2, 3, 4, 5, 6, 7, 8, 9], 0.5, .lerpInverseEmpiricalCDF, 4.5),
9+
([1, 4, 2, 6, 4, 3], 0, .inverseEmpiricalCDF, 1),
10+
([1, 4, 2, 6, 4, 3], 0.25, .inverseEmpiricalCDF, 3),
11+
([1, 4, 2, 6, 4, 3], 0.5, .inverseEmpiricalCDF, 4),
12+
([1, 4, 2, 6, 4, 3], 0.75, .inverseEmpiricalCDF, 4),
13+
([1, 4, 2, 6, 4, 3], 1, .inverseEmpiricalCDF, 6),
14+
([1, 4, 2, 6, 4, 3], 0, .averagedInverseEmpiricalCDF, 1),
15+
([1, 4, 2, 6, 4, 3], 0.25, .averagedInverseEmpiricalCDF, 2.5),
16+
([1, 4, 2, 6, 4, 3], 0.5, .averagedInverseEmpiricalCDF, 3.5),
17+
([1, 4, 2, 6, 4, 3], 0.75, .averagedInverseEmpiricalCDF, 4),
18+
([1, 4, 2, 6, 4, 3], 1, .averagedInverseEmpiricalCDF, 6),
19+
([1, 4, 2, 6, 4, 3], 0, .closest, 1),
20+
([1, 4, 2, 6, 4, 3], 0.25, .closest, 2),
21+
([1, 4, 2, 6, 4, 3], 0.5, .closest, 3),
22+
([1, 4, 2, 6, 4, 3], 0.75, .closest, 4),
23+
([1, 4, 2, 6, 4, 3], 1, .closest, 6),
24+
([1, 4, 2, 6, 4, 3], 0, .lerpInverseEmpiricalCDF, 1),
25+
([1, 4, 2, 6, 4, 3], 0.25, .lerpInverseEmpiricalCDF, 2.25),
26+
([1, 4, 2, 6, 4, 3], 0.5, .lerpInverseEmpiricalCDF, 3.5),
27+
([1, 4, 2, 6, 4, 3], 0.75, .lerpInverseEmpiricalCDF, 4),
28+
([1, 4, 2, 6, 4, 3], 1, .lerpInverseEmpiricalCDF, 6),
29+
([7, 3, 10, 1, 8, 4, 9, 2, 5, 6], 0.5, .inverseEmpiricalCDF, 6),
30+
([5, 9, 2, 7, 1, 8, 3, 6, 4], 0.5, .inverseEmpiricalCDF, 5),
31+
([7, 3, 10, 1, 8, 4, 9, 2, 5, 6], 0.5, .averagedInverseEmpiricalCDF, 5.5),
32+
([5, 9, 2, 7, 1, 8, 3, 6, 4], 0.5, .averagedInverseEmpiricalCDF, 5),
33+
([7, 3, 10, 1, 8, 4, 9, 2, 5, 6], 0.5, .closest, 5),
34+
([5, 9, 2, 7, 1, 8, 3, 6, 4], 0.5, .closest, 5),
35+
([7, 3, 10, 1, 8, 4, 9, 2, 5, 6], 0.5, .lerpInverseEmpiricalCDF, 5.5),
36+
([5, 9, 2, 7, 1, 8, 3, 6, 4], 0.5, .lerpInverseEmpiricalCDF, 5),
37+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.25, .inverseEmpiricalCDF, 6),
38+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.5, .inverseEmpiricalCDF, 11),
39+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.75, .inverseEmpiricalCDF, 16),
40+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.25, .averagedInverseEmpiricalCDF, 5.5),
41+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.5, .averagedInverseEmpiricalCDF, 10.5),
42+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.75, .averagedInverseEmpiricalCDF, 15.5),
43+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.25, .closest, 6),
44+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.5, .closest, 11),
45+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.75, .closest, 15),
46+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.25, .lerpInverseEmpiricalCDF, 5.75),
47+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.5, .lerpInverseEmpiricalCDF, 10.5),
48+
([15, 3, 8, 20, 1, 12, 7, 18, 5, 10, 14, 2, 9, 16, 4, 11, 19, 6, 13, 17], 0.75, .lerpInverseEmpiricalCDF, 15.25),
49+
([1, 2], 0, .inverseEmpiricalCDF, 1),
50+
([1, 2], 1, .inverseEmpiricalCDF, 2),
51+
([1, 2], 0.25, .inverseEmpiricalCDF, 2),
52+
([1, 2], 0.75, .inverseEmpiricalCDF, 2),
53+
([1, 2], 0, .averagedInverseEmpiricalCDF, 1),
54+
([1, 2], 1, .averagedInverseEmpiricalCDF, 2),
55+
([1, 2], 0.25, .averagedInverseEmpiricalCDF, 1.5),
56+
([1, 2], 0.75, .averagedInverseEmpiricalCDF, 1.5),
57+
([1, 2], 0, .closest, 1),
58+
([1, 2], 1, .closest, 2),
59+
([1, 2], 0.25, .closest, 1),
60+
([1, 2], 0.75, .closest, 2),
61+
([1, 2], 0, .lerpInverseEmpiricalCDF, 1),
62+
([1, 2], 1, .lerpInverseEmpiricalCDF, 2),
63+
([1, 2], 0.25, .lerpInverseEmpiricalCDF, 1.25),
64+
([1, 2], 0.75, .lerpInverseEmpiricalCDF, 1.75),
65+
([5], 0.5, .inverseEmpiricalCDF, 5),
66+
([5], 0.5, .averagedInverseEmpiricalCDF, 5),
67+
([5], 0.5, .closest, 5),
68+
([5], 0.5, .lerpInverseEmpiricalCDF, 5),
2569
] as [([Int], Double, QuantileEstimationMethod, Double)]
2670
)
2771
func validData(data: [Int], probability: Double, method: QuantileEstimationMethod, expectedQuantile: Double) async {
@@ -33,11 +77,59 @@ struct QuantileTests {
3377
arguments: [
3478
([], .inverseEmpiricalCDF),
3579
([], .averagedInverseEmpiricalCDF),
36-
([], .closestOrOddIndexed),
80+
([], .closest),
3781
([], .lerpInverseEmpiricalCDF),
3882
] as [([Int], QuantileEstimationMethod)]
3983
)
4084
func invalidData(data: [Int], method: QuantileEstimationMethod) async {
4185
#expect(data.quantile(probability: 1, of: \.self).isNaN)
4286
}
87+
88+
@Test(
89+
"Collection with NaN values returns NaN",
90+
arguments: [
91+
([.nan, 1.0, 2.0], 0.5, .inverseEmpiricalCDF),
92+
([1.0, .nan, 2.0], 0.5, .inverseEmpiricalCDF),
93+
([1.0, 2.0, .nan], 0.5, .inverseEmpiricalCDF),
94+
([.nan, 1.0, 2.0], 0.5, .averagedInverseEmpiricalCDF),
95+
([.nan, 1.0, 2.0], 0.5, .closest),
96+
([.nan, 1.0, 2.0], 0.5, .lerpInverseEmpiricalCDF),
97+
([.nan, 1.0, 2.0], 0.0, .inverseEmpiricalCDF),
98+
([.nan, 1.0, 2.0], 1.0, .inverseEmpiricalCDF),
99+
] as [([Double], Double, QuantileEstimationMethod)]
100+
)
101+
func nanInCollection(data: [Double], probability: Double, method: QuantileEstimationMethod) async {
102+
#expect(data.quantile(probability: probability, of: \.self, method: method).isNaN)
103+
}
104+
105+
@Test(
106+
"Infinity values without conflicting arithmetic return the correct quantile",
107+
arguments: [
108+
([1.0, 2.0, .infinity], 0.75, .inverseEmpiricalCDF, Double.infinity),
109+
([1.0, .infinity], 0.5, .inverseEmpiricalCDF, Double.infinity),
110+
([-.infinity, 1.0, 2.0], 0.0, .inverseEmpiricalCDF, -.infinity),
111+
([1.0, 2.0, .infinity], 1.0, .inverseEmpiricalCDF, Double.infinity),
112+
([1.0, .infinity], 0.75, .closest, Double.infinity),
113+
([-.infinity, 1.0], 0.25, .closest, -.infinity),
114+
([1.0, .infinity], 0.5, .averagedInverseEmpiricalCDF, Double.infinity),
115+
([-.infinity, 1.0], 0.5, .averagedInverseEmpiricalCDF, -.infinity),
116+
([1.0, .infinity], 0.5, .lerpInverseEmpiricalCDF, Double.infinity),
117+
] as [([Double], Double, QuantileEstimationMethod, Double)]
118+
)
119+
func infinityWellBehaved(data: [Double], probability: Double, method: QuantileEstimationMethod, expectedQuantile: Double) async {
120+
let result = data.quantile(probability: probability, of: \.self, method: method)
121+
#expect(result == expectedQuantile)
122+
}
123+
124+
@Test(
125+
"Conflicting infinities under arithmetic return NaN",
126+
arguments: [
127+
([-.infinity, .infinity], 0.5, .averagedInverseEmpiricalCDF),
128+
([-.infinity, 1.0], 0.5, .lerpInverseEmpiricalCDF),
129+
([-.infinity, .infinity], 0.5, .lerpInverseEmpiricalCDF),
130+
] as [([Double], Double, QuantileEstimationMethod)]
131+
)
132+
func conflictingInfinitiesReturnNaN(data: [Double], probability: Double, method: QuantileEstimationMethod) async {
133+
#expect(data.quantile(probability: probability, of: \.self, method: method).isNaN)
134+
}
43135
}

0 commit comments

Comments
 (0)