Edge cases for npqreg

examples/mcycle.csv

@@ -0,0 +1,134 @@
+"","times","accel"
+"1",2.4,0
+"2",2.6,-1.3
+"3",3.2,-2.7
+"4",3.6,0
+"5",4,-2.7
+"6",6.2,-2.7
+"7",6.6,-2.7
+"8",6.8,-1.3
+"9",7.8,-2.7
+"10",8.2,-2.7
+"11",8.8,-1.3
+"12",8.8,-2.7
+"13",9.6,-2.7
+"14",10,-2.7
+"15",10.2,-5.4
+"16",10.6,-2.7
+"17",11,-5.4
+"18",11.4,0
+"19",13.2,-2.7
+"20",13.6,-2.7
+"21",13.8,0
+"22",14.6,-13.3
+"23",14.6,-5.4
+"24",14.6,-5.4
+"25",14.6,-9.3
+"26",14.6,-16
+"27",14.6,-22.8
+"28",14.8,-2.7
+"29",15.4,-22.8
+"30",15.4,-32.1
+"31",15.4,-53.5
+"32",15.4,-54.9
+"33",15.6,-40.2
+"34",15.6,-21.5
+"35",15.8,-21.5
+"36",15.8,-50.8
+"37",16,-42.9
+"38",16,-26.8
+"39",16.2,-21.5
+"40",16.2,-50.8
+"41",16.2,-61.7
+"42",16.4,-5.4
+"43",16.4,-80.4
+"44",16.6,-59
+"45",16.8,-71
+"46",16.8,-91.1
+"47",16.8,-77.7
+"48",17.6,-37.5
+"49",17.6,-85.6
+"50",17.6,-123.1
+"51",17.6,-101.9
+"52",17.8,-99.1
+"53",17.8,-104.4
+"54",18.6,-112.5
+"55",18.6,-50.8
+"56",19.2,-123.1
+"57",19.4,-85.6
+"58",19.4,-72.3
+"59",19.6,-127.2
+"60",20.2,-123.1
+"61",20.4,-117.9
+"62",21.2,-134
+"63",21.4,-101.9
+"64",21.8,-108.4
+"65",22,-123.1
+"66",23.2,-123.1
+"67",23.4,-128.5
+"68",24,-112.5
+"69",24.2,-95.1
+"70",24.2,-81.8
+"71",24.6,-53.5
+"72",25,-64.4
+"73",25,-57.6
+"74",25.4,-72.3
+"75",25.4,-44.3
+"76",25.6,-26.8
+"77",26,-5.4
+"78",26.2,-107.1
+"79",26.2,-21.5
+"80",26.4,-65.6
+"81",27,-16
+"82",27.2,-45.6
+"83",27.2,-24.2
+"84",27.2,9.5
+"85",27.6,4
+"86",28.2,12
+"87",28.4,-21.5
+"88",28.4,37.5
+"89",28.6,46.9
+"90",29.4,-17.4
+"91",30.2,36.2
+"92",31,75
+"93",31.2,8.1
+"94",32,54.9
+"95",32,48.2
+"96",32.8,46.9
+"97",33.4,16
+"98",33.8,45.6
+"99",34.4,1.3
+"100",34.8,75
+"101",35.2,-16
+"102",35.2,-54.9
+"103",35.4,69.6
+"104",35.6,34.8
+"105",35.6,32.1
+"106",36.2,-37.5
+"107",36.2,22.8
+"108",38,46.9
+"109",38,10.7
+"110",39.2,5.4
+"111",39.4,-1.3
+"112",40,-21.5
+"113",40.4,-13.3
+"114",41.6,30.8
+"115",41.6,-10.7
+"116",42.4,29.4
+"117",42.8,0
+"118",42.8,-10.7
+"119",43,14.7
+"120",44,-1.3
+"121",44.4,0
+"122",45,10.7
+"123",46.6,10.7
+"124",47.8,-26.8
+"125",47.8,-14.7
+"126",48.8,-13.3
+"127",50.6,0
+"128",52,10.7
+"129",53.2,-14.7
+"130",55,-2.7
+"131",55,10.7
+"132",55.4,-2.7
+"133",57.6,10.7

examples/npqreg.jl

@@ -1,5 +1,5 @@
 using QuantileRegressions, Distributions, CSV
-mycyle = CSV.read("/home/pkofod/mcycle.csv")
+mycyle = CSV.read("mcycle.csv")
 x = mycyle.times
 y = mycyle.accel
 

src/InteriorPoint.jl

