search

FK83 / scoringRules

scoring rules to evaluate probabilistic forecasts
53 stars 16 forks source link

fix `crps_nbinom` #51

Closed aijordan closed 2 years ago

aijordan commented 2 years ago

This PR changes the automatic dispatch of hypergeo::hypergeo to relevant subroutines that extend the range of size in crps_nbinom that returns accurate results.

Fixes #50

Before merging, we should probably do some additional testing.

aijordan commented 2 years ago

At least the original example from #50 seems to work now:

dat <- list(
  observed = c(
    65L, 62L, 48L, 47L, 60L, 75L, 62L, 50L, 49L, 
    54L, 62L, 51L, 57L, 56L, 48L, 47L, 43L, 35L, 53L, 34L, 58L, 38L, 
    56L, 42L, 58L, 55L, 51L, 54L, 64L, 51L, 45L, 59L, 45L, 46L, 67L, 
    44L, 51L, 47L, 50L, 56L, 69L, 56L, 35L, 39L, 50L, 45L, 43L, 51L, 
    47L, 53L, 36L, 30L
  ), 
  pred_means = c(
    48.3099006480877, 48.1856853530793, 
    48.0249249252662, 47.8581559836325, 47.7155235009397, 47.6268375218476, 
    47.6218179266744, 47.724978436259, 47.9185569151005, 48.1648674172695, 
    48.4256057306771, 48.6616929438018, 48.833317728955, 48.9007310037522, 
    48.8620539654412, 48.7821309824475, 48.7324016131518, 48.7840863337008, 
    49.0091496602714, 49.4823260083412, 50.2675286039073, 51.3049811917389, 
    52.4702988591609, 53.6283020516166, 54.6318059743995, 55.3238719063732, 
    55.5453815997929, 55.2375969406439, 54.5183299747122, 53.5268537918377, 
    52.3989708593686, 51.2620683576382, 50.2333662682769, 49.4087504392848, 
    48.7853128755403, 48.3188875640644, 47.9682808060746, 47.6944186604847, 
    47.4595844049793, 47.2271670054595, 46.9786501957888, 46.7263257219719, 
    46.4853115157984, 46.2704191049954, 46.0962135356727, 45.9771198543237, 
    45.9245710565769, 45.9269448553235, 45.9619878567896, 46.0074298534593, 
    46.0409403503701, 46.0401323330852
  ), 
  size = 249.137522940415
)

scoringRules::crps(
  dat$observed,
  family = "negative-binomial",
  mu = dat$pred_means,
  size = dat$size
)
#>  [1] 12.524908  9.811189  1.764168  1.779142  8.441514 23.132536 10.355357
#>  [8]  2.083245  1.850250  3.582602  9.582249  2.116417  5.081308  4.352307
#> [15]  1.800542  1.903845  3.316265  9.631608  2.681845 11.223429  4.766892
#> [22]  9.090524  2.548535  7.520443  2.523148  1.915409  2.806529  1.961657
#> [29]  5.981568  2.137097  4.305351  4.764429  3.066962  2.312115 13.983689
#> [36]  2.635610  2.310317  1.766204  2.150713  5.545419 17.828083  5.938077
#> [43]  7.604773  4.256486  2.611230  1.748063  2.109028  3.158801  1.803828
#> [50]  4.294543  6.368936 11.923287
aijordan commented 2 years ago

I am as confident as I can be that crps_nbinom now returns accurate values for almost any parameter combination. The following code compares the cprs_nbinom implementation against the crps implementation of a finite support approximation.

size <- rep(each = 9, c(
  10^(0:5), # range of integer values
  2 * 10^(-2:5) + 10^(pmin(-3:4, -1)))) # range of noninteger values
prob <- rep(length.out = length(size), c(
  0.0001, 0.1, 0.5, 0.82, 0.83, 0.90, 0.95, 0.99, 0.9999))
obs <- qnbinom(0.5, size, prob = prob)
dat <- purrr::pmap(
  list(size, prob),
  \(s, p) qnbinom(c(1e-10, 1 - 1e-10), s, prob = p) |> (\(x) x[1]:x[2])())
df <- tibble::tibble(
  size = size,
  prob = prob,
  obs = obs,
  crps_edf = purrr::pmap_dbl(
    list(obs, dat, size, prob),
    \(y, dat, s, p) scoringRules::crps_sample(y, dat, w = dnbinom(dat, s, prob = p))),
  crps_nbinom = scoringRules::crps_nbinom(obs, size, prob = prob))
all.equal(df$crps_edf, df$crps_nbinom)
#> [1] TRUE

For reference:

