Skip to content

Commit 215685a

Browse files
committed
improve coverage losses.jl
1 parent 9796dad commit 215685a

File tree

2 files changed

+89
-87
lines changed

2 files changed

+89
-87
lines changed

src/losses.jl

Lines changed: 53 additions & 67 deletions
Original file line numberDiff line numberDiff line change
@@ -29,18 +29,18 @@ function weight end
2929
function estimator_norm end
3030

3131
"The limit at ∞ of the loss function. Used for scale estimation of bounded loss."
32-
estimator_bound(f::LossFunction) = Inf
32+
estimator_bound(::LossFunction) = Inf
3333

3434
"The tuning constant of the loss function, can be optimized to get efficient or robust estimates."
3535
tuning_constant(loss::L) where {L<:LossFunction} = loss.c
3636

3737
"Boolean if the estimator or loss function is convex"
38-
isconvex(f::LossFunction) = false
39-
isconvex(f::ConvexLossFunction) = true
38+
isconvex(::LossFunction) = false
39+
isconvex(::ConvexLossFunction) = true
4040

4141
"Boolean if the estimator or loss function is bounded"
42-
isbounded(f::LossFunction) = false
43-
isbounded(f::BoundedLossFunction) = true
42+
isbounded(::LossFunction) = false
43+
isbounded(::BoundedLossFunction) = true
4444

4545
"The tuning constant associated to the loss that gives an efficient M-estimator."
4646
estimator_high_breakdown_point_constant(::Type{L}) where {L<:LossFunction} = 1
@@ -205,10 +205,10 @@ struct HuberLoss <: ConvexLossFunction
205205
HuberLoss() = new(estimator_high_efficiency_constant(HuberLoss))
206206
end
207207

