Forecast input modes

ScoringRules.jl accepts four kinds of forecast input. The same scoring function names (crps, logs, dss) dispatch on the type of the first argument, so you rarely need to think about which method you are calling.

Parametric forecasts (Distributions.jl)

Pass any UnivariateDistribution from Distributions.jl as the first argument. The package uses a closed-form CRPS where one exists; it falls back to adaptive quadrature otherwise.

using ScoringRules, Distributions

crps(Normal(1.0, 2.0), 0.5)

0.5169996257988083

crps(Gamma(2.0, 1.5), 3.0)

0.49902339883935243

logs(LogNormal(0.0, 0.5), 1.2)

0.47459520958222445

dss(TDist(5.0), 1.0)

1.1108256237659906

Truncated and censored variants from Distributions.jl are also supported:

d_trunc = truncated(Normal(0.0, 1.0); lower=0.0)
crps(d_trunc, 0.7)

0.14313834014692794

d_cens = Distributions.censored(Normal(0.0, 1.0); lower=0.0)
crps(d_cens, 0.7)

0.30472168144590944

Ensemble / sample forecasts

When a parametric form is not available, pass a vector of draws. The CRPS is computed via the empirical distribution function (EDF) approximation by default; a kernel density estimate is available via method = :kde.

draws = randn(1000)

crps(draws, 0.5)

0.32020030190906623

crps(draws, 0.5; method = :kde)   # Gaussian KDE with Silverman bandwidth

0.3267087546419093

logs(draws, 0.5)   # always uses KDE

1.0253268536783762

dss(draws, 0.5)    # uses sample mean and population variance

0.2995925181818994

Weighted ensemble members

Importance weights can be passed to crps (EDF method) and dss:

w = abs.(randn(1000)) .+ 0.01   # arbitrary positive weights (normalised internally)
crps(draws, 0.5; w = w)

0.34756829706297426

dss(draws, 0.5; w = w)    # weighted mean and variance

0.36398821700022965

Moment-based forecasts

When you know only the mean and variance of the forecast — not its full distribution — use dss_moments. The Dawid–Sebastiani score depends on the forecast through its first two moments alone, so no distribution needs to be specified.

dss_moments(0.5, 1.0, 4.0)   # observation y=0.5, mean=1.0, variance=4.0

1.4487943611198906

The error-spread score of Christensen, Moroz and Palmer (2015) also takes moment forecasts, extending them with the skewness. It assesses whether the spread and skewness of an ensemble forecast are consistent with its error:

ess_moments(0.5, 1.0, 4.0, 0.3)   # as above, with skewness 0.3

11.902500000000002

This is useful when forecasts arrive as published summary statistics rather than as distributional objects.

Quantile forecasts

For quantile-format forecasts supply vectors of quantile levels and the corresponding quantile values:

levels = [0.1, 0.25, 0.5, 0.75, 0.9]
qs     = [-1.28, -0.67, 0.0, 0.67, 1.28]   # standard-normal quantiles

scores = quantile_score(levels, qs, 0.5)

5-element Vector{Float64}:
 0.17800000000000002
 0.2925
 0.25
 0.04250000000000001
 0.07799999999999999

mean(scores)   # mean quantile score

0.1682

For a single quantile level from an ensemble:

quantile_score(draws, 0.5; alpha = 0.9)   # 90th-percentile score

0.0792247201019944

Prediction intervals can be scored with interval_score:

# 90% prediction interval for N(0,1): approximately (−1.645, 1.645)
interval_score(-1.645, 1.645, 0.5, 0.9)

3.2900000000000014

interval_score(draws, 0.5; level = 0.9)   # same, estimated from the ensemble

3.3634020012093537

Multivariate ensembles

For a d-dimensional ensemble of m members, represent the forecast as a d × m matrix (each column is one member) and the observation as a length-d vector:

d, m = 3, 200
X = randn(d, m)   # 3-dimensional, 200-member ensemble
y = [0.1, -0.2, 0.5]

es(X, y)    # energy score

0.5125407162020792

vs(X, y)    # variogram score (default order p = 0.5)

0.5810756128962667

mmds(X, y)  # MMD score (Gaussian kernel, σ = 1)

-0.2519133002100307

See Scoring rules reference for the weighted variants of these scores.