print(df, n = nrow(df))
#> # A tibble: 126 x 5
#>           size   prob        obs crps_edf crps_nbinom
#>          <dbl>  <dbl>      <dbl>    <dbl>       <dbl>
#>   1      1     0.0001       6931 1.93e+ 3    1.93e+ 3
#>   2      1     0.1             6 1.83e+ 0    1.83e+ 0
#>   3      1     0.5             0 3.33e- 1    3.33e- 1
#>   4      1     0.82            0 3.35e- 2    3.35e- 2
#>   5      1     0.83            0 2.98e- 2    2.98e- 2
#>   6      1     0.9             0 1.01e- 2    1.01e- 2
#>   7      1     0.95            0 2.51e- 3    2.51e- 3
#>   8      1     0.99            0 1.00e- 4    1.00e- 4
#>   9      1     1.00            0 1.00e- 8    1.00e- 8
#>  10     10     0.0001      96677 7.26e+ 3    7.26e+ 3
#>  11     10     0.1            87 6.89e+ 0    6.89e+ 0
#>  12     10     0.5             9 1.04e+ 0    1.04e+ 0
#>  13     10     0.82            2 3.52e- 1    3.52e- 1
#>  14     10     0.83            2 3.46e- 1    3.46e- 1
#>  15     10     0.9             1 2.27e- 1    2.27e- 1
#>  16     10     0.95            0 1.72e- 1    1.72e- 1
#>  17     10     0.99            0 9.17e- 3    9.17e- 3
#>  18     10     1.00            0 9.99e- 7    9.99e- 7
#>  19    100     0.0001     996569 2.33e+ 4    2.33e+ 4
#>  20    100     0.1           897 2.21e+ 1    2.21e+ 1
#>  21    100     0.5            99 3.30e+ 0    3.30e+ 0
#>  22    100     0.82           22 1.21e+ 0    1.21e+ 0
#>  23    100     0.83           20 1.16e+ 0    1.16e+ 0
#>  24    100     0.9            11 8.09e- 1    8.09e- 1
#>  25    100     0.95            5 5.32e- 1    5.32e- 1
#>  26    100     0.99            1 2.13e- 1    2.13e- 1
#>  27    100     1.00            0 9.90e- 5    9.90e- 5
#>  28   1000     0.0001    9995667 7.39e+ 4    7.39e+ 4
#>  29   1000     0.1          8997 7.01e+ 1    7.01e+ 1
#>  30   1000     0.5           999 1.05e+ 1    1.05e+ 1
#>  31   1000     0.82          219 3.82e+ 0    3.82e+ 0
#>  32   1000     0.83          205 3.67e+ 0    3.67e+ 0
#>  33   1000     0.9           111 2.59e+ 0    2.59e+ 0
#>  34   1000     0.95           52 1.74e+ 0    1.74e+ 0
#>  35   1000     0.99           10 7.34e- 1    7.34e- 1
#>  36   1000     1.00            0 9.08e- 3    9.08e- 3
#>  37  10000     0.0001   99986667 2.34e+ 5    2.34e+ 5
#>  38  10000     0.1         89997 2.22e+ 2    2.22e+ 2
#>  39  10000     0.5          9999 3.30e+ 1    3.30e+ 1
#>  40  10000     0.82         2195 1.21e+ 1    1.21e+ 1
#>  41  10000     0.83         2048 1.16e+ 1    1.16e+ 1
#>  42  10000     0.9          1111 8.21e+ 0    8.21e+ 0
#>  43  10000     0.95          526 5.50e+ 0    5.50e+ 0
#>  44  10000     0.99          101 2.36e+ 0    2.36e+ 0
#>  45  10000     1.00            1 2.12e- 1    2.12e- 1
#>  46 100000     0.0001  999896667 7.39e+ 5    7.39e+ 5
#>  47 100000     0.1        899997 7.01e+ 2    7.01e+ 2
#>  48 100000     0.5         99999 1.05e+ 2    1.05e+ 2
#>  49 100000     0.82        21951 3.82e+ 1    3.82e+ 1
#>  50 100000     0.83        20482 3.67e+ 1    3.67e+ 1
#>  51 100000     0.9         11111 2.60e+ 1    2.60e+ 1
#>  52 100000     0.95         5263 1.74e+ 1    1.74e+ 1
#>  53 100000     0.99         1010 7.46e+ 0    7.46e+ 0
#>  54 100000     1.00           10 7.29e- 1    7.29e- 1
#>  55      0.021 0.0001          0 5.88e+ 0    5.88e+ 0
#>  56      0.021 0.1             0 4.25e- 3    4.25e- 3
#>  57      0.021 0.5             0 2.28e- 4    2.28e- 4
#>  58      0.021 0.82            0 1.75e- 5    1.75e- 5
#>  59      0.021 0.83            0 1.54e- 5    1.54e- 5
#>  60      0.021 0.9             0 4.90e- 6    4.90e- 6
#>  61      0.021 0.95            0 1.16e- 6    1.16e- 6
#>  62      0.021 0.99            0 4.45e- 8    4.45e- 8
#>  63      0.021 1.00            0 4.41e-12    4.41e-12
#>  64      0.21  0.0001        247 3.98e+ 2    3.98e+ 2
#>  65      0.21  0.1             0 3.44e- 1    3.44e- 1
#>  66      0.21  0.5             0 2.07e- 2    2.07e- 2
#>  67      0.21  0.82            0 1.69e- 3    1.69e- 3
#>  68      0.21  0.83            0 1.49e- 3    1.49e- 3
#>  69      0.21  0.9             0 4.81e- 4    4.81e- 4
#>  70      0.21  0.95            0 1.15e- 4    1.15e- 4
#>  71      0.21  0.99            0 4.45e- 6    4.45e- 6
#>  72      0.21  1.00            0 4.41e-10    4.41e-10
#>  73      2.1   0.0001      17775 3.11e+ 3    3.11e+ 3
#>  74      2.1   0.1            16 2.94e+ 0    2.94e+ 0
#>  75      2.1   0.5             2 4.56e- 1    4.56e- 1
#>  76      2.1   0.82            0 1.25e- 1    1.25e- 1
#>  77      2.1   0.83            0 1.12e- 1    1.12e- 1
#>  78      2.1   0.9             0 4.03e- 2    4.03e- 2
#>  79      2.1   0.95            0 1.05e- 2    1.05e- 2
#>  80      2.1   0.99            0 4.36e- 4    4.36e- 4
#>  81      2.1   1.00            0 4.41e- 8    4.41e- 8
#>  82     20.1   0.0001     197657 1.04e+ 4    1.04e+ 4
#>  83     20.1   0.1           178 9.85e+ 0    9.85e+ 0
#>  84     20.1   0.5            20 1.47e+ 0    1.47e+ 0
#>  85     20.1   0.82            4 5.24e- 1    5.24e- 1
#>  86     20.1   0.83            4 5.00e- 1    5.00e- 1
#>  87     20.1   0.9             2 3.39e- 1    3.39e- 1
#>  88     20.1   0.95            1 2.18e- 1    2.18e- 1
#>  89     20.1   0.99            0 3.38e- 2    3.38e- 2
#>  90     20.1   1.00            0 4.03e- 6    4.03e- 6
#>  91    200.    0.0001    1997468 3.30e+ 4    3.30e+ 4
#>  92    200.    0.1          1798 3.13e+ 1    3.13e+ 1
#>  93    200.    0.5           200 4.67e+ 0    4.67e+ 0
#>  94    200.    0.82           44 1.71e+ 0    1.71e+ 0
#>  95    200.    0.83           41 1.64e+ 0    1.64e+ 0
#>  96    200.    0.9            22 1.15e+ 0    1.15e+ 0
#>  97    200.    0.95           10 7.79e- 1    7.79e- 1
#>  98    200.    0.99            2 3.13e- 1    3.13e- 1
#>  99    200.    1.00            0 3.93e- 4    3.93e- 4
#> 100   2000.    0.0001   19995667 1.04e+ 5    1.04e+ 5
#> 101   2000.    0.1         17998 9.91e+ 1    9.91e+ 1
#> 102   2000.    0.5          2000 1.48e+ 1    1.48e+ 1
#> 103   2000.    0.82          439 5.41e+ 0    5.41e+ 0
#> 104   2000.    0.83          409 5.19e+ 0    5.19e+ 0
#> 105   2000.    0.9           222 3.67e+ 0    3.67e+ 0
#> 106   2000.    0.95          105 2.46e+ 0    2.46e+ 0
#> 107   2000.    0.99           20 1.05e+ 0    1.05e+ 0
#> 108   2000.    1.00            0 3.32e- 2    3.32e- 2
#> 109  20000.    0.0001  199977667 3.30e+ 5    3.30e+ 5
#> 110  20000.    0.1        179998 3.14e+ 2    3.14e+ 2
#> 111  20000.    0.5         20000 4.67e+ 1    4.67e+ 1
#> 112  20000.    0.82         4390 1.71e+ 1    1.71e+ 1
#> 113  20000.    0.83         4096 1.64e+ 1    1.64e+ 1
#> 114  20000.    0.9          2222 1.16e+ 1    1.16e+ 1
#> 115  20000.    0.95         1052 7.78e+ 0    7.78e+ 0
#> 116  20000.    0.99          202 3.34e+ 0    3.34e+ 0
#> 117  20000.    1.00            2 3.11e- 1    3.11e- 1
#> 118 200000.    0.0001 1999797667 1.05e+ 6    1.05e+ 6
#> 119 200000.    0.1       1799998 9.91e+ 2    9.91e+ 2
#> 120 200000.    0.5        200000 1.48e+ 2    1.48e+ 2
#> 121 200000.    0.82        43902 5.41e+ 1    5.41e+ 1
#> 122 200000.    0.83        40964 5.19e+ 1    5.19e+ 1
#> 123 200000.    0.9         22222 3.67e+ 1    3.67e+ 1
#> 124 200000.    0.95        10526 2.46e+ 1    2.46e+ 1
#> 125 200000.    0.99         2020 1.06e+ 1    1.06e+ 1
#> 126 200000.    1.00           20 1.04e+ 0    1.04e+ 0

Created on 2022-01-12 by the reprex package (v2.0.1)

mbojan commented 2 years ago

Splendid!