Along the lines of https://github.com/raybellwaves/xskillscore/issues/33 as the package is growing I wonder how we can streamline the package. Right now the code base is getting larger and larger. To navigate around it might help to consolidate.

Some ideas:

• Testing vectorized metrics: Have a general function that checks whether _metric from np_deterministic gives the same results as metric from deterministic
• Testing dask: Have a general function that checks whether functions in deterministic give the same results lazy as for .load()

Challenges:

• Can have think of a template for new metrics? Many metrics are very similar (the math is one line of code) but the code to implement a metric is more than 100 lines of code

There is a ton of warnings. They appear in pytest and the demo. So far I mostly suppressed all warnings. But there are catchers for warnings, which should be implemented.

seeing warnings in pytest:

pytest -r w

I had to uninstall pytest-tldr and pytest-sugar.

how to reduce: https://docs.pytest.org/en/stable/warnings.html

...and I think should be range(3) rather than len(3)?

Library looks exciting though!

I wonder if it's good to have a util that allows the user to apply multiple metrics at once in parallel?
Like:
xs.score(a, b, 'time', ['rmse', 'mae', 'pearson_r', 'pearson_r_p_val'])

https://black.readthedocs.io/en/stable/

Docs state scalars are ok

    mu : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
     scalars, Mix of labeled and/or unlabeled forecasts mean arrays.
    sig : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or
     scalars, Mix of labeled and/or unlabeled forecasts mean arrays.

However,

