r - 如何计算 DRPS(离散秩概率分数)

Question

我正在复制论文 Forecasting the intermittent demand for slow-moving inventories: A modelling approach 中的评分规则
论文对评分规则的描述如下:



这是我的尝试

y <- rpois(n = 100, lambda = 10) # forecasted distribution
x <- 10 # actual value

drps_score <- function(x = value, y = q){
  # x = actual value (single observation); y = quantile forecasted value (vector)
  Fy = ecdf(y) # cdf function
  indicator <- ifelse(y - x > 0, 1, 0) # Heaviside
  score <- sum((indicator - Fy(y))^2)
  return(score)
}

> drps_score(x = x, y = y)
[1] 53.028

这似乎工作得很好，直到我提供如下的 0 向量:

y <- rep(x = 0, 100)
> drps_score(x = x, y = y)
[1] 0

我知道他们在本文中使用的一种方法是 0 预测，他们的结果并未显示 DRPS 为 0。这让我觉得计算是关闭的。

Answer 1

我认为这里有几个问题在起作用。
首先，我认为您没有在评分函数中计算正确的总和。该分数要求您对 y 的所有可能值(即所有正整数)而不是 y 的所有预测样本进行求和。
其次，我不认为上面的定义给出了想要的结果，当 y=x 时\hat F (y) 定义为 0 那么对于点质量为真实值的预测，你不会得到零分. (是的，我是说来源是“错误的”，或者至少有一个没有给出预期结果的定义。)这是一个重新制定的函数，我认为它可以解决这两个问题:

x <- 10 # actual value

drps_score <- function(x = value, y = q, nsum=100){
    # x = actual value (single observation); y = quantile forecasted value (vector)
    Fy = ecdf(y) # cdf function
    ysum <- 0:nsum
    indicator <- ifelse(ysum - x >= 0, 1, 0) # Heaviside
    score <- sum((indicator - Fy(ysum))^2)
    return(score)
}



> drps_score(x = x, y = rpois(n = 1000, lambda = 8))
[1] 1.248676
> drps_score(x = x, y = rpois(n = 1000, lambda = 9))
[1] 0.878183
> drps_score(x = x, y = rpois(n = 1000, lambda = 10))
[1] 0.692667
> drps_score(x = x, y = rep(10, 100))
[1] 0
> drps_score(x = x, y = rpois(n = 1000, lambda = 11))
[1] 0.883333

上面显示以真值 (lambda=10) 为中心的分布对于不是点质量的分布的得分最低。