`crps_ct()` and `crps_tt()` return `NaN` for `df = Inf`

[sbfnk]

Both return NaN for infinite df, where the limiting normal value is wanted. The code to compute that value is present but sits in a branch ordinary input never reaches.

gradcrps_ct() and gradcrps_tt() accept df = Inf and agree with their normal counterparts to full precision, so the gradient functions show the behaviour expected here.

Reprex

library(scoringRules)

for (sc in c(0.5, 1, 2, 3)) {
  cat(
    "scale", sc,
    " crps_ct", crps_ct(0.5, df = Inf, location = 0, scale = sc,
                        lower = -1, upper = 2),
    " crps_cnorm", crps_cnorm(0.5, location = 0, scale = sc,
                              lower = -1, upper = 2),
    "\n"
  )
}

scale 0.5  crps_ct NaN  crps_cnorm 0.3011697
scale 1  crps_ct NaN  crps_cnorm 0.3240666
scale 2  crps_ct NaN  crps_cnorm 0.4337524
scale 3  crps_ct NaN  crps_cnorm 0.509778

crps_tt() behaves the same way against crps_tnorm():

scale 0.5  crps_tt NaN  crps_tnorm 0.2905801
scale 1  crps_tt NaN  crps_tnorm 0.237287
scale 2  crps_tt NaN  crps_tnorm 0.2382557
scale 3  crps_tt NaN  crps_tnorm 0.2436877

Large finite df converges to the value the Inf call should return, so the target is not in doubt:

for (d in c(1e3, 1e5, 1e7)) {
  cat("df", d, ":", crps_ct(0.5, df = d, location = 0, scale = 2,
                            lower = -1, upper = 2), "\n")
}

df 1000 : 0.4338173
df 1e+05 : 0.433753
df 1e+07 : 0.4337524

The gradient functions accept the same arguments:

gradcrps_ct(0.5, df = Inf, location = 0, scale = 2, lower = -1, upper = 2)
gradcrps_cnorm(0.5, location = 0, scale = 2, lower = -1, upper = 2)

           dloc     dscale
[1,] -0.1273887 0.09475383
           dloc     dscale
[1,] -0.1273887 0.09475383

Cause

crps_ct does handle df == Inf, in the inner else of its location-scale branch:

if (all(scale > 0, na.rm = TRUE)) {
  scale * crps_ct(y / scale, df, lower = lower, upper = upper)
} else {
  out <- scale * crps_ct(y / scale, df, lower = lower, upper = upper)
  ind1 <- df == Inf
  ...
  out[ind1] <- rep_len(scale * crps_cnorm(y / scale, lower = lower,
                                          upper = upper), length(out))[ind1]

That else runs only when some scale is zero or negative. With every scale positive, which is the ordinary case, control takes the branch above it and ind1 is never evaluated. Supplying a non-positive scale alongside a positive one reaches the code and shows it is correct:

crps_ct(c(0.5, 0.5), df = Inf, location = 0, scale = c(2, 0),
        lower = -1, upper = 2)

[1] 0.4337524 0.5000000

The first element is now the crps_cnorm value.

At scale = 1 the function does not reach that branch at all. The dispatch tests only the scale:

if (identical(scale, 1)) {

so df = Inf enters the standardised branch, where G_z evaluates (df + z^2)/(df - 1) as Inf/Inf. bfrac is indeterminate twice over, 2*sqrt(df)/(df - 1) giving Inf/Inf and beta(0.5, df - 0.5)/beta(0.5, 0.5*df)^2 giving 0/0:

crps_ct(0.5, df = Inf, lower = -1, upper = 2)

[1] NaN

gradcrps_ct guards its standardised branch on the degrees of freedom as well as the scale:

all_df_in_1_to_Inf <- all(is.finite(df) & df > 1)

Infinite df therefore falls through to a general branch that dispatches on is.infinite(df) element-wise. crps_ct and crps_tt have no equivalent guard.

Suggested fix

Give crps_ct the three-way dispatch gradcrps_ct already has: standardised for finite df at unit scale, location-scale for finite df, and a general branch that handles infinite df element-wise. Routing df = Inf through the general branch also avoids the infinite recursion that adding the guard alone would cause, since the location-scale branch recurses with df unchanged.

finite_df <- all(is.finite(df) & df > 1)
if (identical(scale, 1) && finite_df) {
  # standardised branch, unchanged
} else if (finite_df && all(is.finite(scale) & scale > 0)) {
  if (!identical(lower, -Inf)) lower <- lower / scale
  if (!identical(upper,  Inf)) upper <- upper / scale
  scale * crps_ct(y / scale, df, lower = lower, upper = upper)
} else {
  input <- data.frame(z = y, df = df, scale = scale,
                      lower = lower, upper = upper)
  out <- rep(NaN, dim(input)[1L])
  isNaN <- with(input, is.na(z) | is.na(df) | df <= 1 |
                  is.na(scale) | scale < 0)
  ind_zero <- !isNaN & input$scale == 0 & input$lower <= input$upper
  ind_inf  <- !isNaN & !ind_zero & is.infinite(input$df)
  ind_fin  <- !isNaN & !ind_zero & !ind_inf
  if (any(ind_zero)) {
    out[ind_zero] <- with(input[ind_zero, ],
                          abs(z - pmax(lower, 0) - pmin(upper, 0)))
  }
  if (any(ind_inf)) {
    out[ind_inf] <- with(input[ind_inf, ],
                         crps_cnorm(z, scale = scale,
                                    lower = lower, upper = upper))
  }
  if (any(ind_fin)) {
    out[ind_fin] <- with(input[ind_fin, ],
                         scale * crps_ct(z / scale, df,
                                         lower = lower / scale,
                                         upper = upper / scale))
  }
  out
}

crps_tt takes the same shape, with crps_tnorm in place of crps_cnorm and its own zero-scale rule (scale == 0 & lower < 0 & upper > 0, giving abs(z)).

Every df = Inf case then matches the normal version:

scale 0.5  current NaN  fixed 0.3011697  crps_cnorm 0.3011697
scale 1  current NaN  fixed 0.3240666  crps_cnorm 0.3240666
scale 2  current NaN  fixed 0.4337524  crps_cnorm 0.4337524
scale 3  current NaN  fixed 0.509778  crps_cnorm 0.509778

scale 0.5  current NaN  fixed 0.2905801  crps_tnorm 0.2905801
scale 1  current NaN  fixed 0.237287  crps_tnorm 0.237287
scale 2  current NaN  fixed 0.2382557  crps_tnorm 0.2382557
scale 3  current NaN  fixed 0.2436877  crps_tnorm 0.2436877

Finite df is untouched. Over a 216-row grid in y, df, scale, location and two bound configurations, the maximum absolute difference against the current implementation is 0 for both functions, with no change in which entries are NaN.

A mixed df vector routes each element to the right formula, and a zero scale still takes precedence over infinite df. With df = c(Inf, 3, Inf) and scale = c(2, 2, 0):

[1] 0.4337524 0.4548508 0.5000000

matching crps_cnorm, crps_ct at df = 3, and the point-mass value in turn.