Feature request: Skill Scores

ahuang11 commented 4 years ago

https://confluence.ecmwf.int/display/FUG/12.B+Statistical+Concepts+-+Probabilistic+Data I think it's just:

def brier_skill_score(obs, fct, ref, threshold):
    bscore_fct = (xs.brier_score(obs > threshold, (fct > threshold).mean('member')) 
    bscore_ref = (xs.brier_score(obs > threshold, (ref > threshold).mean('member')) 
    bskill = 1 - (bscore_fct / bscore_ref)
    return bskill

where ref can be persistence or climatology (where each year of the climatology is a member)

raybellwaves commented 4 years ago

You can see how @matteodefelice did it in his notebook (https://colab.research.google.com/drive/1wWHz_SMCHNuos5fxWRUJTcB6wqkTJQCR)

brier_score = xs.brier_score(obs > obs.quantile(2/3, dim='time'),
                            (fct > fct.quantile(2/3, dim='time')).mean("member"))
baseline = xs.brier_score(obs > obs.quantile(2/3, dim='time'), 2/3)
bss = 1 - (brier_score / baseline)

ahuang11 commented 4 years ago

Yeah, I think both ways are valid, just a different baseline; maybe his can be used as the default if ref=None

def brier_skill_score(obs, fct, threshold, ref=None):
    bscore_fct = (xs.brier_score(obs > threshold, (fct > threshold).mean('member'))
    if ref is None:
        bscore_ref = xs.brier_score(obs > obs.quantile(2/3, dim='time'), 2/3)
    else:
        bscore_ref = (xs.brier_score(obs > threshold, (ref > threshold).mean('member')) 
    bskill = 1 - (bscore_fct / bscore_ref)
    return bskill

aaronspring commented 4 years ago

I think we should somehow come up with a pattern/wrapper that translates every skill into a skill score.

Brier Score and Brier Skill Score
RMSE and RMSE Skill Score
...

The formula would be from https://www.cawcr.gov.au/projects/verification/

Image 16 08 20 at 12 34

anyone ideas how to do this with few lines of code?

ahuang11 commented 4 years ago

Is it applicable across all scores though? I am only familiar with brier skill score (first time heard of RMSE skill score)

aaronspring commented 4 years ago

I think so. In climpred we have a metrics class with attributes min max perfect. Maybe this would be nice to move to xskillscore

bradyrx commented 4 years ago

I think this sort of thing should be at climpred. a and b in xskillscore are two time series for instance. Let's say a is the forecast and b the verification product. To get BSS, you'd evaluate a against b, and then a persistence derivation of b against b for instance. Since we manage reference forecasts at climpred I'd think over there we would just call xskillscore BS twice.

bradyrx commented 4 years ago

I guess with the addition of the sign test, this could still fit here. If you hand xskillscore two different forecasts (one being the dynamical and the other the reference) you can compute skill scores. climpred would then wrap that to do nice looping through all the leads.

aaronspring commented 4 years ago

this is could be a template: https://github.com/bradyrx/climpred/blob/8a3fc954df2043f998987a1964059a9dc0d2e11c/climpred/metrics.py#L128-L140 were we would need the perfect attribute.

But many we should overdo things. xs is perfect for doing easy functions xs.metric(a,b,dim). skill scores can actually really easy be calculated from this manually. I dont see a large computational benefit of implementing this compared to having this sequentially manual_scoring_function(xs.metric(a,b,dim), reference(*args,**kwargs)).

We at least shouldnt change the existing API when implementing this.

raybellwaves commented 3 years ago

See https://github.com/xarray-contrib/xskillscore/issues/284. Leaving open as good discussion here.

matteodefelice commented 3 years ago

By the way, it seems that the line:

baseline = xs.brier_score(obs_final >  obs_final.quantile(2/3, dim = 'time'), 2/3)

doesn't work any more, I get this error on my Colab:

AttributeError                            Traceback (most recent call last)

<ipython-input-68-ad9dcdfd6ce1> in <module>()
----> 1 baseline = xs.brier_score(obs_final >  obs_final.quantile(2/3, dim = 'time'), 2/3)

/usr/local/lib/python3.8/site-packages/xskillscore/core/probabilistic.py in brier_score(observations, forecasts, member_dim, fair, dim, weights, keep_attrs)
    348                 res = (e / M - o) ** 2 - e * (M - e) / (M ** 2 * (M - 1))
    349         else:
--> 350             if member_dim in forecasts.dims:
    351                 forecasts = forecasts.mean(member_dim)
    352             res = xr.apply_ufunc(

AttributeError: 'float' object has no attribute 'dims'

EDIT: yes, now it requires an xarray object and it doesn't work any more with numbers.