@@ -90,12 +90,14 @@ for it = 1:max_it
         mul!(XtX, Xtmp', X)
         F = cholesky!(Symmetric(XtX))
         mul!(Xtqr, Xtmp', r)
-        ldiv!(dy, F, Xtqr)
+        dy .= F\Xtqr
+        #ldiv!(dy, F, Xtqr)
     else
         Q = Diagonal(sqrt.(q)) # Very efficient to do since Q diagonal
         AQtF = qr(mul!(Xtmp, Q, X)) # PE 2004
         rhs = Q*r        # "
-        ldiv!(dy, AQtF, rhs)
+        dy .= AQtF\rhs
+        #ldiv!(dy, AQtF, rhs)
     end
    # dy .= AQtF\rhs
 
@@ -136,10 +138,12 @@ for it = 1:max_it
             @. rtmp = r + dxdz - dsdw - xi
             mul!(Xtqr, Xtmp', rtmp)
 
-            ldiv!(dy, F, Xtqr)
+            dy .= F\Xtqr
+            #ldiv!(dy, F, Xtqr)
         else
             BLAS.axpy!(1.0, Q * (dxdz .- dsdw .- xi), rhs) # no gemv-wrapper gemv(Q, (dxdz - dsdw - xi), rhs,1,1,n)?
-            ldiv!(dy, AQtF, rhs)
+            dy .= AQtF\rhs
+            #ldiv!(dy, AQtF, rhs)
         end
 
         mul!(tmp, X, dy)

src/npqreg.jl

@@ -1,27 +1,38 @@
 function npqreg(y, x, tau, method=IP(); m=50, h=2, xrange=nothing)
-	T = eltype(x)
-	if xrange === nothing
-		xrange = collect(range(extrema(x)..., length=m))
-	end
-	z = copy(x)
-	Z = hcat(fill(1, length(x)), z)
-	ghat = copy(xrange)
-	dghat = copy(xrange)
-	for i in 1:length(xrange)
-	  z .= x .- xrange[i]
-	  Z[:, 2] .= z
+    T = eltype(x)
+    if xrange === nothing
+        xrange = collect(range(extrema(x)..., length=m))
+    end
+    z = copy(x)
+    Z = hcat(fill(1, length(x)), z)
+    ghat = copy(xrange)
+    dghat = copy(xrange)
+    for i in 1:length(xrange)
+        z .= x .- xrange[i]
+        Z[:, 2] .= z
 
-	  w = pdf.(Normal(), z./h)
-	  widx = findall(x-> x > eps(T), w)
-	  w = w[widx]
-	  wy = w.*y[widx]
-	  wZ = w.*Z[widx, :]
+        w = pdf.(Normal(), z./h)
+        widx = findall(x-> x > eps(T), w)
+        if length(widx) > 1
+            # more than one observation with positive
+            # weights, so we can estimate value and derivative
+            w = w[widx]
+            wy = w.*y[widx]
+            wZ = w.*Z[widx, :]
+            r = qreg_coef(wy, wZ, tau, method)
+            ghat[i] = r[1]
+            dghat[i] = r[2]
+        elseif length(widx) == 0
+            # no positive weights, xrange too narrow
+            ghat[i] = NaN
+            dghat[i] = NaN
+        else
+            # only one subject with positive weights, we
+            # can get the value but not the derivative
+            ghat[i] = y[widx[1]]
+            dghat[i] = NaN
+        end
+    end
 
-	  r = qreg_coef(wy, wZ, tau, method)
-
-	  ghat[i] = r[1]
-	  dghat[i] = r[2]
-	end
-
-	xrange, ghat, dghat
+    xrange, ghat, dghat
 end

test/npqreg.jl

@@ -0,0 +1,23 @@
+using QuantileRegressions, CSV
+@testset "npqreg" begin
+	mycyle = CSV.read("../examples/mcycle.csv")
+	x = mycyle.times
+	y = mycyle.accel
+
+	for t = range(0.01,0.99,length=99)
+		res = QuantileRegressions.npqreg(Array(y), Array(x), t; h=2)
+		@test all(!isnan, res[2])
+		@test all(!isnan, res[3])
+		xrange = sort(unique(x))
+		res = QuantileRegressions.npqreg(Array(y), Array(x), t; xrange=xrange, h=0.5)
+	end
+	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.4; xrange=[x[1]], h=2)
+	@test !isnan(res[2][1])
+	@test isnan(res[3][1])
+	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1]], h=0.0001)
+	@test !isnan(res[2][1])
+	@test isnan(res[3][1])
+	res = QuantileRegressions.npqreg([y[1]], [x[1]], 0.8; xrange=[x[1].+1], h=0.1)
+	@test isnan(res[2][1])
+	@test isnan(res[3][1])
+end

test/runtests.jl

@@ -2,4 +2,5 @@ using QuantileRegressions
 using DataFrames, CSV, LinearAlgebra
 using Test
 
-include("engel_test.jl")
+include("engel_test.jl")
+include("npqreg.jl")