`hesscrps_ct()` and `hesscrps_tt()` return `NaN` for every input when `df = Inf`

[sbfnk]

The branch that handles infinite df, or bounds containing NA, never computes a value. It also returns a length-4 vector whatever the number of rows, and drops the dimnames. A third problem is hidden behind the first: where the branch does scale its result, it scales the wrong way.

gradcrps_ct() and gradcrps_tt() handle df = Inf correctly, so the gradient functions give the shape the hessian ones should have.

Reprex

library(scoringRules)

# 1. no value is ever computed
hesscrps_ct(0.5, df = Inf, location = 0, scale = 2, lower = -1, upper = 2)
hesscrps_tt(0.5, df = Inf, location = 0, scale = 2, lower = -1, upper = 2)

# 2. the same for the other trigger, an NA bound
hesscrps_ct(c(0.5, 0.7), df = 3, location = 0, scale = 2,
            lower = c(-1, NA), upper = 2)

# 3. three rows in, four numbers out, no dimnames
r <- hesscrps_ct(c(0.5, 0.7, -0.2), df = Inf, location = 0, scale = 2,
                 lower = -1, upper = 2)
length(r)
dim(r)
dimnames(r)

[1] NaN NaN NaN NaN
[1] NaN NaN NaN NaN
[1] NaN NaN NaN NaN
[1] 4
NULL
NULL

The gradient functions take the same arguments and work:

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

Since df = Inf is the normal case, hesscrps_cnorm gives the values the first call should have returned:

hesscrps_cnorm(0.5, location = 0, scale = 2, lower = -1, upper = 2)

         d2loc     d2scale dloc.dscale dscale.dloc
[1,] 0.2396528 -0.04137951   0.1125898   0.1125898

Cause

Three separate problems in the hesscrps_ct and hesscrps_tt else branch.

input has no df column, so nothing is ever selected. The data frame is built from z, scale, lower and upper:

input <- data.frame(z = y - location, scale = scale,
                    lower = lower - location, upper = upper - location)
isNaN <- is.na(input$z) | is.na(input$df) | input$df <= 1 |
  is.na(input$scale) | input$scale <= 0
ind2 <- input$df == Inf & !isNaN
ind3 <- !isNaN & !ind2

input$df is therefore NULL. is.na(NULL) and NULL <= 1 are both logical(0), and | with a zero-length operand gives logical(0), so isNaN is empty however many rows were passed. ind2 and ind3 are empty too, both any() guards are FALSE, and out is returned untouched:

input <- data.frame(z = c(0.5, 0.7), scale = 2, lower = -1, upper = 2)
length(is.na(input$z) | is.na(input$df) | input$df <= 1)

[1] 0

rep where matrix was meant. The allocation is

out <- rep(NaN, dim(input)[1L], 4,
           dimnames = list(NULL, c("d2loc", "d2scale",
                                   "dloc.dscale", "dscale.dloc")))

rep takes times and length.out in those positions, so the 4 binds to length.out and pins the result to four values regardless of dim(input)[1L], while dimnames goes into ... and is discarded. Once the selection is repaired the assignments would also need out[ind2, ] rather than out[ind2] to address rows of a matrix.

The scaling goes the wrong way. Both ind2 and ind3 compute scale * hesscrps_...(z/scale, ...), where the second derivative with respect to location calls for /scale. That leaves the branch out by scale^2. hesscrps_cnorm divides in the matching branch. This one is invisible until the missing column is fixed, since no row reaches the computation.

gradcrps_tt builds the same structure correctly, with df = df in the data frame, matrix(NaN, ...), and isNaN computed inside with(input, ...).

hesscrps_cnorm shares the rep/matrix problem. Its isNaN never mentions df, so it does compute values; it returns four of them for a two-row input, with a recycling warning.

Suggested fix

Follow gradcrps_tt: put df in the data frame, allocate with matrix, select rows with is.infinite(df), and divide by scale.

input <- data.frame(z = y - location, df = df, scale = scale,
                    lower = lower - location,
                    upper = upper - location)
out <- matrix(NaN, dim(input)[1L], 4,
              dimnames = list(NULL, c("d2loc", "d2scale",
                                      "dloc.dscale", "dscale.dloc")))
isNaN <- with(input, {
  is.na(z) | is.na(df) | df <= 1 |
    is.na(scale) | scale <= 0 |
    is.na(lower) | is.na(upper)
})
ind2 <- !isNaN & is.infinite(input$df)
ind3 <- !isNaN & !ind2
if (any(ind2)) {
  out[ind2, ] <- with(input[ind2, ],
                      hesscrps_cnorm(z/scale, lower = lower/scale,
                                     upper = upper/scale)/scale)
}
if (any(ind3)) {
  out[ind3, ] <- with(input[ind3, ],
                      hesscrps_ct(z/scale, df, lower = lower/scale,
                                  upper = upper/scale)/scale)
}
out

and the same in hesscrps_tt, with hesscrps_tnorm in place of hesscrps_cnorm.

The df = Inf values then reproduce second differences of crps_cnorm across scales:

scale 0.5  fixed d2loc 0.9629697  d2(crps_cnorm) 0.9629697
scale 1    fixed d2loc 0.6248942  d2(crps_cnorm) 0.6248942
scale 2    fixed d2loc 0.2396528  d2(crps_cnorm) 0.2396527
scale 3    fixed d2loc 0.1155737  d2(crps_cnorm) 0.1155737

The shape and dimnames come back:

         d2loc     d2scale dloc.dscale dscale.dloc
[1,] 0.2396528 -0.04137951   0.1125898   0.1125898
[2,] 0.2282251 -0.01957933   0.1472569   0.1472569
[3,] 0.2499373 -0.06157674  -0.0237725  -0.0237725

Rows with an NA bound stay NaN, and a mixed df vector routes each row to the right formula. With df = c(Inf, 3), scale = 2:

         d2loc     d2scale dloc.dscale dscale.dloc
[1,] 0.2396528 -0.04137951  0.11258979  0.11258979
[2,] 0.2102735 -0.04387819  0.09876162  0.09876162

The second row matches a second difference of crps_ct at df = 3, 0.2102735.

Note that the /scale here is the same correction the location-scale branch needs, see #68.