Tutorial

Read this tutorial to learn how to implement new measures using StatisticalMeasuresBase.jl. Refer to Implementing New Measures and other sections for technical details, bells and whistles.

The goal of this tutorial is to define a Huber loss function, HuberLoss(; delta=1), that consumes vectors of real ground truth observations/predictions, with support for weights, class weights and missing values; and a second measure MultitargetHuberLoss(; delta=1) that consumes matrices, tables or multidimensional arrays. To do so requires only that we implement a single Huber loss for scalars and apply appropriate measure wrappers.

As we'll see, the most important wrapper to understand is the multimeasure wrapper.

The basics

Let's start by defining a scalar Huber loss function:

valley(a, delta=1) = abs(a) <= delta ? a^2/2 : delta*(abs(a) - delta/2)
huber(ŷ, y, delta=1) = valley(ŷ - y, delta)

huber (generic function with 2 methods)

Huber loss type for scalar input

Since the Huber loss has a parameter, we create a struct, and make instances callable:

using StatisticalMeasuresBase
import StatisticalMeasuresBase as API

struct HuberLossOnScalars{T<:Real}
    delta::T
end
(measure::HuberLossOnScalars)(ŷ, y) = huber(ŷ, y, measure.delta)
API.is_measure(::HuberLossOnScalars) = true # recommended if `HuberLossOnScalars` is public-facing

measure = HuberLossOnScalars(0.5)
@assert measure(7, 4) == huber(7, 4, 0.5)

Huber loss for vector input

To get a Huber loss for vectors, we wrap the scalar measure using multimeasure, which broadcasts over MLUtils.eachobs((ŷ, y)) in call invocations. We also include a preliminary wrapping to support missing values:

HuberLoss(delta) =
    multimeasure(supports_missings_measure(HuberLossOnScalars(delta)))
HuberLoss(; delta=1) = HuberLoss(delta)

HuberLoss (generic function with 2 methods)

By default, a new measure has Mean() as its external aggregation mode, which you can inspect using the trait StatisticalMeasuresBase.external_aggregation_mode(measure). Wrappers, like support_missings_measure, forward such traits to the wrapped measure. The multimeasure wrapper uses MLUtils.jl's eachobs method to broadcast the atomic measure over oberservations and, by default, aggregates the result using the atomic measure's aggregation mode.

Demonstration

using Statistics
ŷ = [1, 2, missing]
y = [7, 4, 6]
weights = [1, 3, 2]
class_weights = Dict(7=>3.0, 4=>4.0, 6=>2.0)

wts = weights .* (class_weights[η] for η in y)
@assert HuberLoss(1)(ŷ, y, weights, class_weights) ≈
    sum(wts[1:2] .* huber.(ŷ[1:2], y[1:2])) / length(wts)

Note the divisor length(weights) on the last line: The weighted measurement is not literally the weighted mean but the weighted mean scaled by the average value of the weights. To get the bone fide weighted mean, use multimeasure(..., mode=IMean()) instead. Another option is Sum(); see StatisticalMeasuresBase.AggregationMode for other options.

Multi-target Huber loss

Here's the Huber loss for multi-targets (matrices or tables) which includes strict argument checks, and works for higher dimensional arrays as well (argument checks can be removed with unfussy(measure)):

MultitargetHuberLoss(args...; kwargs...) =
    multimeasure(HuberLoss(args...; kwargs...), transform = vec∘collect) |> fussy_measure

MultitargetHuberLoss (generic function with 1 method)

Demonstration

Multi-targets (as matrices):

y       = rand(3, 20)
ŷ       = rand(3, 20)
weights = rand(20)
ms = weights .* vec(mean(huber.(ŷ, y), dims=1))
@assert measurements(MultitargetHuberLoss(), ŷ, y, weights) ≈ ms
@assert MultitargetHuberLoss()(ŷ, y, weights) ≈  sum(ms)/length(weights)

Note the use of the measurements method, to return one measurement per observation, instead of an aggregate.

Mutli-targets (as tables):

