Differences from R

ScoringRules.jl is a port of the R package scoringRules (version 1.1.3). In most cases the two packages produce identical numerical results, but there are known divergences that users coming from R should be aware of.

CRPS gradients and Hessians

R exports closed-form CRPS derivatives with respect to location and scale (gradcrps_*, hesscrps_*) for the normal, logistic and Student's families and their truncated and censored variants. ScoringRules.jl provides no closed-form derivatives; gradients come from automatic differentiation of crps.

GEV CRPS: Gumbel case (shape ≈ 0)

For a GeneralizedExtremeValue(μ, σ, 0) distribution (the Gumbel limit,  ), the general closed-form expression is numerically unstable. R's implementation falls back to numerical integration in this limit.

ScoringRules.jl instead evaluates the Gumbel limit analytically using the exponential integral from SpecialFunctions.jl:

where    and    is the Euler–Mascheroni constant. This is the same formula as for the general GEV but taken in the limit  . Numerically it is more accurate than the integrate fallback because expinti is implemented in arbitrary precision. Results agree with R's numerical integration to within the integration tolerance in all tested cases.

Ensemble DSS: population vs sample variance

dss(dat, y) for an ensemble dat uses the population variance   , computed as mean(dat.^2) - mean(dat).^2. This matches R's dss_sample / dss_edf, which also uses the biased (population) estimator. Julia's var function uses the   denominator, so passing var(dat) directly would not match R. The implementation avoids Statistics.var deliberately.

Ensemble member weights

The sample crps (method = :edf) and dss accept member weights w, matching R's crps_sample and dss_sample. Member weights must be finite, non-negative and sum to a positive value, otherwise an error is thrown. R also errors on missing, infinite and negative weights but returns NaN for an all-zero weight vector.

Parameter conventions

The following parameter conventions differ from R's function arguments:

Quantity	R convention	Julia / Distributions.jl
LogLogistic	locationlog, scalelog	LogLogistic(α, β): α = exp(locationlog), β = 1/scalelog
GeneralizedPareto	location, scale, shape	GeneralizedPareto(μ, σ, ξ) — same order
NegativeBinomial	size, prob (success prob)	NegativeBinomial(r, p) — same meaning
Hypergeometric	m, n, k (white, black, draws)	Hypergeometric(s, f, n) (successes, failures, draws) — s=m, f=n, n=k
TwoPieceNormal	location, scale1, scale2	TwoPieceNormal(location, scale1, scale2) — same

Reporting discrepancies

If you find a numerical result that differs from R and is not listed above, please open an issue on the GitHub repository. The test suite includes a set of reference values generated from R's scoringRules package, so discrepancies can usually be identified precisely.