Source code for climpred.metrics

import warnings

import numpy as np
import xarray as xr
from scipy.stats import norm
from xskillscore import (
    brier_score,
    crps_ensemble,
    crps_gaussian,
    crps_quadrature,
    mae,
    mse,
    pearson_r,
    pearson_r_p_value,
    rmse,
    threshold_brier_score,
)

from .constants import CLIMPRED_DIMS


def _get_norm_factor(comparison):
    """Get normalization factor for PPP, NMSE, NRMSE, MSSS.

    Used in compute_perfect_model. Comparison 'm2e' gets smaller rmse's than
    'm2m' by design, see Seferian et al. 2018. 'm2m', 'm2c' ensemble variance
    is divided by 2 to get control variance.

    Args:
        comparison (function): comparison function.

    Returns:
        fac (int): normalization factor.

    Raises:
        KeyError: if comparison is not matching.

    """
    comparison_name = comparison.__name__
    if comparison_name in ['_m2e', '_e2c', '_e2r']:
        fac = 1
    elif comparison_name in ['_m2c', '_m2m', '_m2r']:
        fac = 2
    else:
        raise KeyError('specify comparison to get normalization factor.')
    return fac


[docs]def _pearson_r(forecast, reference, dim='svd', **metric_kwargs): """ Calculate the Anomaly Correlation Coefficient (ACC). .. math:: ACC = \\frac{cov(f, o)}{\\sigma_{f}\\cdot\\sigma_{o}} .. note:: Use metric ``pearson_r_p_value`` to get the corresponding pvalue. Range: * perfect: 1 * min: -1 See also: * xskillscore.pearson_r * xskillscore.pearson_r_p_value """ return pearson_r(forecast, reference, dim=dim)
def _pearson_r_p_value(forecast, reference, dim='svd', **metric_kwargs): """ Calculate the probability associated with the ACC not being random. """ # p-value returns a runtime error when working with NaNs, such as on a climate # model grid. We can avoid this annoying output by specifically suppressing # warning here. with warnings.catch_warnings(): warnings.simplefilter('ignore') pval = pearson_r_p_value(forecast, reference, dim=dim) return pval
[docs]def _mse(forecast, reference, dim='svd', **metric_kwargs): """ Calculate the Mean Sqaure Error (MSE). .. math:: MSE = \\overline{(f - o)^{2}} Range: * perfect: 0 * min: 0 * max: ∞ See also: * xskillscore.mse """ return mse(forecast, reference, dim=dim)
[docs]def _rmse(forecast, reference, dim='svd', **metric_kwargs): """ Calculate the Root Mean Sqaure Error (RMSE). .. math:: RMSE = \\sqrt{\\overline{(f - o)^{2}}} Range: * perfect: 0 * min: 0 * max: ∞ See also: * xskillscore.rmse """ return rmse(forecast, reference, dim=dim)
[docs]def _mae(forecast, reference, dim='svd', **metric_kwargs): """ Calculate the Mean Absolute Error (MAE). .. math:: MSE = \\overline{(f - o)^{2}} Range: * perfect: 0 * min: 0 * max: ∞ See also: * xskillscore.mae """ return mae(forecast, reference, dim=dim)
[docs]def _brier_score(forecast, reference, **metric_kwargs): """Calculate Brier score for forecasts on binary reference. ..math: BS(f, o) = (f - o)^2 Args: * forecast (xr.object) * reference (xr.object) * func (function): function to be applied to reference and forecasts and then mean('member') to get forecasts and reference in interval [0,1]. (required to be added via **metric_kwargs) Reference: * Brier, Glenn W. “VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF PROBABILITY.” Monthly Weather Review 78, no. 1 (1950). https://doi.org/10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2. Example: >>> def pos(x): return x > 0 >>> compute_perfect_model(ds, control, metric='brier_score', func=pos) See also: * properscoring.brier_score * xskillscore.brier_score """ if 'func' in metric_kwargs: func = metric_kwargs['func'] else: raise ValueError( 'Please provide a function `func` to be applied to comparison and \ reference to get values in interval [0,1]; \ see properscoring.brier_score.' ) return brier_score(func(reference), func(forecast).mean('member'))
[docs]def _threshold_brier_score(forecast, reference, **metric_kwargs): """ Calculate the Brier scores of an ensemble for exceeding given thresholds. Provide threshold via metric_kwargs. .. math:: CRPS(F, x) = \int_z BS(F(z), H(z - x)) dz Range: * perfect: 0 * min: 0 * max: 1 Args: * forecast (xr.object) * reference (xr.object) * threshold (int, float, xr.object): Threshold to check exceedance, see properscoring.threshold_brier_score (required to be added via **metric_kwargs) References: * Brier, Glenn W. “VERIFICATION OF FORECASTS EXPRESSED IN TERMS OF PROBABILITY.” Monthly Weather Review 78, no. 1 (1950). https://doi.org/10.1175/1520-0493(1950)078<0001:VOFEIT>2.0.CO;2. Example: >>> compute_perfect_model(ds, control, metric='threshold_brier_score', threshold=.5) See also: * properscoring.threshold_brier_score * xskillscore.threshold_brier_score """ if 'threshold' not in metric_kwargs: raise ValueError('Please provide threshold.') else: threshold = metric_kwargs['threshold'] # switch args b/c xskillscore.threshold_brier_score(obs, forecasts) return threshold_brier_score(reference, forecast, threshold)
[docs]def _crps(forecast, reference, **metric_kwargs): """ Continuous Ranked Probability Score (CRPS) is the probabilistic MSE. Range: * perfect: 0 * min: 0 * max: ∞ References: * Matheson, James E., and Robert L. Winkler. “Scoring Rules for Continuous Probability Distributions.” Management Science 22, no. 10 (June 1, 1976): 1087–96. https://doi.org/10/cwwt4g. See also: * properscoring.crps_ensemble * xskillscore.crps_ensemble """ # switch positions because xskillscore.crps_ensemble(obs, forecasts) return crps_ensemble(reference, forecast)
def _crps_gaussian(forecast, mu, sig, **metric_kwargs): """CRPS assuming a gaussian distribution. Helper function for CRPSS. See also: * properscoring.crps_gaussian * xskillscore.crps_gaussian """ return crps_gaussian(forecast, mu, sig) def _crps_quadrature( forecast, cdf_or_dist, xmin=None, xmax=None, tol=1e-6, **metric_kwargs ): """CRPS assuming distribution cdf_or_dist. Helper function for CRPSS. See also: * properscoring.crps_quadrature * xskillscore.crps_quadrature """ return crps_quadrature(forecast, cdf_or_dist, xmin, xmax, tol)
[docs]def _crpss(forecast, reference, **metric_kwargs): """Continuous Ranked Probability Skill Score .. note:: When assuming a gaussian distribution of forecasts, use default gaussian=True. If not gaussian, you may specify the distribution type, xmin/xmax/tolerance for integration (see xskillscore.crps_quadrature). .. math:: CRPSS = 1 - \\frac{CRPS_{init}}{CRPS_{clim}} Args: * forecast (xr.object): * reference (xr.object): * gaussian (bool): Assuming gaussian distribution for baseline skill. Default: True (optional) * cdf_or_dist (scipy.stats): distribution to assume if not gaussian. default: scipy.stats.norm * xmin, xmax, tol: only relevant if not gaussian (see xskillscore.crps_quadrature) Range: * perfect: 1 * pos: better than climatology forecast * neg: worse than climatology forecast References: * Matheson, James E., and Robert L. Winkler. “Scoring Rules for Continuous Probability Distributions.” Management Science 22, no. 10 (June 1, 1976): 1087–96. https://doi.org/10/cwwt4g. * Gneiting, Tilmann, and Adrian E Raftery. “Strictly Proper Scoring Rules, Prediction, and Estimation.” Journal of the American Statistical Association 102, no. 477 (March 1, 2007): 359–78. https://doi.org/10/c6758w. Example: >>> compute_perfect_model(ds, control, metric='crpss') >>> compute_perfect_model(ds, control, metric='crpss', gaussian=False, cdf_or_dist=scipy.stats.norm, xmin=-10, xmax=10, tol=1e-6) See also: * properscoring.crps_ensemble * xskillscore.crps_ensemble """ # available climpred dimensions to take mean and std over rdim = [tdim for tdim in reference.dims if tdim in CLIMPRED_DIMS] mu = reference.mean(rdim) sig = reference.std(rdim) # checking metric_kwargs, if not found use defaults: gaussian, else crps_quadrature if 'gaussian' in metric_kwargs: gaussian = metric_kwargs['gaussian'] else: gaussian = True if gaussian: ref_skill = _crps_gaussian(forecast, mu, sig) else: if 'cdf_or_dist' in metric_kwargs: cdf_or_dist = metric_kwargs['cdf_or_dist'] else: cdf_or_dist = norm if 'xmin' in metric_kwargs: xmin = metric_kwargs['xmin'] else: xmin = None if 'xmax' in metric_kwargs: xmax = metric_kwargs['xmax'] else: xmax = None if 'tol' in metric_kwargs: tol = metric_kwargs['tol'] else: tol = 1e-6 ref_skill = _crps_quadrature(forecast, cdf_or_dist, xmin, xmax, tol) forecast_skill = _crps(forecast, reference) skill_score = 1 - forecast_skill / ref_skill.mean('member') return skill_score
[docs]def _crpss_es(forecast, reference, **metric_kwargs): """CRPSS Ensemble Spread. .. math:: CRPSS = 1 - \\frac{CRPS(\\sigma^2_f)}{CRPS(\\sigma^2_o})) References: * Kadow, Christopher, Sebastian Illing, Oliver Kunst, Henning W. Rust, Holger Pohlmann, Wolfgang A. Müller, and Ulrich Cubasch. “Evaluation of Forecasts by Accuracy and Spread in the MiKlip Decadal Climate Prediction System.” Meteorologische Zeitschrift, December 21, 2016, 631–43. https://doi.org/10/f9jrhw. Range: * perfect: 0 * else: negative """ # helper dim to calc mu rdim = [tdim for tdim in reference.dims if tdim in CLIMPRED_DIMS + ['time']] # inside compute_perfect_model if 'init' in forecast.dims: dim2 = 'init' # inside compute_hindcast elif 'time' in forecast.dims: dim2 = 'time' else: raise ValueError('dim2 not found automatically in ', forecast.dims) mu = reference.mean(rdim) forecast, ref2 = xr.broadcast(forecast, reference) sig_r = _mse(forecast, ref2, dim='member').mean(dim2) sig_h = _mse(forecast.mean(dim2), ref2.mean(dim2), 'member') crps_h = _crps_gaussian(forecast, mu, sig_h) if 'member' in crps_h.dims: crps_h = crps_h.mean('member') crps_r = _crps_gaussian(forecast, mu, sig_r) if 'member' in crps_r.dims: crps_r = crps_r.mean('member') return 1 - crps_h / crps_r
[docs]def _bias(forecast, reference, dim='svd', **metric_kwargs): """Calculate unconditional bias. .. math:: bias = f - o Range: * pos: positive bias * neg: negative bias * perfect: 0 References: * https://www.cawcr.gov.au/projects/verification/ * https://www-miklip.dkrz.de/about/murcss/ """ bias = (forecast - reference).mean(dim) return bias
[docs]def _msss_murphy(forecast, reference, dim='svd', **metric_kwargs): """Calculate Murphy's Mean Square Skill Score (MSSS). .. math:: MSSS_{Murphy} = r_{fo}^2 - [\\text{conditional bias}]^2 -\ [\\frac{\\text{(unconditional) bias}}{\\sigma_o}]^2 References: * https://www-miklip.dkrz.de/about/murcss/ * Murphy, Allan H. “Skill Scores Based on the Mean Square Error and Their Relationships to the Correlation Coefficient.” Monthly Weather Review 116, no. 12 (December 1, 1988): 2417–24. https://doi.org/10/fc7mxd. """ acc = _pearson_r(forecast, reference, dim=dim) conditional_bias = _conditional_bias(forecast, reference, dim=dim) uncond_bias = _bias(forecast, reference, dim=dim) / reference.std(dim) skill = acc ** 2 - conditional_bias ** 2 - uncond_bias ** 2 return skill
[docs]def _conditional_bias(forecast, reference, dim='svd', **metric_kwargs): """Calculate the conditional bias between forecast and reference. .. math:: \\text{conditional bias} = r_{fo} - \\frac{\\sigma_f}{\\sigma_o} References: * https://www-miklip.dkrz.de/about/murcss/ """ acc = _pearson_r(forecast, reference, dim=dim) conditional_bias = acc - _std_ratio(forecast, reference, dim=dim) ** -1 return conditional_bias
[docs]def _std_ratio(forecast, reference, dim='svd', **metric_kwargs): """Calculate the ratio of standard deviations of reference over forecast. .. math:: \\text{std ratio} = \\frac{\\sigma_o}{\\sigma_f} References: * https://www-miklip.dkrz.de/about/murcss/ """ ratio = reference.std(dim) / forecast.std(dim) return ratio
[docs]def _bias_slope(forecast, reference, dim='svd', **metric_kwargs): """Calculate bias slope between reference and forecast standard deviations. .. math:: \\text{bias slope}= r_{fo} \\cdot \\text{std ratio} References: * https://www-miklip.dkrz.de/about/murcss/ """ std_ratio = _std_ratio(forecast, reference, dim=dim) acc = _pearson_r(forecast, reference, dim=dim) b_s = std_ratio * acc return b_s
[docs]def _ppp(forecast, reference, dim='svd', **metric_kwargs): """Prognostic Potential Predictability (PPP) metric. .. math:: PPP = 1 - \\frac{MSE}{ \\sigma^2_{ref} \\cdot fac} Args: * forecast (xr.object) * reference (xr.object) * dim (str): dimension to apply metric to * comparison (str): name comparison needed for normalization factor (required to be added via **metric_kwargs) Range: * 1: perfect forecast * positive: better than climatology forecast * negative: worse than climatology forecast References: * Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4. * Pohlmann, Holger, Michael Botzet, Mojib Latif, Andreas Roesch, Martin Wild, and Peter Tschuck. “Estimating the Decadal Predictability of a Coupled AOGCM.” Journal of Climate 17, no. 22 (November 1, 2004): 4463–72. https://doi.org/10/d2qf62. * Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill. Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq. """ mse_skill = _mse(forecast, reference, dim=dim) var = reference.var(dim) if 'comparison' in metric_kwargs: comparison = metric_kwargs['comparison'] else: raise ValueError( 'Comparison needed to normalize PPP. Not found in', metric_kwargs ) fac = _get_norm_factor(comparison) ppp_skill = 1 - mse_skill / var / fac return ppp_skill
[docs]def _nrmse(forecast, reference, dim='svd', **metric_kwargs): """Normalized Root Mean Square Error (NRMSE) metric. .. math:: NRMSE = \\frac{RMSE}{\\sigma_{o} \\cdot \\sqrt{fac} } = \\sqrt{ \\frac{MSE}{ \\sigma^2_{o} \\cdot fac} } Args: * forecast (xr.object) * reference (xr.object) * dim (str): dimension to apply metric to * comparison (str): name comparison needed for normalization factor (required to be added via **metric_kwargs) Range: * 0: perfect forecast * 0 - 1: better than climatology forecast * > 1: worse than climatology forecast References: * Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill.” Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq. * Hawkins, Ed, Steffen Tietsche, Jonathan J. Day, Nathanael Melia, Keith Haines, and Sarah Keeley. “Aspects of Designing and Evaluating Seasonal-to-Interannual Arctic Sea-Ice Prediction Systems.” Quarterly Journal of the Royal Meteorological Society 142, no. 695 (January 1, 2016): 672–83. https://doi.org/10/gfb3pn. """ rmse_skill = _rmse(forecast, reference, dim=dim) std = reference.std(dim) if 'comparison' in metric_kwargs: comparison = metric_kwargs['comparison'] else: raise ValueError( 'Comparison needed to normalize NRMSE. Not found in', metric_kwargs ) fac = _get_norm_factor(comparison) nrmse_skill = rmse_skill / std / np.sqrt(fac) return nrmse_skill
[docs]def _nmse(forecast, reference, dim='svd', **metric_kwargs): """ Calculate Normalized MSE (NMSE) = Normalized Ensemble Variance (NEV). .. math:: NMSE = NEV = \\frac{MSE}{\\sigma^2_{o} \\cdot fac} Args: * forecast (xr.object) * reference (xr.object) * dim (str): dimension to apply metric to * comparison (str): name comparison needed for normalization factor (required to be added via **metric_kwargs) Range: * 0: perfect forecast: 0 * 0 - 1: better than climatology forecast * > 1: worse than climatology forecast References: * Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4. """ mse_skill = _mse(forecast, reference, dim=dim) var = reference.var(dim) if 'comparison' in metric_kwargs: comparison = metric_kwargs['comparison'] else: raise ValueError( 'Comparison needed to normalize NMSE. Not found in', metric_kwargs ) fac = _get_norm_factor(comparison) nmse_skill = mse_skill / var / fac return nmse_skill
[docs]def _nmae(forecast, reference, dim='svd', **metric_kwargs): """ Normalized Ensemble Mean Absolute Error metric. .. math:: NMAE = \\frac{MAE}{\\sigma_{o} \\cdot fac} Args: * forecast (xr.object) * reference (xr.object) * dim (str): dimension to apply metric to * comparison (str): name comparison needed for normalization factor (required to be added via **metric_kwargs) Range: * 0: perfect forecast: 0 * 0 - 1: better than climatology forecast * > 1: worse than climatology forecast References: * Griffies, S. M., and K. Bryan. “A Predictability Study of Simulated North Atlantic Multidecadal Variability.” Climate Dynamics 13, no. 7–8 (August 1, 1997): 459–87. https://doi.org/10/ch4kc4. """ mae_skill = _mae(forecast, reference, dim=dim) std = reference.std(dim).mean() if 'comparison' in metric_kwargs: comparison = metric_kwargs['comparison'] else: raise ValueError( 'Comparison needed to normalize NMSE. Not found in', metric_kwargs ) fac = _get_norm_factor(comparison) nmae_skill = mae_skill / std / fac return nmae_skill
[docs]def _uacc(forecast, reference, dim='svd', **metric_kwargs): """ Calculate Bushuk's unbiased ACC (uACC). .. math:: uACC = \\sqrt{PPP} = \\sqrt{MSSS} Range: * 1: perfect * 0 - 1: better than climatology References: * Bushuk, Mitchell, Rym Msadek, Michael Winton, Gabriel Vecchi, Xiaosong Yang, Anthony Rosati, and Rich Gudgel. “Regional Arctic Sea–Ice Prediction: Potential versus Operational Seasonal Forecast Skill. Climate Dynamics, June 9, 2018. https://doi.org/10/gd7hfq. """ return _ppp(forecast, reference, dim=dim, **metric_kwargs) ** 0.5