frazane / scoringrules

Scoring rules for probabilistic forecast evaluation.
https://frazane.github.io/scoringrules/
Apache License 2.0
27 stars 4 forks source link

Add k-band and iid energy score estimator #25

Open simon-hirsch opened 1 month ago

simon-hirsch commented 1 month ago

If you have a lot of samples/scenarios/ensembles but little time (or little computational power):

Equations 15 and 16 in Berk & Ziel, 2019

sallen12 commented 1 month ago

The iid energy score estimator has been used in a few studies and would be a nice, straightforward addition to the package. I would avoid including the k-band estimator for now - there are some errors in Equation 16 of Berk & Ziel, 2019, and this approximation isn't used elsewhere. The fraction in Equation 15 should also be 2/M rather than M/2, see e.g. Moller et al., 2013, which is important to remember when we implement it.

simon-hirsch commented 1 month ago

All right @sallen12. - We've used the $k$ - band estimator e.g. here here. The fraction for the k-band should read 1 / (M * K) - There is even a typo in our paper, where a "-1" sneaked in. I agree though, that it's not the most used method to approximate the ES. I'll add the IID estimator only then.

sallen12 commented 1 month ago

Thanks for the link to the paper. It's probably just a notational misunderstanding then, but in your Eq 26, if K < M, then for any m > K, the lower bound on the second summation, k = m, will be higher than the upper bound, K. In this case, is it assumed that the summation is zero?

simon-hirsch commented 1 month ago

Hi, what we did, is essentially taking: $\hat{y}^{[m]} - \hat{y}^{[m+k]}$ for $K$, where $\hat{y}^{[m]}$ is the $m$-th ensemble member. I.e. for $K=1$, you just have $$1 / M \sum^M{m=1} | \hat{y}^{[m]} - \hat{y}^{[m+1]} |$$, for $K=2$ you get $$1 / (2M) (\sum^M{m=1} | \hat{y}^{[m]} - \hat{y}^{[m+1]}| + \sum^M_{m=1} | \hat{y}^{[m]} - \hat{y}^{[m+2]} |)$$ and so on. Does that make sense for you?

sallen12 commented 1 month ago

Thanks for clarifying. So, just to check I've understood, the second summation in Eq 26 starts from $k = 1$ rather than $k=m$? In this case everything makes sense to me, and I agree it could be a useful estimator to include in the package 👍

As an aside, for the implementation, it would be useful to randomise the ensemble members before applying this formula (also for the iid formula). Otherwise, if there is some sort of ordering among the ensemble members, e.g. because of how the ensembles are sampled from the predictive distribution, then only taking differences between nearby members (as is the case when $K=1$) will generally underestimate the underlying expectation.

simon-hirsch commented 1 month ago

Yeah, or it ends at $M+K$ if you want to start at $k=m$. Starting at 1 is less confusing though.

Essentially shuffle the ensembles once? Not sure how this will work with the gufuncs, but generally agree that it makes sense. Will check and also put a seed in the estimator 👍

sallen12 commented 1 month ago

Great, all clear now, thanks a lot!