/mnt/c/Users/Solactus/GOOGLE~1/Bash/xskillscore/xskillscore/core/probabilistic.py in xr_crps_gaussian(observations, mu, sig)
     27     """
     28     # check if same dimensions
---> 29     if mu.dims != observations.dims:
     30         observations, mu = xr.broadcast(observations, mu)
     31     if sig.dims != observations.dims:

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

Probably need to do if not isinstance(mu, xr.DataArray): mu = xr.DataArray(mu)

Check whether deterministic metrics would speed up when omitting np.moveaxis or even writing in xarray functions.

see also http://xarray.pydata.org/en/stable/dask.html

import numpy as np
import xarray as xr
import bottleneck

def covariance_gufunc(x, y):
    return ((x - x.mean(axis=-1, keepdims=True))
            * (y - y.mean(axis=-1, keepdims=True))).mean(axis=-1)

def pearson_correlation_gufunc(x, y):
    return covariance_gufunc(x, y) / (x.std(axis=-1) * y.std(axis=-1))

def spearman_correlation_gufunc(x, y):
    x_ranks = bottleneck.rankdata(x, axis=-1)
    y_ranks = bottleneck.rankdata(y, axis=-1)
    return pearson_correlation_gufunc(x_ranks, y_ranks)

def spearman_correlation(x, y, dim):
    return xr.apply_ufunc(
        spearman_correlation_gufunc, x, y,
        input_core_dims=[[dim], [dim]],
        dask='parallelized',
        output_dtypes=[float])

https://github.com/raybellwaves/xskillscore/blob/c611d7bd722d4f06a8426de3a533cdbc63e234c4/xskillscore/core/deterministic.py#L41L44

via https://github.com/raybellwaves/xskillscore/pull/28#discussion_r338361930

Text: Should this check happen before weights = weights / weights.max()?

The docstrings for all deterministic functions claim that ndarrays (i.e. numpy arrays) can be passed through, but that isn't the case since a lot of the preprocessing and wrapper functions leverage xarray methods to set up the package. Given the prefix x for xskillscore, I think it would be reasonable to enforce that a and b have to be xarray objects.

We do this in climpred with a decorator:

https://github.com/bradyrx/climpred/blob/61e397bc07004e292af5b0798d4eca17be5dd263/climpred/checks.py#L19-L57

Which is then applied as a header to relevant functions:

https://github.com/bradyrx/climpred/blob/61e397bc07004e292af5b0798d4eca17be5dd263/climpred/prediction.py#L32-L41

This would be easy to port over to xskillscore and decorate all of the high-level functions with.

if forecast and obs have the same shape, then results just MAE.
forecast needs to have ensemble dimension extra for crps to be computed probabilistic.

https://github.com/raybellwaves/xskillscore/blob/c611d7bd722d4f06a8426de3a533cdbc63e234c4/xskillscore/core/deterministic.py#L51L52

via https://github.com/raybellwaves/xskillscore/pull/28#discussion_r338362278

Text: Could you combine this with the previous line to prevent looping over twice (although it's probably very minor)?

drop_dims = {k: 0 for k in a.dims if k not in new_dim}

as given here: https://peltarion.com/knowledge-center/documentation/modeling-view/build-an-ai-model/loss-functions/mean-squared-logarithmic-error

I'd like to calculate the spearman correlation.

Why? To calc correlation for non-gaussian quantities like mixed layer depth

Ref: Servonnat, Jérôme, Juliette Mignot, Eric Guilyardi, Didier Swingedouw, Roland Séférian, and Sonia Labetoulle. “Reconstructing the Subsurface Ocean Decadal Variability Using Surface Nudging in a Perfect Model Framework.” Climate Dynamics 44, no. 1–2 (January 1, 2015): 315–38. https://doi.org/10/f6v7kq.

Example code from xarray docu: http://xarray.pydata.org/en/stable/dask.html#automatic-parallelization

Wiki: https://en.wikipedia.org/wiki/Spearman%27s_rank_correlation_coefficient

Ways to implement:

• hack pearson_r with x_ranks = bottleneck.rankdata(x, axis=-1) and then rename to corr(type='pearson')
• write a new function the same way as pearson_r called spearman_r

I think something is misaligned

Check mark means it matches my expectation.

The highlighted part / X mark is fishy.

import xarray as xr
import xskillscore as xs

da = xr.DataArray([[[0, 1, 2], [np.nan, 4, 5]],
                   [[0, 1, 2], [3, np.nan, 5]],
                   [[0, 1, 2], [2, 4, 5]]
                  ], dims=['time', 'lat', 'lon'])
da.isel(time=0)
da.isel(time=1)
da.isel(time=2)

xs.pearson_r(da, da, dim='time', skipna=False)
xs.pearson_r(da, da, dim='time', skipna=True)

xs.pearson_r(da, da, dim=['lat', 'lon'], skipna=False)
xs.pearson_r(da, da, dim=['lat', 'lon'], skipna=True)

Create a requirements-py37.yml and push through to CI

http://xarray.pydata.org/en/stable/examples/apply_ufunc_vectorize_1d.html

I think it would be a nice feature to allow optional keyword arguments (kwarg) to

def mae(a, b, dim):

e.g. something like:

def mae(a, b, dim, kwargs=kwargs):

Which should have to get passed through to

    return xr.apply_ufunc(_mae, a, b,
                          input_core_dims=[dim, dim],
                          kwargs={'axis': axis},
                          dask='parallelized',
output_dtypes=[float])

e.g. something like:

    return xr.apply_ufunc(_mae, a, b,
                          input_core_dims=[dim, dim],
                          kwargs={'axis': axis}.update(kwargs),
                          dask='parallelized',
output_dtypes=[float])

A good start could be brier_score which can be adapted from https://github.com/TheClimateCorporation/properscoring

as was stated in the original conversation of starting up this repo: https://groups.google.com/forum/#!topic/xarray/z8ue0G-BLc8

when playing about around in https://github.com/raybellwaves/xskillscore/pull/66 I tried to get rid of the lines with np.rollaxis in correlations. I think it might not be needed and might then speedup even more, but I cannot make it work to pass the tests.

https://github.com/raybellwaves/xskillscore/blob/dd1400fd96afb82895cab881978d293942230bf4/xskillscore/core/np_deterministic.py#L96

Below some observations:

• somehow it doesnt change testing whether kwargs={"axis": -1, "skipna": skipna} or kwargs={"axis": 0, "skipna": skipna} in deterministic.py for correlations.
• apply_ufunc: Core dimensions are automatically moved to the last axes of input variables before applying func, which facilitates using NumPy style generalized ufuncs. http://xarray.pydata.org/en/stable/generated/xarray.apply_ufunc.html
• I would have guessed that we can spare command rollaxis if axis is axis=0

a = np.ones((3,4,5,6))
np.rollaxis(a,-1).shape
# equivalent
np.rollaxis(a,-1,0).shape

# what we now do
%timeit np.sum(np.rollaxis(a,-1),axis=0).shape 8.59 µs ± 251 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

# what might be faster
%timeit np.sum(a,axis=-1).shape 5.9 µs ± 362 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

do you see a solution to this potential speed up?

currently pandas is used to set the time dimension in testing. xr.cftime_range (http://xarray.pydata.org/en/stable/generated/xarray.cftime_range.html) can also do this.

Not sure if it is worth pursing but I recently came across https://github.com/pandas-profiling/pandas-profiling

It provides a few correlation metrics (e.g. https://pandas-profiling.github.io/pandas-profiling/examples/titanic/titanic_report.html#correlations) and has some categorical correlation metrics. I believe the source is: https://github.com/pandas-profiling/pandas-profiling/blob/master/src/pandas_profiling/model/correlations.py

I ran asv on my machine and published the results to https://raybellwaves.github.io/xskillscore/

@aaronspring sent a link to code that looks though it may be able to help automatically keep track of benchmarks with new commits: https://github.com/rmcgibbo/openmm-bench/blob/77b5ca394a64fcd0bbfe78ce73fe93045d93b482/run.sh

Not sure how to set this up.

It would be helpful to have a formal CHANGELOG, which I imagine would go in with https://github.com/raybellwaves/xskillscore/issues/36. Then it becomes clear to the user what each version adds.

I calc skill on chunked and non-chunked inputs and get tiny O(-16) differences. Is this worrying?

The differences are super tiny, but I would not have expected this.
This happens for all functions with xr.apply_ufunc(). It therefore only doesnt happen for xr_mse?

See demo:
https://gist.github.com/aaronspring/d14df048b89c441a331b1b451678efce

Do you have ideas about this?
Is it due to parallelized?

Looking through https://github.com/Azure/MachineLearningNotebooks/blob/master/tutorials/regression-automated-ml.ipynb

and I noticed they have accuracy = 1 - MAPE

As given here: https://www.kaggle.com/c/rossmann-store-sales/overview/evaluation

xskillscore includes median absolute deviation (MAD). It seems to be returning different results than other packages and is not written as Wikipedia (https://en.wikipedia.org/wiki/Median_absolute_deviation) or scipy (https://scipy.github.io/devdocs/generated/scipy.stats.median_absolute_deviation.html) would suggest it to be. I also can't find any documentation on it in forecast literature/pages.

import numpy as np
import pandas as pd
import xarray as xr
import xskillscore as xs
import astropy.stats
import scipy.stats

np.random.seed(12345)
def generate_data():
    time = pd.date_range("1/1/2000", "1/5/2000", freq="D")
    da = xr.DataArray(np.random.rand(len(time)), dims=["time"], coords=[time])
    return da

a = generate_data()
b = generate_data()

xs.mad(a, b, "time")
<xarray.DataArray ()>
array(0.46925829)

astropy.stats.median_absolute_deviation(a, b)
0.13518689535658013

scipy.stats.median_absolute_deviation(a, b)
0.2004280910556657

The scipy description suggests that MAD is the median of the median error. I.e. np.median(x - np.median(x)). It seems like a statistic for a single time series. Does it make sense to apply this to the absolute error of A and B?

Just checking that there's some evidence of this being used in the field.

CC @raybellwaves, @aaronspring

As discussed in hdrake/cmip6hack-multigen#1
we'd like to be able to pass weights (e.g. grid cell areas) to the various skill score metrics.

I'll hopefully submit a PR for this soon.

Adding flake8 will ensure the code adheres to PEP8 standards.

Some info here
https://simpleisbetterthancomplex.com/packages/2016/08/05/flake8.html

For example, I want to take the RMSE across the lat and lon dimension for one time step.

Maybe it’s easier to understand what skipna does on time, lon, lat data when we have a notebook based on xr.tutorial.rasm dataset with some xskillscore metrics. then we might also see directly how land is treated in correlations which is not too obvious in https://github.com/raybellwaves/xskillscore/pull/38

It would be nice to add proper sphinx documentation for xskillscore hosted on readthedocs. A really barebones version is hosted at https://github.com/bradyrx/esmtools/, which gives you the essentials. We have more expansive documentation at https://github.com/bradyrx/climpred, which shows you how to use notebooks in the docs.

I think just having the API hosted more formally + some example notebooks would be great. @hdrake expressed interest in leading this. However, I think readthedocs needs to be run by the owner of the repo (@raybellwaves). Perhaps the sphinx can be assembled by someone and then the readthedocs settings set up by @raybellwaves?

d1 = xr.tutorial.open_dataset('air_temperature').isel(time=slice(0, 5))
d2 = xr.tutorial.open_dataset('air_temperature').isel(time=slice(5, 10))
d2 = d2.where(d2['air'] > 279)
d2['time'] = d1['time']
xs.pearson_r_p_value(d1, d2, ['lat', 'lon'], skipna=True)

air      (time) float64 nan nan nan nan nan

This works though

xs.spearman_r_p_value(d1, d2, ['lat', 'lon'], skipna=True)

Signature: xs.rmse(a, b, dim)
Docstring:
Root Mean Squared Error.

Parameters

a : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or scalars
Mix of labeled and/or unlabeled arrays to which to apply the function.
b : Dataset, DataArray, GroupBy, Variable, numpy/dask arrays or scalars
Mix of labeled and/or unlabeled arrays to which to apply the function.
dim : str
The dimension to apply the correlation along.

Returns

Single value or tuple of Dataset, DataArray, Variable, dask.array.Array or
numpy.ndarray, the first type on that list to appear on an input.
Root Mean Squared Error.

Currently, the _match_nans function in https://github.com/raybellwaves/xskillscore/blob/master/xskillscore/core/np_deterministic.py requires a and b have identical dimensions, while most operations on xarray.Dataset objects broadcast across non-common dimensions. This means that all dependent skill scores (e.g. _mae) are inconsistent between skipna=True and skipna=False options.

def _match_nans(a, b, weights):
    """
    Considers missing values pairwise. If a value is missing
    in a, the corresponding value in b is turned to nan, and
    vice versa.
    Returns
    -------
    a, b, weights : ndarray
        a, b, and weights (if not None) with nans placed at
        pairwise locations.
    """
    if np.isnan(a).any() or np.isnan(b).any():
        # Find pairwise indices in a and b that have nans.
        idx = np.logical_or(np.isnan(a), np.isnan(b))
        a[idx], b[idx] = np.nan, np.nan
        if weights is not None:
            weights[idx] = np.nan
    return a, b, weights

I think the following modifications could work to address this:

def _match_nans(a, b, weights):
    """
    Considers missing values pairwise. If a value is missing
    in a, the corresponding value in b is turned to nan, and
    vice versa.
    Returns
    -------
    a, b, weights : ndarray
        a, b, and weights (if not None) with nans placed at
        pairwise locations.
    """
    if np.isnan(a).any() or np.isnan(b).any():
        a, b = xr.broadcast(a, b)
        # Find pairwise indices in a and b that have nans.
        idx = np.logical_or(np.isnan(a), np.isnan(b))
        a[idx], b[idx] = np.nan, np.nan
        if weights is not None:
            weights[idx] = np.nan
    return a, b, weights

I'll try to test this myself and put in a Pull Request sometime soon but just wanted to open this issue first.

I removed '('time', 'lat', 'lon')' from the AXES variable in xskillscore/tests/test_deterministic.py

As I got the issue:
test_deterministic.py::test_pearson_r_xr_dask[dim4] Aborted (core dumped)
isn't running.

I'm guessing it has something to do with
https://docs.dask.org/en/latest/changelog.html#id1

I would argue that this probably isn't worth looking into as I there are other packages to do stats on vectors. The package focuses on ndarrays.

Probably need to provide the user with a warning.

>>> import xarray as xr
>>> import pandas as pd
>>> import numpy as np
>>> from scipy.stats import norm
>>> import xskillscore as xs
>>> obs = xr.DataArray(
...     np.random.rand(3, 4, 5),
...     coords=[
...         pd.date_range("1/1/2000", "1/3/2000", freq="D"),
...         np.arange(4),
...         np.arange(5),
...     ],
...     dims=["time", "lat", "lon"],
... )
>>> fct = obs.copy()
>>> fct.values = np.random.rand(3, 4, 5)

>>> weights = np.cos(np.deg2rad(obs.lat))
>>> _, weights = xr.broadcast(obs, weights)
>>> weights = weights.isel(time=0)
>>> r = xs.pearson_r(obs, fct, ["time", "lat", "lon"], weights=weights)
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/home/ray/Documents/PYTHON_dev/xskillscore/xskillscore/core/deterministic.py", line 120, in pearson_r
    weights = weights.stack(**{new_dim: dim})
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/dataarray.py", line 1725, in stack
    ds = self._to_temp_dataset().stack(dimensions, **dimensions_kwargs)
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/dataset.py", line 3234, in stack
    result = result._stack_once(dims, new_dim)
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/dataset.py", line 3181, in _stack_once
    shape = [self.dims[d] for d in vdims]
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/dataset.py", line 3181, in <listcomp>
    shape = [self.dims[d] for d in vdims]
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/utils.py", line 386, in __getitem__
    return self.mapping[key]
  File "/home/ray/local/bin/anaconda3/envs/xss/lib/python3.6/site-packages/xarray/core/utils.py", line 417, in __getitem__
    return self.mapping[key]
KeyError: 'time'

import numpy as np
import xarray as xr
import xskillscore as xs

da = xr.DataArray([0, 1, 2], dims=['time'])

print('-----')

xs.pearson_r(da, da, dim='time', skipna=True)

yields TypeError: int() argument must be a string, a bytes-like object or a number, not 'NoneType' because integers cannot be NaN

Since xskillscore requires properscoring which requires numba, maybe can implement numba for native xskillscore numpy functions, similar to: https://github.com/TheClimateCorporation/properscoring/blob/master/properscoring/_gufuncs.py#L5

https://en.wikipedia.org/wiki/Cram%C3%A9r's_V?wprov=sfti1

Believe it’s used in ML to show correlation between variables. Has range 0 - 1

https://scikit-learn.org/stable/modules/generated/sklearn.feature_selection.chi2.html

Instead of digging through the growing README, I think the following should be displayed up front.

all = [
'_pearson_r',
'_pearson_r_p_value',
'_rmse',
'_mse',
'_mae',
'_mad',
'_smape',
'_mape',
'_spearman_r',
'_spearman_r_p_value',
]
all = [
'brier_score',
'crps_ensemble',
'crps_gaussian',
'crps_quadrature',
'threshold_brier_score',
]

See how it is done in xarray

https://github.com/pydata/xarray/blob/d089df385e737f71067309ff7abae15994d581ec/xarray/_version.py

https://github.com/pydata/xarray/blob/d089df385e737f71067309ff7abae15994d581ec/xarray/__init__.py#L4

@aaronspring and @ahuang11 anything else you want to add before I bump up the version to take into account the accessors addition and bug fix to xr_crps_gaussian?

I’m on the organizing committee for PyData Miami. https://pydata.org/miami2020/

Interested in attending @bradyrx @aaronspring @ahuang11?

The conference is on Wednesday-Thursday. If interested we could look at doing a sprint on the Friday. We could do it at RSMAS so you could also tie in a research visit

I dont want to be the party pooper, do you get an increase in performance with xr.metric over a easily self-written xr function? I don't on my laptop.
For comparison, I also wrapped MSE from sklearn with apply_ufunc.

# test bed for xs performance

import numpy as np
import pandas as pd
import sklearn
import xarray as xr
import xskillscore as xs
from sklearn.metrics import mean_squared_error

np.random.seed(12345)


def generate_small_data():
    time = pd.date_range("1/1/2000", "1/5/2000", freq="D")
    da = xr.DataArray(np.random.rand(len(time)), dims=["time"], coords=[time])
    return da


a = generate_small_data()
b = generate_small_data()


def generate_large_data():
    time = pd.date_range("1/1/2000", "1/5/2000", freq="D")
    nlon, nlat = 200, 200
    da = xr.DataArray(np.random.rand(len(time), nlon, nlat), dims=[
                      "time", "lon", "lat"], coords=[time, np.arange(nlon), np.arange(nlat)])
    return da


#a = generate_large_data()
#b = generate_large_data()


def xr_mse(a, b, dim):
    return ((a-b)**2).mean(dim)

def sk_mse(a, b, dim):
    # doesnt work yet for ndims > 2
    # tried here
    if len(a.dims) > 2:
        a = a.stack(s=a.drop(dim).dims)
        b = b.stack(s=b.drop(dim).dims)
        stacked = True
    else:
        stacked = False
    res = xr.apply_ufunc(
        # sklearn.metrics.mean_squared_error
        mean_squared_error,
        a,
        b,
        # weights,
        input_core_dims=[[dim], [dim]],
        # kwargs={"axis": axis, "skipna": skipna},
        dask="parallelized",
        output_dtypes=[float],
    )
    if stacked:
        res = res.unstack()
    return res

%timeit xs.mse(a, b, "time")
723 µs ± 7.79 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit xr_mse(a, b, "time")
383 µs ± 13.6 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

%timeit sk_mse(a, b, 'time')
619 µs ± 59.4 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

I re-ran timeit several times to avoid cache effects or compilation effects. A better way to test this would be directly with asv. I get similar results for larger datasets.

Why is xr performance higher? (xr).mean, (xr).sum, ... is automatically vectorized and basically from numpy.

Called mape.

See https://github.com/facebook/prophet/blob/c179a28be2a1dc4e698a82a5ca5f2aa7724818bf/python/fbprophet/diagnostics.py#L441L453

and https://stackoverflow.com/a/47648179/6046019

np.mean(np.abs((y - yhat) / y)) * 100

https://azure.microsoft.com/en-us/services/devops/pipelines/

See https://github.com/pydata/xarray/blob/master/azure-pipelines.yml

Add some tick boxes for ensuring tests are written and flake8 has been run on the repo. https://github.com/raybellwaves/xskillscore/pull/28#issuecomment-544214094

There may be a better way to setup flake8 in CI...

When implementing properscoring.brier_score, I saw that we dont test lazy computation, i.e. whether we get chunked results when entering chunked inputs. Therefore I added a few assert statements. But this error also occurs now for test_xr_threshold_brier_score_dask which had the line assert actual.chunks is not None already before I did my changes.

    @pytest.mark.parametrize('dim', AXES)
    def test_pearson_r_p_value_xr_dask(a_dask, b_dask, dim):
        actual = pearson_r_p_value(a_dask, b_dask, dim)
>       assert actual.chunks is not None
E       assert None is not None

Maybe it has something to do with the stacking in _preprocess? @ahuang11

Furthermore I get this warning, maybe there is a problem with np.rollaxis:

======================================= warnings summary ========================================
/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.6/site-packages/_pytest/mark/structures.py:332
  /Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.6/site-packages/_pytest/mark/structures.py:332: PytestUnknownMarkWarning: Unknown pytest.mark.flaky - is this a typo?  You can register custom marks to avoid this warning - for details, see https://docs.pytest.org/en/latest/mark.html
    PytestUnknownMarkWarning,

/Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.6/site-packages/_pytest/mark/structures.py:332
  /Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.6/site-packages/_pytest/mark/structures.py:332: PytestUnknownMarkWarning: Unknown pytest.mark.network - is this a typo?  You can register custom marks to avoid this warning - for details, see https://docs.pytest.org/en/latest/mark.html
    PytestUnknownMarkWarning,

xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[time]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[time]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[lat]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[lat]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[lon]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[lon]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[dim3]
xskillscore/tests/test_deterministic.py::test_pearson_r_xr_dask[dim3]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[time]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[time]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[time]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lat]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lat]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lat]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lon]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lon]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[lon]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[dim3]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[dim3]
xskillscore/tests/test_deterministic.py::test_pearson_r_p_value_xr_dask[dim3]
  /Users/aaron.spring/anaconda3/envs/climpred-dev/lib/python3.6/site-packages/dask/array/core.py:1263: FutureWarning: The `numpy.rollaxis` function is not implemented by Dask array. You may want to use the da.map_blocks function or something similar to silence this warning. Your code may stop working in a future release.
    FutureWarning,

-- Docs: https://docs.pytest.org/en/latest/warnings.html

See https://github.com/raybellwaves/xskillscore/pull/18 for the new tests.

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)