208-
_rho(l::HuberLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : (abs(r) - 1 / 2)
209-
_psi(l::HuberLoss, r::Real) = (abs(r) <= 1) ? r : sign(r)
210-
_psider(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0)
211-
_weight(l::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : 1 / abs(r)
208+
_rho(::HuberLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : (abs(r) - 1 / 2)
209+
_psi(::HuberLoss, r::Real) = (abs(r) <= 1) ? r : sign(r)
210+
_psider(::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0)
211+
_weight(::HuberLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : 1 / abs(r)
212212
function estimator_values(l::HuberLoss, r::Real)
213213
rr = abs(r)
214214
if rr <= l.c
@@ -233,10 +233,10 @@ struct L1L2Loss <: ConvexLossFunction
233233
L1L2Loss(c::Real) = new(c)
234234
L1L2Loss() = new(estimator_high_efficiency_constant(L1L2Loss))
235235
end
236-
_rho(l::L1L2Loss, r::Real) = (sqrt(1 + r^2) - 1)
237-
_psi(l::L1L2Loss, r::Real) = r / sqrt(1 + r^2)
238-
_psider(l::L1L2Loss, r::Real) = 1 / (1 + r^2)^(3 / 2)
239-
_weight(l::L1L2Loss, r::Real) = 1 / sqrt(1 + r^2)
236+
_rho(::L1L2Loss, r::Real) = (sqrt(1 + r^2) - 1)
237+
_psi(::L1L2Loss, r::Real) = r / sqrt(1 + r^2)
238+
_psider(::L1L2Loss, r::Real) = 1 / (1 + r^2)^(3 / 2)
239+
_weight(::L1L2Loss, r::Real) = 1 / sqrt(1 + r^2)
240240
function estimator_values(l::L1L2Loss, r::Real)
241241
sqr = sqrt(1 + (r / l.c)^2)
242242
return ((sqr - 1), r / sqr, 1 / sqr)
@@ -257,10 +257,10 @@ struct FairLoss <: ConvexLossFunction
257257
FairLoss(c::Real) = new(c)
258258
FairLoss() = new(estimator_high_efficiency_constant(FairLoss))
259259
end
260-
_rho(l::FairLoss, r::Real) = abs(r) - log(1 + abs(r))
261-
_psi(l::FairLoss, r::Real) = r / (1 + abs(r))
262-
_psider(l::FairLoss, r::Real) = 1 / (1 + abs(r))^2
263-
_weight(l::FairLoss, r::Real) = 1 / (1 + abs(r))
260+
_rho(::FairLoss, r::Real) = abs(r) - log(1 + abs(r))
261+
_psi(::FairLoss, r::Real) = r / (1 + abs(r))
262+
_psider(::FairLoss, r::Real) = 1 / (1 + abs(r))^2
263+
_weight(::FairLoss, r::Real) = 1 / (1 + abs(r))
264264
function estimator_values(l::FairLoss, r::Real)
265265
ir = 1 / (1 + abs(r / l.c))
266266
return (abs(r) / l.c + log(ir), r * ir, ir)
@@ -279,10 +279,10 @@ struct LogcoshLoss <: ConvexLossFunction
279279
LogcoshLoss(c::Real) = new(c)
280280
LogcoshLoss() = new(estimator_high_efficiency_constant(LogcoshLoss))
281281
end
282-
_rho(l::LogcoshLoss, r::Real) = log(cosh(r))
283-
_psi(l::LogcoshLoss, r::Real) = tanh(r)
284-
_psider(l::LogcoshLoss, r::Real) = 1 / (cosh(r))^2
285-
_weight(l::LogcoshLoss, r::Real) = (abs(r) < DELTA) ? (1 - (r)^2 / 3) : tanh(r) / r
282+
_rho(::LogcoshLoss, r::Real) = log(cosh(r))
283+
_psi(::LogcoshLoss, r::Real) = tanh(r)
284+
_psider(::LogcoshLoss, r::Real) = 1 / (cosh(r))^2
285+
_weight(::LogcoshLoss, r::Real) = (abs(r) < DELTA) ? (1 - (r)^2 / 3) : tanh(r) / r
286286
function estimator_values(l::LogcoshLoss, r::Real)
287287
tr = l.c * tanh(r / l.c)
288288
rr = abs(r / l.c)
@@ -302,10 +302,10 @@ struct ArctanLoss <: ConvexLossFunction
302302
ArctanLoss(c::Real) = new(c)
303303
ArctanLoss() = new(estimator_high_efficiency_constant(ArctanLoss))
304304
end
305-
_rho(l::ArctanLoss, r::Real) = r * atan(r) - 1 / 2 * log(1 + r^2)
306-
_psi(l::ArctanLoss, r::Real) = atan(r)
307-
_psider(l::ArctanLoss, r::Real) = 1 / (1 + r^2)
308-
_weight(l::ArctanLoss, r::Real) = (abs(r) < DELTA) ? (1 - r^2 / 3) : atan(r) / r
305+
_rho(::ArctanLoss, r::Real) = r * atan(r) - 1 / 2 * log(1 + r^2)
306+
_psi(::ArctanLoss, r::Real) = atan(r)
307+
_psider(::ArctanLoss, r::Real) = 1 / (1 + r^2)
308+
_weight(::ArctanLoss, r::Real) = (abs(r) < DELTA) ? (1 - r^2 / 3) : atan(r) / r
309309
function estimator_values(l::ArctanLoss, r::Real)
310310
ar = atan(r / l.c)
311311
rr = abs(r / l.c)
@@ -332,13 +332,13 @@ struct CatoniWideLoss <: ConvexLossFunction
332332
CatoniWideLoss() = new(estimator_high_efficiency_constant(CatoniWideLoss))
333333
end
334334

335-
function _rho(l::CatoniWideLoss, r::Real)
335+
function _rho(::CatoniWideLoss, r::Real)
336336
ar = abs(r)
337337
return (1 + ar) * log(1 + ar + r^2 / 2) - 2 * ar + 2 * atan(1 + r) - π / 2
338338
end
339-
_psi(l::CatoniWideLoss, r::Real) = sign(r) * log(1 + abs(r) + r^2 / 2)
340-
_psider(l::CatoniWideLoss, r::Real) = (1 + abs(r)) / (1 + abs(r) + r^2 / 2)
341-
function _weight(l::CatoniWideLoss, r::Real)
339+
_psi(::CatoniWideLoss, r::Real) = sign(r) * log(1 + abs(r) + r^2 / 2)
340+
_psider(::CatoniWideLoss, r::Real) = (1 + abs(r)) / (1 + abs(r) + r^2 / 2)
341+
function _weight(::CatoniWideLoss, r::Real)
342342
return (abs(r) < DELTA) ? oftype(r, 1) : log(1 + abs(r) + r^2 / 2) / abs(r)
343343
end
344344
function estimator_values(l::CatoniWideLoss, r::Real)
@@ -368,26 +368,26 @@ struct CatoniNarrowLoss <: ConvexLossFunction
368368
CatoniNarrowLoss() = new(estimator_high_efficiency_constant(CatoniNarrowLoss))
369369
end
370370

371-
function _rho(l::CatoniNarrowLoss, r::Real)
371+
function _rho(::CatoniNarrowLoss, r::Real)
372372
ar = abs(r)
373373
if ar >= 1
374374
return (ar - 1) * log(2) + 2 - π / 2
375375
end
376376
return (1 - ar) * log(1 - ar + r^2 / 2) + 2 * ar + 2 * atan(1 - ar) - π / 2
377377
end
378-
function _psi(l::CatoniNarrowLoss, r::Real)
378+
function _psi(::CatoniNarrowLoss, r::Real)
379379
if abs(r) <= 1
380380
return -sign(r) * log(1 - abs(r) + r^2 / 2)
381381
end
382382
return sign(r) * log(2)
383383
end
384-
function _psider(l::CatoniNarrowLoss, r::Real)
384+
function _psider(::CatoniNarrowLoss, r::Real)
385385
if abs(r) <= 1
386386
return (1 - abs(r)) / (1 - abs(r) + r^2 / 2)
387387
end
388388
return 0
389389
end
390-
function _weight(l::CatoniNarrowLoss, r::Real)
390+
function _weight(::CatoniNarrowLoss, r::Real)
391391
if r == 0
392392
return oftype(r, 1)
393393
elseif abs(r) <= 1
@@ -422,18 +422,16 @@ struct CauchyLoss <: LossFunction
422422
CauchyLoss(c::Real) = new(c)
423423
CauchyLoss() = new(estimator_high_efficiency_constant(CauchyLoss))
424424
end
425-
_rho(l::CauchyLoss, r::Real) = log(1 + r^2) # * 1/2 # remove factor 1/2 so the loss has a norm
426-
_psi(l::CauchyLoss, r::Real) = r / (1 + r^2)
427-
_psider(l::CauchyLoss, r::Real) = (1 - r^2) / (1 + r^2)^2
428-
_weight(l::CauchyLoss, r::Real) = 1 / (1 + r^2)
425+
_rho(::CauchyLoss, r::Real) = log(1 + r^2) # * 1/2 # remove factor 1/2 so the loss has a norm
426+
_psi(::CauchyLoss, r::Real) = r / (1 + r^2)
427+
_psider(::CauchyLoss, r::Real) = (1 - r^2) / (1 + r^2)^2
428+
_weight(::CauchyLoss, r::Real) = 1 / (1 + r^2)
429429
function estimator_values(l::CauchyLoss, r::Real)
430430
ir = 1 / (1 + (r / l.c)^2)
431431
return (-log(ir), r * ir, ir)
432432
end
433433
estimator_norm(l::CauchyLoss) = l.c * π
434434
estimator_bound(l::CauchyLoss) = 1.0
435-
isconvex(::CauchyLoss) = false
436-
isbounded(::CauchyLoss) = false
437435

438436
estimator_high_efficiency_constant(::Type{CauchyLoss}) = 2.385
439437
estimator_high_breakdown_point_constant(::Type{CauchyLoss}) = 1.468
@@ -449,19 +447,17 @@ struct GemanLoss <: BoundedLossFunction
449447
GemanLoss(c::Real) = new(c)
450448
GemanLoss() = new(estimator_high_efficiency_constant(GemanLoss))
451449
end
452-
_rho(l::GemanLoss, r::Real) = 1 / 2 * r^2 / (1 + r^2)
453-
_psi(l::GemanLoss, r::Real) = r / (1 + r^2)^2
454-
_psider(l::GemanLoss, r::Real) = (1 - 3 * r^2) / (1 + r^2)^3
455-
_weight(l::GemanLoss, r::Real) = 1 / (1 + r^2)^2
450+
_rho(::GemanLoss, r::Real) = 1 / 2 * r^2 / (1 + r^2)
451+
_psi(::GemanLoss, r::Real) = r / (1 + r^2)^2
452+
_psider(::GemanLoss, r::Real) = (1 - 3 * r^2) / (1 + r^2)^3
453+
_weight(::GemanLoss, r::Real) = 1 / (1 + r^2)^2
456454
function estimator_values(l::GemanLoss, r::Real)
457455
rr2 = (r / l.c)^2
458456
ir = 1 / (1 + rr2)
459457
return (1 / 2 * rr2 * ir, r * ir^2, ir^2)
460458
end
461459
estimator_norm(::GemanLoss) = Inf
462460
estimator_bound(::GemanLoss) = 1 / 2
463-
isconvex(::GemanLoss) = false
464-
isbounded(::GemanLoss) = true
465461

466462
estimator_high_efficiency_constant(::Type{GemanLoss}) = 3.787
467463
estimator_high_breakdown_point_constant(::Type{GemanLoss}) = 0.61200
@@ -478,19 +474,17 @@ struct WelschLoss <: BoundedLossFunction
478474
WelschLoss(c::Real) = new(c)
479475
WelschLoss() = new(estimator_high_efficiency_constant(WelschLoss))
480476
end
481-
_rho(l::WelschLoss, r::Real) = -1 / 2 * Base.expm1(-r^2)
482-
_psi(l::WelschLoss, r::Real) = r * exp(-r^2)
483-
_psider(l::WelschLoss, r::Real) = (1 - 2 * r^2) * exp(-r^2)
484-
_weight(l::WelschLoss, r::Real) = exp(-r^2)
477+
_rho(::WelschLoss, r::Real) = -1 / 2 * Base.expm1(-r^2)
478+
_psi(::WelschLoss, r::Real) = r * exp(-r^2)
479+
_psider(::WelschLoss, r::Real) = (1 - 2 * r^2) * exp(-r^2)
480+
_weight(::WelschLoss, r::Real) = exp(-r^2)
485481
function estimator_values(l::WelschLoss, r::Real)
486482
rr2 = (r / l.c)^2
487483
er = exp(-rr2)
488484
return (-1 / 2 * Base.expm1(-rr2), r * er, er)
489485
end
490486
estimator_norm(::WelschLoss) = Inf
491487
estimator_bound(::WelschLoss) = 1 / 2
492-
isconvex(::WelschLoss) = false
493-
isbounded(::WelschLoss) = true
494488

495489
estimator_high_efficiency_constant(::Type{WelschLoss}) = 2.985
496490
estimator_high_breakdown_point_constant(::Type{WelschLoss}) = 0.8165
@@ -507,18 +501,16 @@ struct TukeyLoss <: BoundedLossFunction
507501
TukeyLoss(c::Real) = new(c)
508502
TukeyLoss() = new(estimator_high_efficiency_constant(TukeyLoss))
509503
end
510-
_rho(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 / 6 * (1 - (1 - r^2)^3) : 1 / 6
511-
_psi(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? r * (1 - r^2)^2 : oftype(r, 0)
512-
_psider(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 - 6 * r^2 + 5 * r^4 : oftype(r, 0)
513-
_weight(l::TukeyLoss, r::Real) = (abs(r) <= 1) ? (1 - r^2)^2 : oftype(r, 0)
504+
_rho(::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 / 6 * (1 - (1 - r^2)^3) : 1 / 6
505+
_psi(::TukeyLoss, r::Real) = (abs(r) <= 1) ? r * (1 - r^2)^2 : oftype(r, 0)
506+
_psider(::TukeyLoss, r::Real) = (abs(r) <= 1) ? 1 - 6 * r^2 + 5 * r^4 : oftype(r, 0)
507+
_weight(::TukeyLoss, r::Real) = (abs(r) <= 1) ? (1 - r^2)^2 : oftype(r, 0)
514508
function estimator_values(l::TukeyLoss, r::Real)
515509
pr = (abs(r) <= l.c) * (1 - (r / l.c)^2)
516510
return (1 / 6 * (1 - pr^3), r * pr^2, pr^2)
517511
end
518512
estimator_norm(::TukeyLoss) = Inf
519513
estimator_bound(::TukeyLoss) = 1 / 6
520-
isconvex(::TukeyLoss) = false
521-
isbounded(::TukeyLoss) = true
522514

523515
estimator_high_efficiency_constant(::Type{TukeyLoss}) = 4.685
524516
estimator_high_breakdown_point_constant(::Type{TukeyLoss}) = 1.5476
@@ -590,8 +582,6 @@ function estimator_values(l::YohaiZamarLoss, r::Real)
590582
end
591583
estimator_norm(::YohaiZamarLoss) = Inf
592584
estimator_bound(::YohaiZamarLoss) = 1
593-
isconvex(::YohaiZamarLoss) = false
594-
isbounded(::YohaiZamarLoss) = true
595585

596586
estimator_high_efficiency_constant(::Type{YohaiZamarLoss}) = 3.1806
597587
estimator_high_breakdown_point_constant(::Type{YohaiZamarLoss}) = 1.2139
@@ -607,9 +597,9 @@ struct HardThresholdLoss <: BoundedLossFunction
607597
HardThresholdLoss(c::Real) = new(c)
608598
HardThresholdLoss() = new(estimator_high_efficiency_constant(HardThresholdLoss))
609599
end
610-
_rho(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : oftype(r, 1 / 2)
611-
_psi(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r : oftype(r, 0)
612-
function _psider(l::HardThresholdLoss, r::Real)
600+
_rho(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r^2 / 2 : oftype(r, 1 / 2)
601+
_psi(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? r : oftype(r, 0)
602+
function _psider(::HardThresholdLoss, r::Real)
613603
if abs(r) <= 1
614604
return oftype(r, 1)
615605
elseif abs(r) < 1 + DELTA
@@ -618,7 +608,7 @@ function _psider(l::HardThresholdLoss, r::Real)
618608
return oftype(r, 0)
619609
end
620610
end
621-
_weight(l::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0)
611+
_weight(::HardThresholdLoss, r::Real) = (abs(r) <= 1) ? oftype(r, 1) : oftype(r, 0)
622612
function estimator_values(l::HardThresholdLoss, r::Real)
623613
ar = abs(r / l.c)
624614
if ar <= 1
@@ -629,8 +619,6 @@ function estimator_values(l::HardThresholdLoss, r::Real)
629619
end
630620
estimator_norm(::HardThresholdLoss) = Inf
631621
estimator_bound(::HardThresholdLoss) = 1 / 2
632-
isconvex(::HardThresholdLoss) = false
633-
isbounded(::HardThresholdLoss) = true
634622

635623
## Computed with `find_zero(c -> I1(c) - 0.95, 2.8, Order1())` after simplification (fun_eff = I1)
636624
estimator_high_efficiency_constant(::Type{HardThresholdLoss}) = 2.795
@@ -709,8 +697,6 @@ end
709697

710698
estimator_norm(::HampelLoss) = Inf
711699
estimator_bound(l::HampelLoss) = (l.ν1 + l.ν2 - 1) / 2
712-
isconvex(::HampelLoss) = false
713-
isbounded(::HampelLoss) = true
714700

715701
# Values of `c` for (ν1, ν2) = (2, 4)
716702
estimator_high_efficiency_constant(::Type{HampelLoss}) = 1.382

test/estimators.jl

Lines changed: 36 additions & 20 deletions
Original file line numberDiff line numberDiff line change
@@ -89,35 +89,51 @@ emp_norm(l::LossFunction) = 2 * quadgk(x -> exp(-RobustModels.rho(l, x)), 0, Inf
8989

9090
# Only for LossFunction
9191
if isnothing(estimator)
92-
if !isbounded(est)
93-
@testset "Estimator norm: $(name)" begin
94-
@test emp_norm(est) RobustModels.estimator_norm(est) rtol = 1e-5
92+
@testset "Loss methods: $(name)" begin
93+
@test loss(est) == est
94+
95+
# estimator_bound
96+
if !isconvex(est)
97+
@test isfinite(RobustModels.estimator_bound(est))
98+
else
99+
@test !isfinite(RobustModels.estimator_bound(est))
95100
end
96-
end
97101

98-
if !in(name, ("L2", "L1"))
99-
@testset "Estimator high efficiency: $(name)" begin
100-
vopt = estimator_high_efficiency_constant(typest)
101-
if name != "HardThreshold"
102-
v = efficiency_tuning_constant(typest; eff=0.95, c0=0.9 * vopt)
103-
@test isapprox(v, vopt; rtol=1e-3)
102+
# tuning_constant
103+
@test isfinite(RobustModels.tuning_constant(est))
104+
105+
@testset "Estimator norm: $(name)" begin
106+
if !isbounded(est)
107+
@test emp_norm(est) RobustModels.estimator_norm(est) rtol = 1e-5
108+
else
109+
@test !isfinite(RobustModels.estimator_norm(est))
104110
end
105111
end
106-
end
107112

108-
if isbounded(est)
109-
@testset "Estimator high breakdown point: $(name)" begin
110-
vopt = estimator_high_breakdown_point_constant(typest)
111-
v = breakdown_point_tuning_constant(typest; bp=0.5, c0=1.1 * vopt)
112-
@test isapprox(v, vopt; rtol=1e-3)
113+
if !in(name, ("L2", "L1"))
114+
@testset "Estimator high efficiency: $(name)" begin
115+
vopt = estimator_high_efficiency_constant(typest)
116+
if name != "HardThreshold"
117+
v = efficiency_tuning_constant(typest; eff=0.95, c0=0.9 * vopt)
118+
@test isapprox(v, vopt; rtol=1e-3)
119+
end
120+
end
113121
end
114122

115-
@testset "τ-Estimator high efficiency: $(name)" begin
116-
vopt = estimator_tau_efficient_constant(typest)
117-
if name != "HardThreshold"
118-
v = tau_efficiency_tuning_constant(typest; eff=0.95, c0=1.1 * vopt)
123+
if isbounded(est)
124+
@testset "Estimator high breakdown point: $(name)" begin
125+
vopt = estimator_high_breakdown_point_constant(typest)
126+
v = breakdown_point_tuning_constant(typest; bp=0.5, c0=1.1 * vopt)
119127
@test isapprox(v, vopt; rtol=1e-3)
120128
end
129+
130+
@testset "τ-Estimator high efficiency: $(name)" begin
131+
vopt = estimator_tau_efficient_constant(typest)
132+
if name != "HardThreshold"
133+
v = tau_efficiency_tuning_constant(typest; eff=0.95, c0=1.1 * vopt)
134+
@test isapprox(v, vopt; rtol=1e-3)
135+
end
136+
end
121137
end
122138
end
123139
end

0 commit comments

Comments
 (0)