using DataFrames
df    = DataFrame(y', :auto)
df̂    = DataFrame(ŷ', :auto)
@assert MultitargetHuberLoss()(df̂, df, weights) ≈ MultitargetHuberLoss()(df̂, df, weights)

Multi-dimensional arrays:

y    = rand(2, 3, 20)
ŷ    = rand(2, 3, 20)
weights = rand(20)
ms = weights .* vec(mean(huber.(ŷ, y), dims=[1, 2]))
@assert measurements(MultitargetHuberLoss(), ŷ, y, weights) ≈ ms
@assert MultitargetHuberLoss()(ŷ, y, weights) ≈  sum(ms)/length(weights)

Invalid arguments:

class_weights = Dict()
MultitargetHuberLoss()(ŷ, y, class_weights)
ERROR: ArgumentError: Class pool or value set is not compatible with the class_weight dictionary keys.

Overloading traits

Here's how to overload measure traits to make additional promises of behavior.

using ScientificTypesBase
@trait HuberLossOnScalars orientation=Loss()

typeof(HuberLoss())

StatisticalMeasuresBase.Multimeasure{StatisticalMeasuresBase.SupportsMissingsMeasure{Main.HuberLossOnScalars{Int64}}, Nothing, Mean, typeof(identity)}

const HuberLossType = API.Multimeasure{
    <:API.SupportsMissingsMeasure{
    <:HuberLossOnScalars
}}
@trait(
    HuberLossType,
    observation_scitype = Union{Continuous,Missing},
    human_name = "Huber loss",
)

API.observation_scitype(HuberLoss())

Union{Missing, ScientificTypesBase.Continuous}

const MultitargetHuberLossType = API.FussyMeasure{<:API.Multimeasure{<:HuberLossType}}
@trait(
    MultitargetHuberLossType,
    observation_scitype = AbstractArray{<:Union{Continuous,Missing}},
    can_consume_tables = true,
    human_name = "multi-target Huber loss",
)

Note

While most wrappers forward atomic measure traits appropriately, StatisticalMeasuresBase.observation_scitype(measure) (default=Union{}) and StatisticalMeasuresBase.can_consume_tables(measure) (default=false) must be explicitly overloaded for measures wrapped using multimeasure.

Improving display of measures

Because most useful measures are wraps of atomic measures, they don't display well:

HuberLoss()

multimeasure(Main.HuberLossOnScalars{Int64}(1))

This can be improved using the @fix_show macro:

@fix_show HuberLoss::HuberLossType
HuberLoss()

HuberLoss(
  delta = 1)

Macro shortcut (experimental)

The entire workflow above is equivalent to the code block below, except that HuberLoss is also fussy and MultitargetHuberLoss gets an additional keyword argument atomic_weights, for specifying weights per component of the vector-valued target (per column if y is a table). The loss functions also accept nothing for weights or class weights.

You will need to try this in a new Julia session or change the loss name. See @combination for details.

valley(a, delta=1) = abs(a) <= delta ? a^2/2 : delta*(abs(a) - delta/2)
huber(ŷ, y, delta=1) = valley(ŷ - y, delta)

using StatisticalMeasuresBase
import StatisticalMeasuresBase as API
using ScientificTypesBase

@combination(
    HuberLoss(; delta=1) = multimeasure(huber),
    orientation =  Loss(),
    observation_scitype = Continuous,
    human_name = "Huber loss",
)

Demonstration

n = 10
y = vcat(fill(1, n)', fill(1, n)', fill(1, n)')
ŷ = vcat(fill(1, n)', fill(2, n)', fill(3, n)')
ŷ - y

3×10 Matrix{Int64}:
 0  0  0  0  0  0  0  0  0  0
 1  1  1  1  1  1  1  1  1  1
 2  2  2  2  2  2  2  2  2  2

delta = 10
valley(1, delta)

0.5

valley(2, delta)

2.0

measure = MultitargetHuberLoss(delta=10, atomic_weights=[5, 4, 1])
@assert measure(ŷ, y, fill(100, n)) ≈ (5*0 + 4*0.5 + 1*2.0)/3*100

API.observation_scitype(measure)

AbstractArray{<:Union{Missing, ScientificTypesBase.Continuous}}