Skip to content

Commit

Permalink
Fixes numerical accuracy issues in quantile.
Browse files Browse the repository at this point in the history
Fixes issue JuliaStats/StatsBase.jl#164, and another when `p < eps()`.
  • Loading branch information
simonbyrne committed May 25, 2016
1 parent 3a346a6 commit 4e81d34
Show file tree
Hide file tree
Showing 2 changed files with 14 additions and 7 deletions.
16 changes: 9 additions & 7 deletions base/statistics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -577,19 +577,21 @@ end
@inline function _quantile(v::AbstractVector, p::Real)
T = float(eltype(v))
isnan(p) && return T(NaN)
0 <= p <= 1 || throw(ArgumentError("input probability out of [0,1] range"))

lv = length(v)
index = 1 + (lv-1)*p
1 <= index <= lv || error("input probability out of [0,1] range")
f0 = (lv-1)*p # 0-based interpolated index
t0 = trunc(f0)
h = f0 - t0

indlo = floor(index)
i = trunc(Int,indlo)
i = trunc(Int,t0) + 1

if index == indlo
if h == 0
return T(v[i])
else
h = T(index - indlo)
return (1-h)*T(v[i]) + h*T(v[i+1])
a = T(v[i])
b = T(v[i+1])
return a + h*(b-a)
end
end

Expand Down
5 changes: 5 additions & 0 deletions test/statistics.jl
Original file line number Diff line number Diff line change
Expand Up @@ -324,6 +324,11 @@ end
@test quantile(0.0:100.0, 0.0:0.1:1.0, sorted=true) == collect(0.0:10.0:100.0)
@test quantile(100f0:-1f0:0.0, 0.0:0.1:1.0) == collect(0f0:10f0:100f0)

@test quantile([0,1],1e-18) == 1e-18

# StatsBase issue 164
y = [0.40003674665581906,0.4085630862624367,0.41662034698690303,0.41662034698690303,0.42189053966652057,0.42189053966652057,0.42553514344518345,0.43985732442991354]
@test issorted(quantile(y, linspace(0.01, 0.99, 17)))

# test invalid hist nbins argument (#9999)
@test_throws ArgumentError hist(Int[], -1)
Expand Down

0 comments on commit 4e81d34

Please sign in to comment.