Calculating the RPS

\[ RPS(P,y) = \sum^\infty_{k=0} \left \{ P(Y\leq k) - 1(y\leq k) \right \}^2\]

rps <- function(y, mu, maxk=1000) {
    k <- 1:maxk
    sum( (ppois(k, mu) - (y<=k))^2 )
}

maxy <- 50
scores <- rep(NA, maxy)
for(i in 1:maxy) {
    scores[i] <- rps(i, mu=20)
}
qplot(1:maxy, scores, geom="line")

Portions of the calculations when \(\mu = y\), i.e. predictive distribution has the right mean.

k <- 1:30
y <- 20
mu <- 20
cumsum((ppois(k, mu) - (y<=k))^2)
##  [1] 1.873524e-15 2.093674e-13 1.047319e-11 2.975975e-10 5.468479e-09
##  [6] 7.055597e-08 6.767585e-07 5.033409e-06 2.998755e-05 1.468808e-04
## [11] 6.042770e-04 2.126213e-03 6.499077e-03 1.749559e-02 4.199196e-02
## [16] 9.086576e-02 1.790916e-01 3.245743e-01 5.457162e-01 7.401156e-01
## [21] 8.670669e-01 9.451250e-01 9.902843e-01 1.014862e+00 1.027447e+00
## [26] 1.033514e+00 1.036268e+00 1.037447e+00 1.037923e+00 1.038104e+00

Portions of the calculations when \(\mu \neq y\) is wrong.

k <- 1:30
y <- 20
mu <- 10
cumsum((ppois(k, mu) - (y<=k))^2)
##  [1] 2.493996e-07 7.918952e-06 1.147529e-04 9.704727e-04 5.470999e-03
##  [6] 2.240779e-02 7.090492e-02 1.816739e-01 3.913735e-01 7.313088e-01
## [11] 1.216806e+00 1.843367e+00 2.590666e+00 3.430715e+00 4.335609e+00
## [16] 5.282257e+00 6.253906e+00 7.239585e+00 8.232688e+00 8.232690e+00
## [21] 8.232691e+00 8.232691e+00 8.232691e+00 8.232691e+00 8.232691e+00
## [26] 8.232691e+00 8.232691e+00 8.232691e+00 8.232691e+00 8.232691e+00