`gradcrps_tt()` disagrees with finite differences of `crps_tt()`

[sbfnk]

gradcrps_tt() returns location and scale derivatives that do not match central differences of crps_tt(). The errors include sign flips and factors of 100, and occur at scale = 1 too.

gradcrps_norm(), gradcrps_logis(), gradcrps_t() and the censored variants agree with finite differences to around 1e-9 on the same grid.

Reprex

library(scoringRules)

fd <- function(f, x, h = 1e-5) {
  (f(x + h) - f(x - h)) / (2 * h)
}

cases <- list(
  list(y = 0.5, df = 3, location = -1, scale = 0.5, lower = 0, upper = Inf),
  list(y = 1.0, df = 5, location = 0, scale = 1.0, lower = -1, upper = 2),
  list(y = -0.5, df = 4, location = 0, scale = 2.0, lower = -2, upper = Inf)
)

for (p in cases) {
  g <- gradcrps_tt(
    p$y,
    df = p$df, location = p$location, scale = p$scale,
    lower = p$lower, upper = p$upper
  )
  dloc <- fd(
    function(m) {
      crps_tt(
        p$y,
        df = p$df, location = m, scale = p$scale,
        lower = p$lower, upper = p$upper
      )
    },
    p$location
  )
  dscale <- fd(
    function(s) {
      crps_tt(
        p$y,
        df = p$df, location = p$location, scale = s,
        lower = p$lower, upper = p$upper
      )
    },
    p$scale
  )
  cat(
    "y", p$y, " df", p$df, " location", p$location,
    " scale", p$scale, " lower", p$lower, " upper", p$upper, "\n"
  )
  cat("  gradcrps_tt  dloc", g[1], " dscale", g[2], "\n")
  cat("  finite diff  dloc", dloc, " dscale", dscale, "\n")
}

y 0.5  df 3  location -1  scale 0.5  lower 0  upper Inf
  gradcrps_tt  dloc 5.447827  dscale 10.99537
  finite diff  dloc -0.03045686  dscale 0.03879918
y 1  df 5  location 0  scale 1  lower -1  upper 2
  gradcrps_tt  dloc -0.3097082  dscale -0.07262911
  finite diff  dloc -0.3906397  dscale -0.234492
y -0.5  df 4  location 0  scale 2  lower -2  upper Inf
  gradcrps_tt  dloc 0.1632217  dscale 0.5088614
  finite diff  dloc 0.3330361  dscale 0.339047

In the first case gradcrps_tt() reports dloc = +5.45 where the function it differentiates has slope -0.030.

Over a 252-row grid in y, df, location, scale, lower and upper, 102 rows disagree. The agreeing rows are those where y sits well clear of a finite truncation bound.

Cause

The sign of the G terms in the standardised branch. gradcrps_tnorm and gradcrps_tt are otherwise the same code, with G in the t version playing the role d plays in the normal version:

# gradcrps_tnorm
term2 <- 2 * d_l/a * ((z * p_l - z * p_z + d_l - d_z)/a + term0)
term3 <- 2 * d_u/a * ((z * p_u - z * p_z + d_u - d_z)/a + term0)

# gradcrps_tt
term2 <- 2 * d_l/a * ((z * p_l - z * p_z - G_l + G_z)/a + term0)
term3 <- 2 * d_u/a * ((z * p_u - z * p_z - G_u + G_z)/a + term0)

G(z) = (df + z^2)/(df - 1) * dt(z, df) is the Student-t counterpart of dnorm: -G is the antiderivative of z * dt(z, df) exactly as -dnorm is of z * dnorm(z), and G(z) tends to dnorm(z) as df grows. So G belongs in those expressions with the signs d carries.

hesscrps_tt defines the same quantity with the opposite sign, G <- -(df + z^2)/(df - 1) * PDF, which may be where the discrepancy came from.

Suggested fix

term2 <- 2 * d_l/a * ((z * p_l - z * p_z + G_l - G_z)/a + term0)
term3 <- 2 * d_u/a * ((z * p_u - z * p_z + G_u - G_z)/a + term0)

The three cases then reproduce the finite differences to seven significant digits:

y 0.5  df 3  location -1  scale 0.5  lower 0  upper Inf
  proposed     dloc -0.03045686  dscale 0.03879918
  finite diff  dloc -0.03045686  dscale 0.03879918
y 1  df 5  location 0  scale 1  lower -1  upper 2
  proposed     dloc -0.3906397  dscale -0.234492
  finite diff  dloc -0.3906397  dscale -0.234492
y -0.5  df 4  location 0  scale 2  lower -2  upper Inf
  proposed     dloc 0.3330361  dscale 0.339047
  finite diff  dloc 0.3330361  dscale 0.339047

Over a 540-row grid in y, df, location, scale and three bound configurations, the rows disagreeing with finite differences go from 397 to 0, at a step size large enough that quadrature noise in crps_tt does not dominate the difference.