Feature request: Skill Scores

[ahuang11]

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]

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]

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]

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/

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

[ahuang11]

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

[aaronspring]

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]

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]

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]

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]

See #284. Leaving open as good discussion here.

[matteodefelice]

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.