Fixes numerical accuracy issues in quantile.

Fixes issue JuliaStats/StatsBase.jl#164, and another when `p < eps()`.

base/statistics.jl

@@ -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
 

test/statistics.jl

@@ -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)