diff --git a/Project.toml b/Project.toml index 842361e..ed02ef1 100644 --- a/Project.toml +++ b/Project.toml @@ -40,7 +40,8 @@ julia = "1.6" [extras] DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" +Suppressor = "fd094767-a336-5f1f-9728-57cf17d0bbfb" Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40" [targets] -test = ["DataFrames", "Test"] +test = ["DataFrames", "Suppressor", "Test"] diff --git a/src/losses.jl b/src/losses.jl index ce9b210..f2e31e3 100644 --- a/src/losses.jl +++ b/src/losses.jl @@ -29,18 +29,18 @@ function weight end function estimator_norm end "The limit at ∞ of the loss function. Used for scale estimation of bounded loss." -estimator_bound(f::LossFunction) = Inf +estimator_bound(::LossFunction) = Inf "The tuning constant of the loss function, can be optimized to get efficient or robust estimates." tuning_constant(loss::L) where {L<:LossFunction} = loss.c "Boolean if the estimator or loss function is convex" -isconvex(f::LossFunction) = false -isconvex(f::ConvexLossFunction) = true +isconvex(::LossFunction) = false +isconvex(::ConvexLossFunction) = true "Boolean if the estimator or loss function is bounded" -isbounded(f::LossFunction) = false -isbounded(f::BoundedLossFunction) = true +isbounded(::LossFunction) = false +isbounded(::BoundedLossFunction) = true "The tuning constant associated to the loss that gives an efficient M-estimator." estimator_high_breakdown_point_constant(::Type{L}) where {L<:LossFunction} = 1 @@ -205,10 +205,10 @@ struct HuberLoss <: ConvexLossFunction HuberLoss() = new(estimator_high_efficiency_constant(HuberLoss)) end -_rho(l::HuberLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : (abs(r) - 1 / 2) -_psi(l::HuberLoss, r::Real) = (abs(r) <= 1) ? r : sign(r) -_psider(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0) -_weight(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : 1 / abs(r) +_rho(::HuberLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : (abs(r) - 1 / 2) +_psi(::HuberLoss, r::Real) = (abs(r) <= 1) ? r : sign(r) +_psider(::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0) +_weight(::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : 1 / abs(r) function estimator_values(l::HuberLoss, r::Real) rr = abs(r) if rr <= l.c @@ -233,10 +233,10 @@ struct L1L2Loss <: ConvexLossFunction L1L2Loss(c::Real) = new(c) L1L2Loss() = new(estimator_high_efficiency_constant(L1L2Loss)) end -_rho(l::L1L2Loss, r::Real) = (sqrt(1 + r^2) - 1) -_psi(l::L1L2Loss, r::Real) = r / sqrt(1 + r^2) -_psider(l::L1L2Loss, r::Real) = 1 / (1 + r^2)^(3 / 2) -_weight(l::L1L2Loss, r::Real) = 1 / sqrt(1 + r^2) +_rho(::L1L2Loss, r::Real) = (sqrt(1 + r^2) - 1) +_psi(::L1L2Loss, r::Real) = r / sqrt(1 + r^2) +_psider(::L1L2Loss, r::Real) = 1 / (1 + r^2)^(3 / 2) +_weight(::L1L2Loss, r::Real) = 1 / sqrt(1 + r^2) function estimator_values(l::L1L2Loss, r::Real) sqr = sqrt(1 + (r / l.c)^2) return ((sqr - 1), r / sqr, 1 / sqr) @@ -257,10 +257,10 @@ struct FairLoss <: ConvexLossFunction FairLoss(c::Real) = new(c) FairLoss() = new(estimator_high_efficiency_constant(FairLoss)) end -_rho(l::FairLoss, r::Real) = abs(r) - log(1 + abs(r)) -_psi(l::FairLoss, r::Real) = r / (1 + abs(r)) -_psider(l::FairLoss, r::Real) = 1 / (1 + abs(r))^2 -_weight(l::FairLoss, r::Real) = 1 / (1 + abs(r)) +_rho(::FairLoss, r::Real) = abs(r) - log(1 + abs(r)) +_psi(::FairLoss, r::Real) = r / (1 + abs(r)) +_psider(::FairLoss, r::Real) = 1 / (1 + abs(r))^2 +_weight(::FairLoss, r::Real) = 1 / (1 + abs(r)) function estimator_values(l::FairLoss, r::Real) ir = 1 / (1 + abs(r / l.c)) return (abs(r) / l.c + log(ir), r * ir, ir) @@ -279,10 +279,10 @@ struct LogcoshLoss <: ConvexLossFunction LogcoshLoss(c::Real) = new(c) LogcoshLoss() = new(estimator_high_efficiency_constant(LogcoshLoss)) end -_rho(l::LogcoshLoss, r::Real) = log(cosh(r)) -_psi(l::LogcoshLoss, r::Real) = tanh(r) -_psider(l::LogcoshLoss, r::Real) = 1 / (cosh(r))^2 -_weight(l::LogcoshLoss, r::Real) = (abs(r) < DELTA) ? (1 - (r)^2 / 3) : tanh(r) / r +_rho(::LogcoshLoss, r::Real) = log(cosh(r)) +_psi(::LogcoshLoss, r::Real) = tanh(r) +_psider(::LogcoshLoss, r::Real) = 1 / (cosh(r))^2 +_weight(::LogcoshLoss, r::Real) = (abs(r) < DELTA) ? (1 - (r)^2 / 3) : tanh(r) / r function estimator_values(l::LogcoshLoss, r::Real) tr = l.c * tanh(r / l.c) rr = abs(r / l.c) @@ -302,10 +302,10 @@ struct ArctanLoss <: ConvexLossFunction ArctanLoss(c::Real) = new(c) ArctanLoss() = new(estimator_high_efficiency_constant(ArctanLoss)) end -_rho(l::ArctanLoss, r::Real) = r * atan(r) - 1 / 2 * log(1 + r^2) -_psi(l::ArctanLoss, r::Real) = atan(r) -_psider(l::ArctanLoss, r::Real) = 1 / (1 + r^2) -_weight(l::ArctanLoss, r::Real) = (abs(r) < DELTA) ? (1 - r^2 / 3) : atan(r) / r +_rho(::ArctanLoss, r::Real) = r * atan(r) - 1 / 2 * log(1 + r^2) +_psi(::ArctanLoss, r::Real) = atan(r) +_psider(::ArctanLoss, r::Real) = 1 / (1 + r^2) +_weight(::ArctanLoss, r::Real) = (abs(r) < DELTA) ? (1 - r^2 / 3) : atan(r) / r function estimator_values(l::ArctanLoss, r::Real) ar = atan(r / l.c) rr = abs(r / l.c) @@ -332,13 +332,13 @@ struct CatoniWideLoss <: ConvexLossFunction CatoniWideLoss() = new(estimator_high_efficiency_constant(CatoniWideLoss)) end -function _rho(l::CatoniWideLoss, r::Real) +function _rho(::CatoniWideLoss, r::Real) ar = abs(r) return (1 + ar) * log(1 + ar + r^2 / 2) - 2 * ar + 2 * atan(1 + r) - π / 2 end -_psi(l::CatoniWideLoss, r::Real) = sign(r) * log(1 + abs(r) + r^2 / 2) -_psider(l::CatoniWideLoss, r::Real) = (1 + abs(r)) / (1 + abs(r) + r^2 / 2) -function _weight(l::CatoniWideLoss, r::Real) +_psi(::CatoniWideLoss, r::Real) = sign(r) * log(1 + abs(r) + r^2 / 2) +_psider(::CatoniWideLoss, r::Real) = (1 + abs(r)) / (1 + abs(r) + r^2 / 2) +function _weight(::CatoniWideLoss, r::Real) return (abs(r) < DELTA) ? oftype(r, 1) : log(1 + abs(r) + r^2 / 2) / abs(r) end function estimator_values(l::CatoniWideLoss, r::Real) @@ -368,26 +368,26 @@ struct CatoniNarrowLoss <: ConvexLossFunction CatoniNarrowLoss() = new(estimator_high_efficiency_constant(CatoniNarrowLoss)) end -function _rho(l::CatoniNarrowLoss, r::Real) +function _rho(::CatoniNarrowLoss, r::Real) ar = abs(r) if ar >= 1 return (ar - 1) * log(2) + 2 - π / 2 end return (1 - ar) * log(1 - ar + r^2 / 2) + 2 * ar + 2 * atan(1 - ar) - π / 2 end -function _psi(l::CatoniNarrowLoss, r::Real) +function _psi(::CatoniNarrowLoss, r::Real) if abs(r) <= 1 return -sign(r) * log(1 - abs(r) + r^2 / 2) end return sign(r) * log(2) end -function _psider(l::CatoniNarrowLoss, r::Real) +function _psider(::CatoniNarrowLoss, r::Real) if abs(r) <= 1 return (1 - abs(r)) / (1 - abs(r) + r^2 / 2) end return 0 end -function _weight(l::CatoniNarrowLoss, r::Real) +function _weight(::CatoniNarrowLoss, r::Real) if r == 0 return oftype(r, 1) elseif abs(r) <= 1 @@ -422,18 +422,16 @@ struct CauchyLoss <: LossFunction CauchyLoss(c::Real) = new(c) CauchyLoss() = new(estimator_high_efficiency_constant(CauchyLoss)) end -_rho(l::CauchyLoss, r::Real) = log(1 + r^2) # * 1/2 # remove factor 1/2 so the loss has a norm -_psi(l::CauchyLoss, r::Real) = r / (1 + r^2) -_psider(l::CauchyLoss, r::Real) = (1 - r^2) / (1 + r^2)^2 -_weight(l::CauchyLoss, r::Real) = 1 / (1 + r^2) +_rho(::CauchyLoss, r::Real) = log(1 + r^2) # * 1/2 # remove factor 1/2 so the loss has a norm +_psi(::CauchyLoss, r::Real) = r / (1 + r^2) +_psider(::CauchyLoss, r::Real) = (1 - r^2) / (1 + r^2)^2 +_weight(::CauchyLoss, r::Real) = 1 / (1 + r^2) function estimator_values(l::CauchyLoss, r::Real) ir = 1 / (1 + (r / l.c)^2) return (-log(ir), r * ir, ir) end estimator_norm(l::CauchyLoss) = l.c * π estimator_bound(l::CauchyLoss) = 1.0 -isconvex(::CauchyLoss) = false -isbounded(::CauchyLoss) = false estimator_high_efficiency_constant(::Type{CauchyLoss}) = 2.385 estimator_high_breakdown_point_constant(::Type{CauchyLoss}) = 1.468 @@ -449,10 +447,10 @@ struct GemanLoss <: BoundedLossFunction GemanLoss(c::Real) = new(c) GemanLoss() = new(estimator_high_efficiency_constant(GemanLoss)) end -_rho(l::GemanLoss, r::Real) = 1 / 2 * r^2 / (1 + r^2) -_psi(l::GemanLoss, r::Real) = r / (1 + r^2)^2 -_psider(l::GemanLoss, r::Real) = (1 - 3 * r^2) / (1 + r^2)^3 -_weight(l::GemanLoss, r::Real) = 1 / (1 + r^2)^2 +_rho(::GemanLoss, r::Real) = 1 / 2 * r^2 / (1 + r^2) +_psi(::GemanLoss, r::Real) = r / (1 + r^2)^2 +_psider(::GemanLoss, r::Real) = (1 - 3 * r^2) / (1 + r^2)^3 +_weight(::GemanLoss, r::Real) = 1 / (1 + r^2)^2 function estimator_values(l::GemanLoss, r::Real) rr2 = (r / l.c)^2 ir = 1 / (1 + rr2) @@ -460,8 +458,6 @@ function estimator_values(l::GemanLoss, r::Real) end estimator_norm(::GemanLoss) = Inf estimator_bound(::GemanLoss) = 1 / 2 -isconvex(::GemanLoss) = false -isbounded(::GemanLoss) = true estimator_high_efficiency_constant(::Type{GemanLoss}) = 3.787 estimator_high_breakdown_point_constant(::Type{GemanLoss}) = 0.61200 @@ -478,10 +474,10 @@ struct WelschLoss <: BoundedLossFunction WelschLoss(c::Real) = new(c) WelschLoss() = new(estimator_high_efficiency_constant(WelschLoss)) end -_rho(l::WelschLoss, r::Real) = -1 / 2 * Base.expm1(-r^2) -_psi(l::WelschLoss, r::Real) = r * exp(-r^2) -_psider(l::WelschLoss, r::Real) = (1 - 2 * r^2) * exp(-r^2) -_weight(l::WelschLoss, r::Real) = exp(-r^2) +_rho(::WelschLoss, r::Real) = -1 / 2 * Base.expm1(-r^2) +_psi(::WelschLoss, r::Real) = r * exp(-r^2) +_psider(::WelschLoss, r::Real) = (1 - 2 * r^2) * exp(-r^2) +_weight(::WelschLoss, r::Real) = exp(-r^2) function estimator_values(l::WelschLoss, r::Real) rr2 = (r / l.c)^2 er = exp(-rr2) @@ -489,8 +485,6 @@ function estimator_values(l::WelschLoss, r::Real) end estimator_norm(::WelschLoss) = Inf estimator_bound(::WelschLoss) = 1 / 2 -isconvex(::WelschLoss) = false -isbounded(::WelschLoss) = true estimator_high_efficiency_constant(::Type{WelschLoss}) = 2.985 estimator_high_breakdown_point_constant(::Type{WelschLoss}) = 0.8165 @@ -507,18 +501,16 @@ struct TukeyLoss <: BoundedLossFunction TukeyLoss(c::Real) = new(c) TukeyLoss() = new(estimator_high_efficiency_constant(TukeyLoss)) end -_rho(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 / 6 * (1 - (1 - r^2)^3) : 1 / 6 -_psi(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? r * (1 - r^2)^2 : oftype(r, 0) -_psider(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 - 6 * r^2 + 5 * r^4 : oftype(r, 0) -_weight(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? (1 - r^2)^2 : oftype(r, 0) +_rho(::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 / 6 * (1 - (1 - r^2)^3) : 1 / 6 +_psi(::TukeyLoss, r::Real) = (abs(r) <= 1) ? r * (1 - r^2)^2 : oftype(r, 0) +_psider(::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 - 6 * r^2 + 5 * r^4 : oftype(r, 0) +_weight(::TukeyLoss, r::Real) = (abs(r) <= 1) ? (1 - r^2)^2 : oftype(r, 0) function estimator_values(l::TukeyLoss, r::Real) pr = (abs(r) <= l.c) * (1 - (r / l.c)^2) return (1 / 6 * (1 - pr^3), r * pr^2, pr^2) end estimator_norm(::TukeyLoss) = Inf estimator_bound(::TukeyLoss) = 1 / 6 -isconvex(::TukeyLoss) = false -isbounded(::TukeyLoss) = true estimator_high_efficiency_constant(::Type{TukeyLoss}) = 4.685 estimator_high_breakdown_point_constant(::Type{TukeyLoss}) = 1.5476 @@ -590,8 +582,6 @@ function estimator_values(l::YohaiZamarLoss, r::Real) end estimator_norm(::YohaiZamarLoss) = Inf estimator_bound(::YohaiZamarLoss) = 1 -isconvex(::YohaiZamarLoss) = false -isbounded(::YohaiZamarLoss) = true estimator_high_efficiency_constant(::Type{YohaiZamarLoss}) = 3.1806 estimator_high_breakdown_point_constant(::Type{YohaiZamarLoss}) = 1.2139 @@ -607,9 +597,9 @@ struct HardThresholdLoss <: BoundedLossFunction HardThresholdLoss(c::Real) = new(c) HardThresholdLoss() = new(estimator_high_efficiency_constant(HardThresholdLoss)) end -_rho(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : oftype(r, 1 / 2) -_psi(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r : oftype(r, 0) -function _psider(l::HardThresholdLoss, r::Real) +_rho(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : oftype(r, 1 / 2) +_psi(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r : oftype(r, 0) +function _psider(::HardThresholdLoss, r::Real) if abs(r) <= 1 return oftype(r, 1) elseif abs(r) < 1 + DELTA @@ -618,7 +608,7 @@ function _psider(l::HardThresholdLoss, r::Real) return oftype(r, 0) end end -_weight(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0) +_weight(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0) function estimator_values(l::HardThresholdLoss, r::Real) ar = abs(r / l.c) if ar <= 1 @@ -629,8 +619,6 @@ function estimator_values(l::HardThresholdLoss, r::Real) end estimator_norm(::HardThresholdLoss) = Inf estimator_bound(::HardThresholdLoss) = 1 / 2 -isconvex(::HardThresholdLoss) = false -isbounded(::HardThresholdLoss) = true ## Computed with `find_zero(c -> I1(c) - 0.95, 2.8, Order1())` after simplification (fun_eff = I1) estimator_high_efficiency_constant(::Type{HardThresholdLoss}) = 2.795 @@ -709,8 +697,6 @@ end estimator_norm(::HampelLoss) = Inf estimator_bound(l::HampelLoss) = (l.ν1 + l.ν2 - 1) / 2 -isconvex(::HampelLoss) = false -isbounded(::HampelLoss) = true # Values of `c` for (ν1, ν2) = (2, 4) estimator_high_efficiency_constant(::Type{HampelLoss}) = 1.382 diff --git a/test/estimators.jl b/test/estimators.jl index 3cff609..594e6e5 100644 --- a/test/estimators.jl +++ b/test/estimators.jl @@ -19,129 +19,154 @@ using RobustModels: emp_norm(l::LossFunction) = 2 * quadgk(x -> exp(-RobustModels.rho(l, x)), 0, Inf)[1] +@testset "Losses and estimators" begin -@testset "Methods loss functions: $(name)" for name in losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - l = typeloss() - - @testset "Methods estimators: $(estimator)" for estimator in ( - nothing, "M", "S", "MM", "Tau", "GeneralizedQuantile" - ) - # Check LossFunction methods - if isnothing(estimator) - estimator_name = "Loss function" - typest = typeloss - - # Check AbstractEstimator methods - else - estimator_name = "$(estimator) Estimator" - T = getproperty(RobustModels, Symbol(estimator * "Estimator")) - if estimator in ("S", "MM", "Tau") - if !in(name, bounded_losses) - @test_throws TypeError T{typeloss} - continue + @testset "loss functions: $(name)" for name in losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + l = typeloss() + + @testset "Methods estimators: $(estimator)" for estimator in ( + nothing, "M", "S", "MM", "Tau", "GeneralizedQuantile" + ) + # Check LossFunction methods + if isnothing(estimator) + estimator_name = "Loss function" + typest = typeloss + + # Check AbstractEstimator methods + else + estimator_name = "$(estimator) Estimator" + T = getproperty(RobustModels, Symbol(estimator * "Estimator")) + if estimator in ("S", "MM", "Tau") + if !in(name, bounded_losses) + @test_throws TypeError T{typeloss} + continue + end end + typest = T{typeloss} end - typest = T{typeloss} - end - est = typest() - @test_nowarn println(est) - - if !isnothing(estimator) - if estimator == "Tau" - # @test isa(loss(est), Tuple{BoundedLossFunction, BoundedLossFunction}) - @test isa(loss(est), CompositeLossFunction) - @test typeof(first(loss(est))) == typeloss - @test typeof(last(loss(est))) == typeloss - else - @test typeof(loss(est)) == typeloss + est = typest() + @test_nowarn show(devnull, est) + + if !isnothing(estimator) + if estimator == "Tau" + # @test isa(loss(est), Tuple{BoundedLossFunction, BoundedLossFunction}) + @test isa(loss(est), CompositeLossFunction) + @test typeof(first(loss(est))) == typeloss + @test typeof(last(loss(est))) == typeloss + else + @test typeof(loss(est)) == typeloss + end end - end - @testset "Bounded $(estimator_name): $(name)" begin - if name in bounded_losses - @test isbounded(est) - else - @test !isbounded(est) + @testset "Bounded $(estimator_name): $(name)" begin + if name in bounded_losses + @test isbounded(est) + else + @test !isbounded(est) + end end - end - @testset "Convex $(estimator_name): $(name)" begin - if name in convex_losses - @test isconvex(est) - else - @test !isconvex(est) + @testset "Convex $(estimator_name): $(name)" begin + if name in convex_losses + @test isconvex(est) + else + @test !isconvex(est) + end end - end - - @testset "$(estimator_name) values: $(name)" begin - ρ = rho(est, 1) - ψ = psi(est, 1) - ψp = psider(est, 1) - w = weight(est, 1) - - vals = estimator_values(est, 1) - @test length(vals) == 3 - @test vals[1] ≈ ρ rtol = 1e-6 - @test vals[2] ≈ ψ rtol = 1e-6 - @test vals[3] ≈ w rtol = 1e-6 - end - # Only for LossFunction - if isnothing(estimator) - if !isbounded(est) - @testset "Estimator norm: $(name)" begin - @test emp_norm(est) ≈ RobustModels.estimator_norm(est) rtol = 1e-5 - end + @testset "$(estimator_name) values: $(name)" begin + ρ = rho(est, 1) + ψ = psi(est, 1) + ψp = psider(est, 1) + w = weight(est, 1) + + vals = estimator_values(est, 1) + @test length(vals) == 3 + @test vals[1] ≈ ρ rtol = 1e-6 + @test vals[2] ≈ ψ rtol = 1e-6 + @test vals[3] ≈ w rtol = 1e-6 end - if !in(name, ("L2", "L1")) - @testset "Estimator high efficiency: $(name)" begin - vopt = estimator_high_efficiency_constant(typest) - if name != "HardThreshold" - v = efficiency_tuning_constant(typest; eff=0.95, c0=0.9 * vopt) - @test isapprox(v, vopt; rtol=1e-3) + # Only for LossFunction + if isnothing(estimator) + @testset "Loss methods: $(name)" begin + @test loss(est) == est + + # estimator_bound + if !isconvex(est) + @test isfinite(RobustModels.estimator_bound(est)) + else + @test !isfinite(RobustModels.estimator_bound(est)) end - end - end - if isbounded(est) - @testset "Estimator high breakdown point: $(name)" begin - vopt = estimator_high_breakdown_point_constant(typest) - v = breakdown_point_tuning_constant(typest; bp=0.5, c0=1.1 * vopt) - @test isapprox(v, vopt; rtol=1e-3) - end + # tuning_constant + @test isfinite(RobustModels.tuning_constant(est)) + + @testset "Estimator norm: $(name)" begin + est_norm = RobustModels.estimator_norm(est) + if !isbounded(est) + @test emp_norm(est) ≈ est_norm rtol = 1e-5 + else + @test !isfinite(est_norm) + end + end - @testset "τ-Estimator high efficiency: $(name)" begin - vopt = estimator_tau_efficient_constant(typest) - if name != "HardThreshold" - v = tau_efficiency_tuning_constant(typest; eff=0.95, c0=1.1 * vopt) - @test isapprox(v, vopt; rtol=1e-3) + if !in(name, ("L2", "L1")) + @testset "Estimator high efficiency: $(name)" begin + vopt = estimator_high_efficiency_constant(typest) + if name != "HardThreshold" + v = efficiency_tuning_constant( + typest; eff=0.95, c0=0.9 * vopt + ) + @test isapprox(v, vopt; rtol=1e-3) + end + end + end + + if isbounded(est) + @testset "Estimator high breakdown point: $(name)" begin + vopt = estimator_high_breakdown_point_constant(typest) + v = breakdown_point_tuning_constant( + typest; bp=0.5, c0=1.1 * vopt + ) + @test isapprox(v, vopt; rtol=1e-3) + end + + @testset "τ-Estimator high efficiency: $(name)" begin + vopt = estimator_tau_efficient_constant(typest) + if name != "HardThreshold" + v = tau_efficiency_tuning_constant( + typest; eff=0.95, c0=1.1 * vopt + ) + @test isapprox(v, vopt; rtol=1e-3) + end + end end end end - end - if typest <: AbstractQuantileEstimator - @testset "MQuantile: $(name)" begin - τ = 0.5 - qest1 = GeneralizedQuantileEstimator(l, τ) - qest2 = GeneralizedQuantileEstimator{typeloss}(τ) - @test qest1 == qest2 - if name == "L2" - qest4 = ExpectileEstimator(τ) - @test qest1 == qest4 - elseif name == "L1" - qest4 = RobustModels.QuantileEstimator(τ) - @test qest1 == qest4 - end + if typest <: AbstractQuantileEstimator + @testset "MQuantile: $(name)" begin + τ = 0.5 + qest1 = GeneralizedQuantileEstimator(l, τ) + qest2 = GeneralizedQuantileEstimator{typeloss}(τ) + @test qest1 == qest2 + if name == "L2" + qest4 = ExpectileEstimator(τ) + @test qest1 == qest4 + elseif name == "L1" + qest4 = RobustModels.QuantileEstimator(τ) + @test qest1 == qest4 + end - @testset "Pass through of methods" for fun in (rho, psi, psider, weight) - ρ1 = fun(est, 1) - ρ2 = fun(qest1, 1) - @test ρ1 == ρ2 - end + @testset "Pass through of methods" for fun in (rho, psi, psider, weight) + ρ1 = fun(est, 1) + ρ2 = fun(qest1, 1) + @test ρ1 == ρ2 + end + end end end end diff --git a/test/interface.jl b/test/interface.jl index ac7be31..83790b4 100644 --- a/test/interface.jl +++ b/test/interface.jl @@ -1,6 +1,7 @@ using Tables using Missings +using Suppressor: @capture_err m1 = fit(LinearModel, form, data) λlm = dispersion(m1) @@ -11,15 +12,17 @@ est1 = MEstimator(loss1) est2 = MEstimator(loss2) -@testset "linear: L2 estimator" begin +@testset "Interface" begin VERBOSE && println("\n\t\u25CF Estimator: L2") # OLS VERBOSE && println(m1) VERBOSE && println(" lm : ", coef(m1)) + #! format: off # Formula, dense and sparse entry and methods :cg and :chol - @testset "$(typeof(A)),\t$(method)" for (A, b) in data_tuples, method in nopen_methods + @testset "interface: $(typeof(A)),\t$(method)" for (A, b) in data_tuples, method in nopen_methods + #! format: on name = if (A == form) "formula" @@ -152,7 +155,13 @@ est2 = MEstimator(loss2) else b_mod = nothing end - @test_throws ArgumentError fit(RobustLinearModel, A, b_mod, est1) + output = @capture_err begin + @test_throws ArgumentError fit(RobustLinearModel, A, b_mod, est1) + end + if VERBOSE && length(output) > 0 + println(output) + end + if typemod == "missing" @test_nowarn fit(RobustLinearModel, A, b_mod, est1; dropmissing=true) end @@ -170,11 +179,17 @@ est2 = MEstimator(loss2) RobustLinearModel, A_mod, b_mod, est1; dropmissing=true ) else - @test_throws MethodError fit(RobustLinearModel, A_mod, b, est1) - @test_throws MethodError fit(RobustLinearModel, A, b_mod, est1) - @test_throws MethodError fit(RobustLinearModel, A_mod, b_mod, est1) + output = @capture_err begin + @test_throws MethodError fit(RobustLinearModel, A_mod, b, est1) + @test_throws MethodError fit(RobustLinearModel, A, b_mod, est1) + @test_throws MethodError fit(RobustLinearModel, A_mod, b_mod, est1) + end + if VERBOSE && length(output) > 0 + println(output) + end end end end end + end diff --git a/test/linearfit.jl b/test/linearfit.jl index 241abe7..c0c06d6 100644 --- a/test/linearfit.jl +++ b/test/linearfit.jl @@ -8,254 +8,259 @@ loss2 = RobustModels.TukeyLoss() est1 = MEstimator(loss1) est2 = MEstimator(loss2) -@testset "linear: M-estimator $(name)" for name in losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - l = typeloss() - est = MEstimator(l) - est3 = MEstimator{typeloss}() - kwargs = (; initial_scale=:mad) - if name == "L1" - kwargs = (; initial_scale=:mad, maxiter=100) - end - m2 = fit(RobustLinearModel, form, data, est; verbose=false, kwargs...) - @test all(isfinite.(coef(m2))) +@testset "Estimators fit" begin - if name != "L1" - @test_nowarn show(devnull, m2) - # else - # @test_warn L1_warning show(devnull, m2) - end + @testset "linear: M-estimator $(name)" for name in losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + l = typeloss() + est = MEstimator(l) + est3 = MEstimator{typeloss}() + kwargs = (; initial_scale=:mad) + if name == "L1" + kwargs = (; initial_scale=:mad, maxiter=100) + end - VERBOSE && println("\n\t\u25CF Estimator: $(name)") - VERBOSE && println(" lm : ", coef(m1)) - VERBOSE && println("rlm(:auto) : ", coef(m2)) + m2 = fit(RobustLinearModel, form, data, est; verbose=false, kwargs...) + @test all(isfinite.(coef(m2))) - @testset "method: $(method)" for method in nopen_methods - m3 = fit(RobustLinearModel, form, data, est; method=method, kwargs...) - m4 = rlm(form, data, est; method=method, kwargs...) - m5 = rlm(form, data, est3; method=method, kwargs...) + if name != "L1" + @test_nowarn show(devnull, m2) + # else + # @test_warn L1_warning show(devnull, m2) + end - VERBOSE && println("rlm($(method)) : ", coef(m3)) - VERBOSE && println("rlm2($(method)): ", coef(m4)) - VERBOSE && println("rlm3($(method)): ", coef(m5)) + VERBOSE && println("\n\t\u25CF Estimator: $(name)") + VERBOSE && println(" lm : ", coef(m1)) + VERBOSE && println("rlm(:auto) : ", coef(m2)) - if name != "L1" - @test isapprox(coef(m2), coef(m3); rtol=1e-5) + @testset "method: $(method)" for method in nopen_methods + m3 = fit(RobustLinearModel, form, data, est; method=method, kwargs...) + m4 = rlm(form, data, est; method=method, kwargs...) + m5 = rlm(form, data, est3; method=method, kwargs...) + + VERBOSE && println("rlm($(method)) : ", coef(m3)) + VERBOSE && println("rlm2($(method)): ", coef(m4)) + VERBOSE && println("rlm3($(method)): ", coef(m5)) + + if name != "L1" + @test isapprox(coef(m2), coef(m3); rtol=1e-5) + end end + + # refit + β2 = copy(coef(m2)) + refit!(m2, y; wts=weights(m2), verbose=false, kwargs...) + @test all(coef(m2) .== β2) + + # empty refit + refit!(m2; kwargs...) + @test all(coef(m2) .== β2) end - # refit - β2 = copy(coef(m2)) - refit!(m2, y; wts=weights(m2), verbose=false, kwargs...) - @test all(coef(m2) .== β2) + not_bounded_losses = setdiff(Set(losses), Set(bounded_losses)) - # empty refit - refit!(m2; kwargs...) - @test all(coef(m2) .== β2) -end + @testset "linear: S-estimator" begin + VERBOSE && println("\n\t\u25CF Estimator type: S-estimators") + VERBOSE && println(" lm : ", coef(m1)) + @testset "Not bounded: $(name)" for name in not_bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + @test_throws TypeError SEstimator{typeloss}() + end -not_bounded_losses = setdiff(Set(losses), Set(bounded_losses)) + @testset "Bounded: $(name)" for name in bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + est = SEstimator{typeloss}() + m = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + ) + VERBOSE && println("rlm($name) : ", coef(m)) + ## TODO: find better test + @test size(coef(m), 1) == 2 -@testset "linear: S-estimator" begin - VERBOSE && println("\n\t\u25CF Estimator type: S-estimators") - VERBOSE && println(" lm : ", coef(m1)) - @testset "Not bounded: $(name)" for name in not_bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - @test_throws TypeError SEstimator{typeloss}() + # Resampling + rng = MersenneTwister(1) + opts = Dict(:Nsteps_β => 3, :rng => rng) + m2 = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + resample=true, + resampling_options=opts, + ) + @test (scale(m2) / scale(m) - 1) <= 1e-2 + end + + @testset "Infinite loop: issue #32" begin + # Infinite loop: issue #32 + Xs = reshape( + [ + 0.001481, + 0.0017, + 0.00133, + 0.001853, + 0.002086, + 0.003189, + 0.001161, + 0.002441, + 0.001133, + 0.001308, + 0.001, + 0.009309, + 0.1456, + 0.3127, + 0.2627, + 0.1704, + 0.101, + 0.06855, + 0.02578, + ], + :, + 1, + ) + ys = [ + 1.222, + 1.599, + 2.238, + 2.233, + 2.668, + 2.637, + 3.177, + 2.539, + 2.339, + 1.481, + 1.733, + 0.04986, + 0.0812, + 0.1057, + 0.1197, + 0.1348, + 0.1006, + 0.1021, + 0.08278, + ] + @test_throws Exception rlm(Xs, ys, SEstimator{TukeyLoss}(), initial_scale=:mad) + end end - @testset "Bounded: $(name)" for name in bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - est = SEstimator{typeloss}() - m = fit( - RobustLinearModel, - form, - data, - est; - method=:cg, - verbose=false, - initial_scale=:mad, - ) - VERBOSE && println("rlm($name) : ", coef(m)) - ## TODO: find better test - @test size(coef(m), 1) == 2 + @testset "linear: MM-estimator" begin + VERBOSE && println("\n\t\u25CF Estimator type: MM-estimators") + VERBOSE && println(" lm : ", coef(m1)) + @testset "Not bounded: $(name)" for name in not_bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + @test_throws TypeError MMEstimator{typeloss}() + end - # Resampling - rng = MersenneTwister(1) - opts = Dict(:Nsteps_β => 3, :rng => rng) - m2 = fit( - RobustLinearModel, - form, - data, - est; - method=:cg, - verbose=false, - initial_scale=:mad, - resample=true, - resampling_options=opts, - ) - @test (scale(m2) / scale(m) - 1) <= 1e-2 - end + @testset "Bounded: $(name)" for name in bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + est = MMEstimator{typeloss}() + m = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + ) + VERBOSE && println("rlm($name) : ", coef(m)) + ## TODO: find better test + @test size(coef(m), 1) == 2 - @testset "Infinite loop: issue #32" begin - # Infinite loop: issue #32 - Xs = reshape( - [ - 0.001481, - 0.0017, - 0.00133, - 0.001853, - 0.002086, - 0.003189, - 0.001161, - 0.002441, - 0.001133, - 0.001308, - 0.001, - 0.009309, - 0.1456, - 0.3127, - 0.2627, - 0.1704, - 0.101, - 0.06855, - 0.02578, - ], - :, - 1, - ) - ys = [ - 1.222, - 1.599, - 2.238, - 2.233, - 2.668, - 2.637, - 3.177, - 2.539, - 2.339, - 1.481, - 1.733, - 0.04986, - 0.0812, - 0.1057, - 0.1197, - 0.1348, - 0.1006, - 0.1021, - 0.08278, - ] - @test_throws Exception rlm(Xs, ys, SEstimator{TukeyLoss}(), initial_scale=:mad) + # Resampling + rng = MersenneTwister(1) + opts = Dict(:Npoints => 10, :Nsteps_β => 3, :Nsteps_σ => 3, :rng => rng) + m2 = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + resample=true, + resampling_options=opts, + ) + @test (scale(m2) / scale(m) - 1) <= 1e-2 + end end -end -@testset "linear: MM-estimator" begin - VERBOSE && println("\n\t\u25CF Estimator type: MM-estimators") - VERBOSE && println(" lm : ", coef(m1)) - @testset "Not bounded: $(name)" for name in not_bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - @test_throws TypeError MMEstimator{typeloss}() - end + @testset "linear: τ-estimator" begin + VERBOSE && println("\n\t\u25CF Estimator type: τ-estimators") + VERBOSE && println(" lm : ", coef(m1)) + @testset "Not bounded: $(name)" for name in not_bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + @test_throws TypeError TauEstimator{typeloss}() + end - @testset "Bounded: $(name)" for name in bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - est = MMEstimator{typeloss}() - m = fit( - RobustLinearModel, - form, - data, - est; - method=:cg, - verbose=false, - initial_scale=:mad, - ) - VERBOSE && println("rlm($name) : ", coef(m)) - ## TODO: find better test - @test size(coef(m), 1) == 2 + @testset "Bounded: $(name)" for name in bounded_losses + typeloss = getproperty(RobustModels, Symbol(name * "Loss")) + est = TauEstimator{typeloss}() - # Resampling - rng = MersenneTwister(1) - opts = Dict(:Npoints => 10, :Nsteps_β => 3, :Nsteps_σ => 3, :rng => rng) - m2 = fit( - RobustLinearModel, - form, - data, - est; - method=:cg, - verbose=false, - initial_scale=:mad, - resample=true, - resampling_options=opts, - ) - @test (scale(m2) / scale(m) - 1) <= 1e-2 - end -end + m = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + ) + VERBOSE && println("rlm($name) : ", coef(m)) + ## TODO: find better test + @test size(coef(m), 1) == 2 -@testset "linear: τ-estimator" begin - VERBOSE && println("\n\t\u25CF Estimator type: τ-estimators") - VERBOSE && println(" lm : ", coef(m1)) - @testset "Not bounded: $(name)" for name in not_bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - @test_throws TypeError TauEstimator{typeloss}() + # Resampling + rng = MersenneTwister(1) + opts = Dict(:Nsteps_β => 3, :Nsteps_σ => 3, :rng => rng) + m2 = fit( + RobustLinearModel, + form, + data, + est; + method=:cg, + verbose=false, + initial_scale=:mad, + resample=true, + resampling_options=opts, + ) + @test (tauscale(m2) / tauscale(m) - 1) <= 1e-2 + end end - @testset "Bounded: $(name)" for name in bounded_losses - typeloss = getproperty(RobustModels, Symbol(name * "Loss")) - est = TauEstimator{typeloss}() - - m = fit( + @testset "linear: leverage correction" begin + m2 = fit( RobustLinearModel, form, data, - est; + est2; method=:cg, - verbose=false, initial_scale=:mad, + correct_leverage=false, ) - VERBOSE && println("rlm($name) : ", coef(m)) - ## TODO: find better test - @test size(coef(m), 1) == 2 - - # Resampling - rng = MersenneTwister(1) - opts = Dict(:Nsteps_β => 3, :Nsteps_σ => 3, :rng => rng) - m2 = fit( + m3 = fit( RobustLinearModel, form, data, - est; + est2; method=:cg, - verbose=false, initial_scale=:mad, - resample=true, - resampling_options=opts, + correct_leverage=true, ) - @test (tauscale(m2) / tauscale(m) - 1) <= 1e-2 + VERBOSE && println("rlm(without leverage correction) : ", coef(m2)) + VERBOSE && println("rlm( with leverage correction) : ", coef(m3)) + ## TODO: find better test + @test isapprox(coef(m2), coef(m3); rtol=0.1) end -end -@testset "linear: leverage correction" begin - m2 = fit( - RobustLinearModel, - form, - data, - est2; - method=:cg, - initial_scale=:mad, - correct_leverage=false, - ) - m3 = fit( - RobustLinearModel, - form, - data, - est2; - method=:cg, - initial_scale=:mad, - correct_leverage=true, - ) - VERBOSE && println("rlm(without leverage correction) : ", coef(m2)) - VERBOSE && println("rlm( with leverage correction) : ", coef(m3)) - ## TODO: find better test - @test isapprox(coef(m2), coef(m3); rtol=0.1) end diff --git a/test/mquantile.jl b/test/mquantile.jl index fe2d71f..e98bf3c 100644 --- a/test/mquantile.jl +++ b/test/mquantile.jl @@ -4,7 +4,7 @@ loss1 = RobustModels.L2Loss() loss2 = RobustModels.TukeyLoss() -@testset "linear: Expectile estimators" begin +@testset "Expectile estimators" begin m = fit( RobustLinearModel, form, @@ -91,7 +91,7 @@ loss2 = RobustModels.TukeyLoss() end -@testset "linear: M-Quantile estimators" begin +@testset "M-Quantile estimators" begin @testset "$(name) estimator" for name in ("Huber", "L1L2", "Fair", "Logcosh", "Arctan") #, "Cauchy", "Geman", "Welsch", "Tukey", "YohaiZamar") VERBOSE && println("\n\t\u25CF M-Quantile Estimator: $name") loss = getproperty(RobustModels, Symbol(name * "Loss"))() diff --git a/test/qreg.jl b/test/qreg.jl index 932dcee..3fd4292 100644 --- a/test/qreg.jl +++ b/test/qreg.jl @@ -1,113 +1,120 @@ -@testset "Quantile regression: low-level function" begin - τs = range(0.1, 0.9; step=0.1) - βs = hcat(map(τ -> RobustModels.interiormethod(X, y, τ)[1], τs)...) - VERBOSE && println("Coefficients: $(vcat(τs', βs))") - @test size(βs) == (size(X, 2), length(τs)) -end +@testset "Quantile regression" begin -@testset "Quantile regression: fit method" begin - τ = 0.5 - # Formula, dense and sparse entry and methods :cg and :chol - @testset "Argument type: $(typeof(A))" for (A, b) in data_tuples - m1 = fit(QuantileRegression, A, b; quantile=τ, verbose=false) - m2 = quantreg(A, b; quantile=τ, verbose=false) - @test_nowarn show(devnull, m2) - @test all(coef(m1) .== coef(m2)) - - # make sure that it is not a TableRegressionModel - @test !isa(m1, TableRegressionModel) - @test !isa(m2, TableRegressionModel) - - # refit - β = copy(coef(m2)) - refit!(m2, y; quantile=τ, verbose=false) - @test all(coef(m2) .== β) - - # interface - @testset "interface method: $(f)" for f in interface_methods - # make sure the method is defined - @test_nowarn robvar = f(m1) - end + @testset "Quantile regression: low-level function" begin + τs = range(0.1, 0.9; step=0.1) + βs = hcat(map(τ -> RobustModels.interiormethod(X, y, τ)[1], τs)...) + VERBOSE && println("Coefficients: $(vcat(τs', βs))") + @test size(βs) == (size(X, 2), length(τs)) + end - # later fit! - m3 = fit(QuantileRegression, A, b; quantile=τ, dofit=false) - @test all(0 .== coef(m3)) - fit!(m3; verbose=false) - @test all(β .== coef(m3)) - - # leverage weights - @test_nowarn refit!(m3; correct_leverage=true) - - # handling of missing values - @testset "Handling of missing values" begin - # check that Missing eltype is routed correctly - if isa(A, FormulaTerm) - if isa(b, NamedTuple) - b_missing = NamedTuple(k => allowmissing(v) for (k, v) in pairs(b)) - elseif isa(b, DataFrame) - b_missing = allowmissing(b) + + @testset "Quantile regression: fit method" begin + τ = 0.5 + # Formula, dense and sparse entry and methods :cg and :chol + @testset "Argument type: $(typeof(A))" for (A, b) in data_tuples + m1 = fit(QuantileRegression, A, b; quantile=τ, verbose=false) + m2 = quantreg(A, b; quantile=τ, verbose=false) + @test_nowarn show(devnull, m2) + @test all(coef(m1) .== coef(m2)) + + # make sure that it is not a TableRegressionModel + @test !isa(m1, TableRegressionModel) + @test !isa(m2, TableRegressionModel) + + # refit + β = copy(coef(m2)) + refit!(m2, y; quantile=τ, verbose=false) + @test all(coef(m2) .== β) + + # interface + @testset "interface method: $(f)" for f in interface_methods + # make sure the method is defined + @test_nowarn robvar = f(m1) + end + + # later fit! + m3 = fit(QuantileRegression, A, b; quantile=τ, dofit=false) + @test all(0 .== coef(m3)) + fit!(m3; verbose=false) + @test all(β .== coef(m3)) + + # leverage weights + @test_nowarn refit!(m3; correct_leverage=true) + + # handling of missing values + @testset "Handling of missing values" begin + # check that Missing eltype is routed correctly + if isa(A, FormulaTerm) + if isa(b, NamedTuple) + b_missing = NamedTuple(k => allowmissing(v) for (k, v) in pairs(b)) + elseif isa(b, DataFrame) + b_missing = allowmissing(b) + else + b_missing = nothing + end + @test_throws ArgumentError fit(QuantileRegression, A, b_missing) + @test_nowarn fit(QuantileRegression, A, b_missing; dropmissing=true) else - b_missing = nothing + A_missing = allowmissing(A) + b_missing = allowmissing(b) + @test_throws ArgumentError fit(QuantileRegression, A_missing, b) + @test_throws ArgumentError fit(QuantileRegression, A, b_missing) + @test_throws ArgumentError fit(QuantileRegression, A_missing, b_missing) + + @test_nowarn fit(QuantileRegression, A_missing, b; dropmissing=true) + @test_nowarn fit(QuantileRegression, A, b_missing; dropmissing=true) + @test_nowarn fit( + QuantileRegression, A_missing, b_missing; dropmissing=true + ) end - @test_throws ArgumentError fit(QuantileRegression, A, b_missing) - @test_nowarn fit(QuantileRegression, A, b_missing; dropmissing=true) - else - A_missing = allowmissing(A) - b_missing = allowmissing(b) - @test_throws ArgumentError fit(QuantileRegression, A_missing, b) - @test_throws ArgumentError fit(QuantileRegression, A, b_missing) - @test_throws ArgumentError fit(QuantileRegression, A_missing, b_missing) - - @test_nowarn fit(QuantileRegression, A_missing, b; dropmissing=true) - @test_nowarn fit(QuantileRegression, A, b_missing; dropmissing=true) - @test_nowarn fit(QuantileRegression, A_missing, b_missing; dropmissing=true) end end end -end -@testset "Quantile regression: different quantiles" begin - τs = range(0.1, 0.9; step=0.1) - m2 = fit(QuantileRegression, form, data; quantile=0.5, verbose=false) - - @testset "$(τ) quantile" for τ in τs - m1 = fit(QuantileRegression, form, data; quantile=τ, verbose=false) - @test_nowarn show(devnull, m1) - β = coef(m1) - res = residuals(m1) - ## The quantile regression line exactly passes through p points, with p number of columns of X. - @test count(x -> isapprox(x, 0; atol=1e-7), res) == length(β) - # @test count(iszero, res) == length(β) - - # refit with new quantile - refit!(m2; quantile=τ) - @test all(coef(m1) .== coef(m2)) - end -end - -@testset "Quantile regression: sparsity estimation" begin - m2 = fit(QuantileRegression, form, data; quantile=0.25, verbose=false) - s = RobustModels.location_variance(m2, false) + @testset "Quantile regression: different quantiles" begin + τs = range(0.1, 0.9; step=0.1) + m2 = fit(QuantileRegression, form, data; quantile=0.5, verbose=false) - τs = range(0.25, 0.75; step=0.25) - @testset "(q, method, kernel): $(τ), $(method), $(kernel)" for τ in τs, - method in (:jones, :bofinger, :hall_sheather), - kernel in (:epanechnikov, :triangle, :window) + @testset "$(τ) quantile" for τ in τs + m1 = fit(QuantileRegression, form, data; quantile=τ, verbose=false) + @test_nowarn show(devnull, m1) + β = coef(m1) + res = residuals(m1) + ## The quantile regression line exactly passes through p points, with p number of columns of X. + @test count(x -> isapprox(x, 0; atol=1e-7), res) == length(β) + # @test count(iszero, res) == length(β) - if τ != m2.τ + # refit with new quantile refit!(m2; quantile=τ) - s = RobustModels.location_variance(m2, false) + @test all(coef(m1) .== coef(m2)) end + end + + @testset "Quantile regression: sparsity estimation" begin + m2 = fit(QuantileRegression, form, data; quantile=0.25, verbose=false) + s = RobustModels.location_variance(m2, false) + + τs = range(0.25, 0.75; step=0.25) + @testset "(q, method, kernel): $(τ), $(method), $(kernel)" for τ in τs, + method in (:jones, :bofinger, :hall_sheather), + kernel in (:epanechnikov, :triangle, :window) - si = RobustModels.location_variance( - m2, false; bw_method=method, α=0.05, kernel=kernel - ) - if method == :jones && kernel == :epanechnikov - @test isapprox(s, si; rtol=1e-4) - else - @test !isapprox(s, si; rtol=1e-4) + if τ != m2.τ + refit!(m2; quantile=τ) + s = RobustModels.location_variance(m2, false) + end + + si = RobustModels.location_variance( + m2, false; bw_method=method, α=0.05, kernel=kernel + ) + if method == :jones && kernel == :epanechnikov + @test isapprox(s, si; rtol=1e-4) + else + @test !isapprox(s, si; rtol=1e-4) + end end end + end diff --git a/test/robustridge.jl b/test/robustridge.jl index 1115a03..aa6496d 100644 --- a/test/robustridge.jl +++ b/test/robustridge.jl @@ -10,110 +10,121 @@ loss2 = RobustModels.TukeyLoss() est1 = MEstimator(loss1) est2 = MEstimator(loss2) -@testset "linear: Ridge M-estimator $(lossname)" for lossname in ("L2", "Huber", "Tukey") - typeloss = getproperty(RobustModels, Symbol(lossname * "Loss")) - l = typeloss() - est = MEstimator(typeloss()) - - # Formula, dense and sparse entry and methods :cg and :chol - @testset "$(typeof(A)),\t$(method)" for (A, b) in data_tuples, method in nopen_methods - - aspace = (method in (:cg, :qr)) ? " " : " " - name = "MEstimator($(typeloss)),\t" - name *= if (A == form) - "formula" - elseif (A == X) - "dense " - else - "sparse " + +@testset "Ridge M-estimator " begin + + #! format: off + @testset "linear: Ridge M-estimator $(lossname)" for lossname in ("L2", "Huber", "Tukey") + #! format: on + typeloss = getproperty(RobustModels, Symbol(lossname * "Loss")) + l = typeloss() + est = MEstimator(typeloss()) + + #! format: off + # Formula, dense and sparse entry and methods :cg and :chol + @testset "$(typeof(A)),\t$(method)" for (A, b) in data_tuples, method in nopen_methods + #! format: on + + aspace = (method in (:cg, :qr)) ? " " : " " + name = "MEstimator($(typeloss)),\t" + name *= if (A == form) + "formula" + elseif (A == X) + "dense " + else + "sparse " + end + name *= (method in (:cg, :qr)) ? ", " : "," + name *= "$(method)" + + kwargs = (; method=method, initial_scale=:L1) + # use the dispersion from GLM to ensure that the loglikelihood is correct + m2 = fit(RobustLinearModel, A, b, est; kwargs...) + m3 = fit(RobustLinearModel, A, b, est; ridgeλ=10, kwargs...) + m4 = fit( + RobustLinearModel, A, b, est; ridgeλ=10, ridgeG=float([0 0; 0 1]), kwargs... + ) + m5 = fit(RobustLinearModel, A, b, est; ridgeλ=0.1, ridgeG=[0, 1], kwargs...) + m6 = rlm(A, b, est; ridgeλ=0.1, ridgeG=[0, 1], dropcollinear=true, kwargs...) + m7 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) + + VERBOSE && println("\n\t\u25CF Estimator: $(name)") + VERBOSE && println(" lm$(aspace) : ", coef(m1)) + VERBOSE && println("rlm($(method)) : ", coef(m2)) + VERBOSE && println("ridge λ=10 rlm3($(method)): ", coef(m3)) + VERBOSE && println("ridge λ=10 rlm4($(method)): ", coef(m4)) + VERBOSE && println("ridge λ=0.1 rlm5($(method)): ", coef(m5)) + VERBOSE && println("ridge λ=0.1 dropcollinear=true rlm6($(method)): ", coef(m6)) + VERBOSE && println("ridge λ=10 βprior=[0,2] rlm7($(method)): ", coef(m7)) + @test isapprox(coef(m3), coef(m4); rtol=1e-5) + + # Test printing the model + @test_nowarn show(devnull, m3) + + # refit + refit!(m5; ridgeλ=10, initial_scale=:L1) + @test isapprox(m5.pred.λ, 10.0; rtol=1e-5) + @test isapprox(coef(m3), coef(m5); rtol=1e-6) + + refit!(m6; ridgeλ=10, initial_scale=:L1) + @test isapprox(m6.pred.λ, 10.0; rtol=1e-5) + @test isapprox(coef(m3), coef(m6); rtol=1e-6) end - name *= (method in (:cg, :qr)) ? ", " : "," - name *= "$(method)" - - kwargs = (; method=method, initial_scale=:L1) - # use the dispersion from GLM to ensure that the loglikelihood is correct - m2 = fit(RobustLinearModel, A, b, est; kwargs...) - m3 = fit(RobustLinearModel, A, b, est; ridgeλ=10, kwargs...) - m4 = fit( - RobustLinearModel, A, b, est; ridgeλ=10, ridgeG=float([0 0; 0 1]), kwargs... - ) - m5 = fit(RobustLinearModel, A, b, est; ridgeλ=0.1, ridgeG=[0, 1], kwargs...) - m6 = rlm(A, b, est; ridgeλ=0.1, ridgeG=[0, 1], dropcollinear=true, kwargs...) - m7 = rlm(A, b, est; ridgeλ=10, ridgeG=[0, 1], βprior=[0.0, 2.0], kwargs...) - - VERBOSE && println("\n\t\u25CF Estimator: $(name)") - VERBOSE && println(" lm$(aspace) : ", coef(m1)) - VERBOSE && println("rlm($(method)) : ", coef(m2)) - VERBOSE && println("ridge λ=10 rlm3($(method)): ", coef(m3)) - VERBOSE && println("ridge λ=10 rlm4($(method)): ", coef(m4)) - VERBOSE && println("ridge λ=0.1 rlm5($(method)): ", coef(m5)) - VERBOSE && println("ridge λ=0.1 dropcollinear=true rlm6($(method)): ", coef(m6)) - VERBOSE && println("ridge λ=10 βprior=[0,2] rlm7($(method)): ", coef(m7)) - @test isapprox(coef(m3), coef(m4); rtol=1e-5) - - # Test printing the model - @test_nowarn show(devnull, m3) - - # refit - refit!(m5; ridgeλ=10, initial_scale=:L1) - @test isapprox(m5.pred.λ, 10.0; rtol=1e-5) - @test isapprox(coef(m3), coef(m5); rtol=1e-6) - - refit!(m6; ridgeλ=10, initial_scale=:L1) - @test isapprox(m6.pred.λ, 10.0; rtol=1e-5) - @test isapprox(coef(m3), coef(m6); rtol=1e-6) end -end -@testset "linear: Ridge L2 estimator methods" begin - m2 = fit(RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1) - m3 = fit(RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1, ridgeλ=1) + @testset "linear: Ridge L2 estimator methods" begin + m2 = fit(RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1) + m3 = fit( + RobustLinearModel, form, data, est1; method=:chol, initial_scale=:L1, ridgeλ=1 + ) - @testset "method: $(f)" for f in interface_methods - # make sure the interfaces for RobustLinearModel are well defined - @test_nowarn f(m3) - end + @testset "method: $(f)" for f in interface_methods + # make sure the interfaces for RobustLinearModel are well defined + @test_nowarn f(m3) + end + + # Check coefficients and dof change + λs = vcat([0], 10 .^ range(-2, 1; length=4)) - # Check coefficients and dof change - λs = vcat([0], 10 .^ range(-2, 1; length=4)) + L = length(coef(m3)) + βs = zeros(L, 5) + dofs = zeros(5) + for (i, λ) in enumerate(λs) + refit!(m3; ridgeλ=λ, initial_scale=:L1) + βs[:, i] = coef(m3) + dofs[i] = dof(m2) + end - L = length(coef(m3)) - βs = zeros(L, 5) - dofs = zeros(5) - for (i, λ) in enumerate(λs) - refit!(m3; ridgeλ=λ, initial_scale=:L1) - βs[:, i] = coef(m3) - dofs[i] = dof(m2) + @test isapprox(dofs[1], dof(m2); rtol=1e-5) + @test isapprox(βs[:, 1], coef(m2); rtol=1e-5) + @test sort(dofs; rev=true) == dofs + @test all(sort(βs[r, :]; rev=true) == βs[r, :] for r in 2:L) # not the intercept end - @test isapprox(dofs[1], dof(m2); rtol=1e-5) - @test isapprox(βs[:, 1], coef(m2); rtol=1e-5) - @test sort(dofs; rev=true) == dofs - @test all(sort(βs[r, :]; rev=true) == βs[r, :] for r in 2:L) # not the intercept -end + @testset "linear: Ridge L2 exact solution" begin + rng = MersenneTwister(seed) -@testset "linear: Ridge L2 exact solution" begin - rng = MersenneTwister(seed) + n = 10_000 + p = 4 + σ = 0.5 + λ = 1000 + Xt = randn(rng, n, p) + βt = randn(rng, p) + yt = Xt * βt + σ * randn(rng, n) - n = 10_000 - p = 4 - σ = 0.5 - λ = 1000 - Xt = randn(rng, n, p) - βt = randn(rng, p) - yt = Xt * βt + σ * randn(rng, n) + vc = inv(Hermitian(Xt'Xt + λ * I(p))) - vc = inv(Hermitian(Xt'Xt + λ * I(p))) + βsol = vc * (Xt' * yt) + dofsol = tr(Xt * vc * Xt') + stdβsol = σ * √(n / (n - p)) * .√(diag(vc * Xt'Xt * vc')) - βsol = vc * (Xt' * yt) - dofsol = tr(Xt * vc * Xt') - stdβsol = σ * √(n / (n - p)) * .√(diag(vc * Xt'Xt * vc')) + m = rlm(Xt, yt, est1; method=:chol, ridgeλ=λ) - m = rlm(Xt, yt, est1; method=:chol, ridgeλ=λ) + @test coef(m) ≈ βsol + @test dof(m) ≈ dofsol + @test vcov(m) ≈ vc + @test_skip stderror(m) ≈ stdβsol + end - @test coef(m) ≈ βsol - @test dof(m) ≈ dofsol - @test vcov(m) ≈ vc - @test_skip stderror(m) ≈ stdβsol end diff --git a/test/runtests.jl b/test/runtests.jl index d5a97ca..1c2d883 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -1,4 +1,4 @@ -if @isdefined(ARGS) +if @isdefined(ARGS) && length(ARGS) > 0 println("Test arguments: $(ARGS)") end diff --git a/test/univariate.jl b/test/univariate.jl index 148f847..05b6f6c 100644 --- a/test/univariate.jl +++ b/test/univariate.jl @@ -2,128 +2,136 @@ using Statistics: median using RobustModels: mean_and_sem, compatdims -funcs = (:mean, :std, :var, :sem, :mean_and_std, :mean_and_var, :mean_and_sem) -seed = 306 +@testset "Univariate statistics" begin -xorig = randn(MersenneTwister(seed), 100) -x = copy(xorig) -x[1:10] .= 50 + funcs = (:mean, :std, :var, :sem, :mean_and_std, :mean_and_var, :mean_and_sem) + seed = 306 + xorig = randn(MersenneTwister(seed), 100) + x = copy(xorig) + x[1:10] .= 50 -@testset "robust univariate statistics: MEstimator{$(lossname)}" for lossname in ( - "L2Loss", "L1Loss", "HuberLoss", "TukeyLoss" -) - typeloss = getproperty(RobustModels, Symbol(lossname)) - est = MEstimator{typeloss}() - @testset "method: $(name)" for name in funcs - func = getproperty(RobustModels, Symbol(name)) - if lossname == "L1Loss" - @test_throws ArgumentError func(est, x) - continue - end + @testset "robust univariate statistics: MEstimator{$(lossname)}" for lossname in ( + "L2Loss", "L1Loss", "HuberLoss", "TukeyLoss" + ) + typeloss = getproperty(RobustModels, Symbol(lossname)) + est = MEstimator{typeloss}() + + @testset "method: $(name)" for name in funcs + func = getproperty(RobustModels, Symbol(name)) + if lossname == "L1Loss" + @test_throws ArgumentError func(est, x) + continue + end - # Compute robust statistics - s = func(est, x) - if lossname == "L2Loss" - @test all(isapprox.(s, func(x); rtol=1e-7)) - elseif lossname == "TukeyLoss" - @test all(s .<= func(x)) - # the estimate is better than when the outliers are removed... - @test_skip all(func(xorig) .<= s) - else - @test all(func(xorig) .<= s .<= func(x)) + # Compute robust statistics + s = func(est, x) + if lossname == "L2Loss" + @test all(isapprox.(s, func(x); rtol=1e-7)) + elseif lossname == "TukeyLoss" + @test all(s .<= func(x)) + # the estimate is better than when the outliers are removed... + @test_skip all(func(xorig) .<= s) + else + @test all(func(xorig) .<= s .<= func(x)) + end end end -end -@testset "robust univariate statistics: Bounded estimator $(typeest)" for typeest in ( - SEstimator, MMEstimator, TauEstimator -) - est = typeest{TukeyLoss}() + @testset "robust univariate statistics: Bounded estimator $(typeest)" for typeest in ( + SEstimator, MMEstimator, TauEstimator + ) + est = typeest{TukeyLoss}() - @testset "method: $(name)" for name in funcs - func = getproperty(RobustModels, Symbol(name)) + @testset "method: $(name)" for name in funcs + func = getproperty(RobustModels, Symbol(name)) - resampling_options = Dict(:rng => MersenneTwister(seed)) - s = func(est, x; resampling_options=resampling_options) - VERBOSE && println( - "statistics $(name): $(round.(s; digits=4)) ≈ $(round.(func(xorig); digits=4)) (with outliers removed)", - ) - @test all(isapprox.(s, func(xorig); rtol=2)) + resampling_options = Dict(:rng => MersenneTwister(seed)) + s = func(est, x; resampling_options=resampling_options) + VERBOSE && println( + "statistics $(name): $(round.(s; digits=4)) ≈ $(round.(func(xorig); digits=4)) (with outliers removed)", + ) + @test all(isapprox.(s, func(xorig); rtol=2)) + end end -end -######################################################################## -## With iterables -######################################################################## + ######################################################################## + ## With iterables + ######################################################################## -d = Base.values(Dict(:a => 1, :b => 2, :c => 3)) -g = (i for i in 1:3) -est = MEstimator{L2Loss}() + d = Base.values(Dict(:a => 1, :b => 2, :c => 3)) + g = (i for i in 1:3) + est = MEstimator{L2Loss}() -@testset "robust univariate statistics: iterable $(typeof(a).name.name)" for a in (d, g) - @test all(isapprox.(mean(est, a), mean(a); rtol=1e-7)) -end + @testset "robust univariate statistics: iterable $(typeof(a).name.name)" for a in (d, g) + @test all(isapprox.(mean(est, a), mean(a); rtol=1e-7)) + end -######################################################################## -## Arrays and dims -######################################################################## - -# use the L2Loss to be directly compared to the non-robust functions -# it should give exactly the same results -est = MEstimator{L2Loss}() - -yorig = randn(MersenneTwister(seed), 306) -y1 = reshape(yorig, (306, 1)) -y2 = reshape(yorig, (1, 306)) -y3 = reshape(yorig, (1, 3, 102)) -y4 = reshape(yorig, (17, 18, 1)) - -@testset "robust univariate statistics: Array size: $(size(a))" for a in (y, y1, y2, y3, y4) - @testset "dims=$(dims)" for dims in (1, 2, (1,), (1, 2), (3, 1), 4, (:)) - ## Mean - m = @test_nowarn mean(est, a; dims=dims) - - # non-robust mean - res = mean(a; dims=dims) - @test all(isapprox.(m, res; rtol=1e-7, nans=true)) - - ## Dispersion: std, var, sem - for disp_name in (:std, :var, :sem) - func = getproperty(RobustModels, Symbol(disp_name)) - - ## Test only the dispersion - s = @test_nowarn func(est, a; dims=dims) - - ## Test `mean_and_ == (mean, )` - func_tup = getproperty(RobustModels, Symbol("mean_and_" * String(disp_name))) - ms = @test_nowarn func_tup(est, a; dims=dims) - @test length(ms) == 2 - @test ms[1] ≈ m - @test ms[2] ≈ s nans = true - - # non-robust dispersion - if dims === Colon() - # apply on a flatten array - res = func(vec(a)) - elseif disp_name == :sem - # `StatsBase.sem` does not allow the `dims` keyword - dims = compatdims(ndims(a), dims) - if isnothing(dims) - continue + ######################################################################## + ## Arrays and dims + ######################################################################## + + # use the L2Loss to be directly compared to the non-robust functions + # it should give exactly the same results + est = MEstimator{L2Loss}() + + yorig = randn(MersenneTwister(seed), 306) + y1 = reshape(yorig, (306, 1)) + y2 = reshape(yorig, (1, 306)) + y3 = reshape(yorig, (1, 3, 102)) + y4 = reshape(yorig, (17, 18, 1)) + + #! format: off + @testset "robust univariate statistics: Array size: $(size(a))" for a in (y, y1, y2, y3, y4) + #! format: on + @testset "dims=$(dims)" for dims in (1, 2, (1,), (1, 2), (3, 1), 4, (:)) + ## Mean + m = @test_nowarn mean(est, a; dims=dims) + + # non-robust mean + res = mean(a; dims=dims) + @test all(isapprox.(m, res; rtol=1e-7, nans=true)) + + ## Dispersion: std, var, sem + for disp_name in (:std, :var, :sem) + func = getproperty(RobustModels, Symbol(disp_name)) + + ## Test only the dispersion + s = @test_nowarn func(est, a; dims=dims) + + ## Test `mean_and_ == (mean, )` + func_tup = getproperty( + RobustModels, Symbol("mean_and_" * String(disp_name)) + ) + ms = @test_nowarn func_tup(est, a; dims=dims) + @test length(ms) == 2 + @test ms[1] ≈ m + @test ms[2] ≈ s nans = true + + # non-robust dispersion + if dims === Colon() + # apply on a flatten array + res = func(vec(a)) + elseif disp_name == :sem + # `StatsBase.sem` does not allow the `dims` keyword + dims = compatdims(ndims(a), dims) + if isnothing(dims) + continue + end + res = mapslices(func, a; dims=dims) + else + # call the non-robust version with dims + res = func(a; dims=dims) end - res = mapslices(func, a; dims=dims) - else - # call the non-robust version with dims - res = func(a; dims=dims) - end - @test all(isapprox.(s, res; rtol=1e-7, nans=true)) + @test all(isapprox.(s, res; rtol=1e-7, nans=true)) + end end end + end diff --git a/test/weights.jl b/test/weights.jl index 5638b07..010797e 100644 --- a/test/weights.jl +++ b/test/weights.jl @@ -1,107 +1,127 @@ -seed = 123987 - -rng = MersenneTwister(seed) - -n = 100 -r = 30 -σ = 10 -xlin = collect(1.0:n) -Xlin = hcat(ones(Float64, n), xlin) -β1lin = [0, -5] -β2lin = [20, -1] - -ylin = vcat(Xlin[1:r, :] * β1lin, Xlin[(r + 1):end, :] * β2lin) -ylin += σ * randn(rng, n) - -wlin = float((1:n) .<= r) -wwlin = 1 .- wlin - - -estimator_list = ( - MEstimator{L2Loss}(), - MEstimator{HuberLoss}(), - MEstimator{TukeyLoss}(), - MMEstimator{TukeyLoss}(), - TauEstimator{TukeyLoss}(), -) - -@testset "weights: $(est)" for est in estimator_list - mtot = rlm(Xlin, ylin, est) - - m1 = rlm(Xlin[1:r, :], ylin[1:r], est) - m1w = rlm(Xlin, ylin, est; wts=wlin) - - m2 = rlm(Xlin[(r + 1):end, :], ylin[(r + 1):end], est) - m2w = rlm(Xlin, ylin, est; wts=wwlin) - - # Check identical values - @testset "$(func)" for func in interface_methods - f1 = func(m1) - f1w = func(m1w) - f2 = func(m2) - f2w = func(m2w) - if func == weights - @test isempty(f1) - @test isempty(f2) - @test !isempty(f1w) - @test !isempty(f2w) - elseif f1 isa AbstractArray && size(f1) != size(f1w) - if ndims(f1) == 1 - @test f1 ≈ f1w[1:r] - @test f2 ≈ f2w[(r + 1):end] - elseif ndims(f1) == 2 - if size(f1, 2) == size(f1w, 2) - @test f1 ≈ f1w[1:r, :] - @test f2 ≈ f2w[(r + 1):end, :] + +@testset "Fit with weights" begin + + seed = 123987 + + rng = MersenneTwister(seed) + + n = 100 + r = 30 + σ = 10 + xlin = collect(1.0:n) + Xlin = hcat(ones(Float64, n), xlin) + β1lin = [0, -5] + β2lin = [20, -1] + + ylin = vcat(Xlin[1:r, :] * β1lin, Xlin[(r + 1):end, :] * β2lin) + ylin += σ * randn(rng, n) + + wlin = float((1:n) .<= r) + wwlin = 1 .- wlin + + datalin = DataFrame(; x=xlin, y=ylin) + formlin = @formula(y ~ 1 + x) + + + estimator_list = ( + MEstimator{L2Loss}(), + MEstimator{HuberLoss}(), + MEstimator{TukeyLoss}(), + MMEstimator{TukeyLoss}(), + TauEstimator{TukeyLoss}(), + ) + + @testset "weights: $(est)" for est in estimator_list + mtot = rlm(Xlin, ylin, est) + + @testset "Model definition: $(model_def)" for model_def in ("arrays", "formula") + with_arrays = model_def == "arrays" + if with_arrays + m1 = rlm(Xlin[1:r, :], ylin[1:r], est) + m1w = rlm(Xlin, ylin, est; wts=wlin) + + m2 = rlm(Xlin[(r + 1):end, :], ylin[(r + 1):end], est) + m2w = rlm(Xlin, ylin, est; wts=wwlin) + + else + m1 = rlm(formlin, datalin[1:r, :], est) + m1w = rlm(formlin, datalin, est; wts=wlin) + + m2 = rlm(formlin, datalin[(r + 1):end, :], est) + m2w = rlm(formlin, datalin, est; wts=wwlin) + end + + # Check identical values + @testset "$(func)" for func in interface_methods + f1 = func(m1) + f1w = func(m1w) + f2 = func(m2) + f2w = func(m2w) + if func == weights + @test isempty(f1) + @test isempty(f2) + @test !isempty(f1w) + @test !isempty(f2w) + elseif f1 isa AbstractArray && size(f1) != size(f1w) + if ndims(f1) == 1 + @test f1 ≈ f1w[1:r] + @test f2 ≈ f2w[(r + 1):end] + elseif ndims(f1) == 2 + if size(f1, 2) == size(f1w, 2) + @test f1 ≈ f1w[1:r, :] + @test f2 ≈ f2w[(r + 1):end, :] + else + @test f1 ≈ f1w[1:r, 1:r] + @test f2 ≈ f2w[(r + 1):end, (r + 1):end] + end + end else - @test f1 ≈ f1w[1:r, 1:r] - @test f2 ≈ f2w[(r + 1):end, (r + 1):end] + @test f1 ≈ f1w + @test f2 ≈ f2w end end - else - @test f1 ≈ f1w - @test f2 ≈ f2w end end -end -@testset "weights: QuantileRegression)" begin - mtot = quantreg(Xlin, ylin) - - m1 = quantreg(Xlin[1:r, :], ylin[1:r]) - m1w = quantreg(Xlin, ylin; wts=wlin) - - m2 = quantreg(Xlin[(r + 1):end, :], ylin[(r + 1):end]) - m2w = quantreg(Xlin, ylin; wts=wwlin) - - # Check identical values - @testset "$(func)" for func in interface_methods - f1 = func(m1) - f1w = func(m1w) - f2 = func(m2) - f2w = func(m2w) - if func == weights - @test isempty(f1) - @test isempty(f2) - @test !isempty(f1w) - @test !isempty(f2w) - elseif f1 isa AbstractArray && size(f1) != size(f1w) - if ndims(f1) == 1 - @test f1 ≈ f1w[1:r] - @test f2 ≈ f2w[(r + 1):end] - elseif ndims(f1) == 2 - if size(f1, 2) == size(f1w, 2) - @test f1 ≈ f1w[1:r, :] - @test f2 ≈ f2w[(r + 1):end, :] - else - @test f1 ≈ f1w[1:r, 1:r] - @test f2 ≈ f2w[(r + 1):end, (r + 1):end] + @testset "weights: QuantileRegression)" begin + mtot = quantreg(Xlin, ylin) + + m1 = quantreg(Xlin[1:r, :], ylin[1:r]) + m1w = quantreg(Xlin, ylin; wts=wlin) + + m2 = quantreg(Xlin[(r + 1):end, :], ylin[(r + 1):end]) + m2w = quantreg(Xlin, ylin; wts=wwlin) + + # Check identical values + @testset "$(func)" for func in interface_methods + f1 = func(m1) + f1w = func(m1w) + f2 = func(m2) + f2w = func(m2w) + if func == weights + @test isempty(f1) + @test isempty(f2) + @test !isempty(f1w) + @test !isempty(f2w) + elseif f1 isa AbstractArray && size(f1) != size(f1w) + if ndims(f1) == 1 + @test f1 ≈ f1w[1:r] + @test f2 ≈ f2w[(r + 1):end] + elseif ndims(f1) == 2 + if size(f1, 2) == size(f1w, 2) + @test f1 ≈ f1w[1:r, :] + @test f2 ≈ f2w[(r + 1):end, :] + else + @test f1 ≈ f1w[1:r, 1:r] + @test f2 ≈ f2w[(r + 1):end, (r + 1):end] + end end + else + @test f1 ≈ f1w + @test f2 ≈ f2w end - else - @test f1 ≈ f1w - @test f2 ≈ f2w end end + end