`hesscrps_ct()` and `hesscrps_tt()` are too large by a factor of `scale`

[sbfnk]

Both return exactly scale times the true second derivative. At scale = 1 they are correct, which makes the error easy to miss. The normal and logistic counterparts are unaffected.

Reprex

The check uses a second difference of the CRPS itself, so it does not depend on the gradient functions. hesscrps_tt() in particular cannot be checked by differencing gradcrps_tt(), which has a separate problem (filed separately).

library(scoringRules)

d2 <- function(f, x, h = 1e-4) {
  (f(x + h) - 2 * f(x) + f(x - h)) / h^2
}

for (sc in c(0.5, 1, 2, 3)) {
  h <- hesscrps_ct(
    0.5,
    df = 3, location = 0, scale = sc, lower = -1, upper = 2
  )
  fd <- d2(
    function(m) {
      crps_ct(0.5, df = 3, location = m, scale = sc, lower = -1, upper = 2)
    },
    0
  )
  cat("ct  scale", sc, " hesscrps", h[1], " d2(crps)", fd, " ratio", h[1] / fd, "\n")
}

for (sc in c(0.5, 1, 2, 3)) {
  h <- hesscrps_tt(
    0.5,
    df = 3, location = 0, scale = sc, lower = -1, upper = 2
  )
  fd <- d2(
    function(m) {
      crps_tt(0.5, df = 3, location = m, scale = sc, lower = -1, upper = 2)
    },
    0
  )
  cat("tt  scale", sc, " hesscrps", h[1], " d2(crps)", fd, " ratio", h[1] / fd, "\n")
}

The ratio equals scale in every row:

ct  scale 0.5  hesscrps 0.4038342  d2(crps) 0.8076683  ratio 0.5
ct  scale 1  hesscrps 0.5361169  d2(crps) 0.5361169  ratio 1
ct  scale 2  hesscrps 0.4205469  d2(crps) 0.2102735  ratio 2
ct  scale 3  hesscrps 0.3075585  d2(crps) 0.1025196  ratio 2.999997
tt  scale 0.5  hesscrps 0.3018874  d2(crps) 0.6037747  ratio 0.5000001
tt  scale 1  hesscrps 0.2529334  d2(crps) 0.2529334  ratio 1
tt  scale 2  hesscrps 0.08822629  d2(crps) 0.04411293  ratio 2.00001
tt  scale 3  hesscrps 0.03453879  d2(crps) 0.01151389  ratio 2.99975

For contrast, the normal version has ratio 1 throughout:

for (sc in c(0.5, 2)) {
  h <- hesscrps_cnorm(0.5, location = 0, scale = sc, lower = -1, upper = 2)
  fd <- d2(
    function(m) {
      crps_cnorm(0.5, location = m, scale = sc, lower = -1, upper = 2)
    },
    0
  )
  cat("cnorm  scale", sc, " ratio", h[1] / fd, "\n")
}

cnorm  scale 0.5  ratio 1
cnorm  scale 2  ratio 1

Cause

The location-scale branch. hesscrps_cnorm ends it with a division:

hesscrps_cnorm((y - location)/scale, lower = lower, upper = upper)/scale

hesscrps_ct and hesscrps_tt have the same branch without the /scale:

hesscrps_ct((y - location)/scale, df, lower = lower, upper = upper)

The CRPS is scale-equivariant, CRPS(y; μ, σ) = σ · CRPS(z; 0, 1) with z = (y - μ)/σ, so the second derivative with respect to location picks up σ · (1/σ)^2 = 1/σ. The gradient needs no factor, which is why gradcrps_ct and gradcrps_tt are right to have none there.

Suggested fix

hesscrps_ct((y - location)/scale, df, lower = lower, upper = upper)/scale

and the same in hesscrps_tt. Every row above then matches:

ct  scale 0.5  fixed 0.8076683  d2(crps) 0.8076683
tt  scale 0.5  fixed 0.6037747  d2(crps) 0.6037747
ct  scale 1  fixed 0.5361169  d2(crps) 0.5361169
tt  scale 1  fixed 0.2529334  d2(crps) 0.2529334
ct  scale 2  fixed 0.2102735  d2(crps) 0.2102735
tt  scale 2  fixed 0.04411315  d2(crps) 0.04411293
ct  scale 3  fixed 0.1025195  d2(crps) 0.1025196
tt  scale 3  fixed 0.01151293  d2(crps) 0.01151389

The general branch of the same two functions has a separate set of problems, filed separately.