Metrics

modelskill.metrics

The metrics module contains different skill metrics for evaluating the difference between a model and an observation.

• bias
• max_error
• root_mean_squared_error (rmse)
• urmse
• mean_absolute_error (mae)
• mean_absolute_percentage_error (mape)
• kling_gupta_efficiency (kge)
• nash_sutcliffe_efficiency (nse)
• r2 (r2=nse)
• model_efficiency_factor (mef)
• wilmott
• scatter_index (si)
• scatter_index2
• corrcoef (cc)
• spearmanr (rho)
• lin_slope
• hit_ratio
• explained_variance (ev)
• peak_ratio (pr)

Circular metrics (for directional data with units in degrees):

• c_bias
• c_max_error
• c_mean_absolute_error (c_mae)
• c_root_mean_squared_error (c_rmse)
• c_unbiased_root_mean_squared_error (c_urmse)

The names in parentheses are shorthand aliases for the different metrics.

Examples:

>>> obs = np.array([0.3, 2.1, -1.0])
>>> mod = np.array([0.0, 2.3, 1.0])
>>> bias(obs, mod)
np.float64(0.6333333333333332)
>>> max_error(obs, mod)
np.float64(2.0)
>>> rmse(obs, mod)
np.float64(1.173314393786536)
>>> urmse(obs, mod)
np.float64(0.9877021593352702)
>>> mae(obs, mod)
np.float64(0.8333333333333331)
>>> mape(obs, mod)
np.float64(103.17460317460316)
>>> nse(obs, mod)
np.float64(0.14786795048143053)
>>> r2(obs, mod)
np.float64(0.14786795048143053)
>>> mef(obs, mod)
np.float64(0.9231099877688299)
>>> si(obs, mod)
np.float64(0.8715019052958266)
>>> spearmanr(obs, mod)
np.float64(0.5)
>>> willmott(obs, mod)
np.float64(0.7484604452865941)
>>> hit_ratio(obs, mod, a=0.5)
np.float64(0.6666666666666666)
>>> ev(obs, mod)
np.float64(0.39614855570839064)

add_metric

add_metric(metric, has_units=False)

Adds a metric to the metric list. Useful for custom metrics.

Some metrics are dimensionless, others have the same dimension as the observations.

Parameters:

Name	Type	Description	Default
metric	str or callable	Metric name or function	required
has_units	bool	True if metric has a dimension, False otherwise. Default:False	False

Returns:

Type	Description
None

Examples:

>>> add_metric(hit_ratio)
>>> add_metric(rmse,True)

Source code in modelskill/metrics.py

def add_metric(metric: Callable, has_units: bool = False) -> None:
    """Adds a metric to the metric list. Useful for custom metrics.

    Some metrics are dimensionless, others have the same dimension as the observations.

    Parameters
    ----------
    metric : str or callable
        Metric name or function
    has_units : bool
        True if metric has a dimension, False otherwise. Default:False

    Returns
    -------
    None

    Examples
    --------
    >>> add_metric(hit_ratio)
    >>> add_metric(rmse,True)
    """
    defined_metrics.add(metric.__name__)
    if has_units:
        METRICS_WITH_DIMENSION.add(metric.__name__)

    # add the function to the module
    setattr(sys.modules[__name__], metric.__name__, metric)

bias

bias(obs, model)

Bias (mean error)

\[ bias=\frac{1}{n}\sum_{i=1}^n (model_i - obs_i) \]

Range: \((-\infty, \infty)\); Best: 0

Source code in modelskill/metrics.py

def bias(obs, model) -> Any:
    r"""Bias (mean error)

    $$
    bias=\frac{1}{n}\sum_{i=1}^n (model_i - obs_i)
    $$

    Range: $(-\infty, \infty)$; Best: 0
    """

    assert obs.size == model.size
    return np.mean(model - obs)

c_bias

c_bias(obs, model)

Circular bias (mean error)

Parameters:

Name	Type	Description	Default
obs	ndarray	Observation in degrees (0, 360)	required
model	ndarray	Model in degrees (0, 360)	required
Range			required

Returns:

Type	Description
float	Circular bias

Examples:

>>> obs = np.array([10., 355., 170.])
>>> mod = np.array([20., 5., -180.])
>>> c_bias(obs, mod)
np.float64(10.0)

Source code in modelskill/metrics.py

def c_bias(obs: np.ndarray, model: np.ndarray) -> Any:
    """Circular bias (mean error)

    Parameters
    ----------
    obs : np.ndarray
        Observation in degrees (0, 360)
    model : np.ndarray
        Model in degrees (0, 360)

    Range: [-180., 180.]; Best: 0.

    Returns
    -------
    float
        Circular bias

    Examples
    --------
    >>> obs = np.array([10., 355., 170.])
    >>> mod = np.array([20., 5., -180.])
    >>> c_bias(obs, mod)
    np.float64(10.0)
    """
    from scipy.stats import circmean

    resi = _c_residual(obs, model)
    return circmean(resi, low=-180.0, high=180.0)

c_mae

c_mae(obs, model, weights=None)

alias for circular mean absolute error

Source code in modelskill/metrics.py

def c_mae(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """alias for circular mean absolute error"""
    return c_mean_absolute_error(obs, model, weights)

c_max_error

c_max_error(obs, model)

Circular max error

Parameters:

Name	Type	Description	Default
obs	ndarray	Observation in degrees (0, 360)	required
model	ndarray	Model in degrees (0, 360)	required
Range			required

Returns:

Type	Description
float	Circular max error

Examples:

>>> obs = np.array([10., 350., 10.])
>>> mod = np.array([20., 10., 350.])
>>> c_max_error(obs, mod)
np.float64(20.0)

Source code in modelskill/metrics.py

def c_max_error(obs: np.ndarray, model: np.ndarray) -> Any:
    """Circular max error

    Parameters
    ----------
    obs : np.ndarray
        Observation in degrees (0, 360)
    model : np.ndarray
        Model in degrees (0, 360)

    Range: :math:`[0, \\infty)`; Best: 0

    Returns
    -------
    float
        Circular max error

    Examples
    --------
    >>> obs = np.array([10., 350., 10.])
    >>> mod = np.array([20., 10., 350.])
    >>> c_max_error(obs, mod)
    np.float64(20.0)
    """

    resi = _c_residual(obs, model)

    # Compute the absolute differences and then
    # find the shortest distance between angles
    abs_diffs = np.abs(resi)
    circular_diffs = np.minimum(abs_diffs, 360 - abs_diffs)
    return np.max(circular_diffs)

c_mean_absolute_error

c_mean_absolute_error(obs, model, weights=None)

Circular mean absolute error

Parameters:

Name	Type	Description	Default
obs	ndarray	Observation in degrees (0, 360)	required
model	ndarray	Model in degrees (0, 360)	required
weights	ndarray	Weights, by default None	None
Range			required

Returns:

Type	Description
float	Circular mean absolute error

Source code in modelskill/metrics.py

def c_mean_absolute_error(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """Circular mean absolute error

    Parameters
    ----------
    obs : np.ndarray
        Observation in degrees (0, 360)
    model : np.ndarray
        Model in degrees (0, 360)
    weights : np.ndarray, optional
        Weights, by default None

    Range: [0, 180]; Best: 0

    Returns
    -------
    float
        Circular mean absolute error
    """

    resi = _c_residual(obs, model)
    return np.average(np.abs(resi), weights=weights)

c_rmse

c_rmse(obs, model, weights=None)

alias for circular root mean squared error

Source code in modelskill/metrics.py

def c_rmse(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """alias for circular root mean squared error"""
    return c_root_mean_squared_error(obs, model, weights)

c_root_mean_squared_error

c_root_mean_squared_error(obs, model, weights=None)

Circular root mean squared error

Parameters:

Name	Type	Description	Default
obs	ndarray	Observation in degrees (0, 360)	required
model	ndarray	Model in degrees (0, 360)	required
weights	ndarray	Weights, by default None	None
Range			required

Returns:

Type	Description
float	Circular root mean squared error

Source code in modelskill/metrics.py

def c_root_mean_squared_error(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """Circular root mean squared error

    Parameters
    ----------
    obs : np.ndarray
        Observation in degrees (0, 360)
    model : np.ndarray
        Model in degrees (0, 360)
    weights : np.ndarray, optional
        Weights, by default None

    Range: [0, 180]; Best: 0

    Returns
    -------
    float
        Circular root mean squared error
    """
    residual = _c_residual(obs, model)
    return np.sqrt(np.average(residual**2, weights=weights))

c_unbiased_root_mean_squared_error

c_unbiased_root_mean_squared_error(obs, model, weights=None)

Circular unbiased root mean squared error

Parameters:

Name	Type	Description	Default
obs	ndarray	Observation in degrees (0, 360)	required
model	ndarray	Model in degrees (0, 360)	required
weights	ndarray	Weights, by default None	None
Range			required

Returns:

Type	Description
float	Circular unbiased root mean squared error

Source code in modelskill/metrics.py

def c_unbiased_root_mean_squared_error(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """Circular unbiased root mean squared error

    Parameters
    ----------
    obs : np.ndarray
        Observation in degrees (0, 360)
    model : np.ndarray
        Model in degrees (0, 360)
    weights : np.ndarray, optional
        Weights, by default None

    Range: [0, 180]; Best: 0

    Returns
    -------
    float
        Circular unbiased root mean squared error
    """
    from scipy.stats import circmean

    residual = _c_residual(obs, model)
    residual = residual - circmean(residual, low=-180.0, high=180.0)
    return np.sqrt(np.average(residual**2, weights=weights))

c_urmse

c_urmse(obs, model, weights=None)

alias for circular unbiased root mean squared error

Source code in modelskill/metrics.py

def c_urmse(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
) -> Any:
    """alias for circular unbiased root mean squared error"""
    return c_unbiased_root_mean_squared_error(obs, model, weights)

cc

cc(obs, model, weights=None)

alias for corrcoef

Source code in modelskill/metrics.py

def cc(obs: np.ndarray, model: np.ndarray, weights=None) -> Any:
    """alias for corrcoef"""
    return corrcoef(obs, model, weights)

corrcoef

corrcoef(obs, model, weights=None)

Pearson’s Correlation coefficient (CC)

\[ CC = \frac{\sum_{i=1}^n (model_i - \overline{model})(obs_i - \overline{obs}) } {\sqrt{\sum_{i=1}^n (model_i - \overline{model})^2} \sqrt{\sum_{i=1}^n (obs_i - \overline{obs})^2} } \]

Range: [-1, 1]; Best: 1

See Also

spearmanr np.corrcoef

Source code in modelskill/metrics.py

def corrcoef(obs, model, weights=None) -> Any:
    r"""Pearson’s Correlation coefficient (CC)

    $$
    CC = \frac{\sum_{i=1}^n (model_i - \overline{model})(obs_i - \overline{obs}) }
                   {\sqrt{\sum_{i=1}^n (model_i - \overline{model})^2}
                    \sqrt{\sum_{i=1}^n (obs_i - \overline{obs})^2} }
    $$

    Range: [-1, 1]; Best: 1

    See Also
    --------
    spearmanr
    np.corrcoef
    """
    assert obs.size == model.size
    if len(obs) <= 1:
        return np.nan

    if weights is None:
        return np.corrcoef(obs, model)[0, 1]
    else:
        C = np.cov(obs, model, fweights=weights)
        return C[0, 1] / np.sqrt(C[0, 0] * C[1, 1])

ev

ev(obs, model)

alias for explained_variance

Source code in modelskill/metrics.py

def ev(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for explained_variance"""
    assert obs.size == model.size
    return explained_variance(obs, model)

explained_variance

explained_variance(obs, model)

EV: Explained variance

EV is the explained variance and measures the proportion [0 - 1] to which the model accounts for the variation (dispersion) of the observations.

In cases with no bias, EV is equal to r2

\[ \frac{ \sum_{i=1}^n (obs_i - \overline{obs})^2 - \sum_{i=1}^n \left( (obs_i - \overline{obs}) - (model_i - \overline{model}) \right)^2}{\sum_{i=1}^n (obs_i - \overline{obs})^2} \]

Range: [0, 1]; Best: 1

See Also

r2

Source code in modelskill/metrics.py

def explained_variance(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""EV: Explained variance

     EV is the explained variance and measures the proportion
     [0 - 1] to which the model accounts for the variation
     (dispersion) of the observations.

     In cases with no bias, EV is equal to r2

    $$
    \frac{ \sum_{i=1}^n (obs_i - \overline{obs})^2 -
    \sum_{i=1}^n \left( (obs_i - \overline{obs}) -
    (model_i - \overline{model}) \right)^2}{\sum_{i=1}^n
    (obs_i - \overline{obs})^2}
    $$

    Range: [0, 1]; Best: 1

    See Also
    --------
    r2
    """

    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan

    nominator: Any = np.sum((obs - obs.mean()) ** 2) - np.sum(  # type: ignore
        ((obs - obs.mean()) - (model - model.mean())) ** 2
    )
    denominator: Any = np.sum((obs - obs.mean()) ** 2)

    return nominator / denominator

hit_ratio

hit_ratio(obs, model, a=0.1)

Fraction within obs ± acceptable deviation

\[ HR = \frac{1}{n}\sum_{i=1}^n I_{|(model_i - obs_i)|} < a \]

Range: [0, 1]; Best: 1

Examples:

>>> obs = np.array([1.0, 1.1, 1.2, 1.3, 1.4, 1.4, 1.3])
>>> model = np.array([1.02, 1.16, 1.3, 1.38, 1.49, 1.45, 1.32])
>>> hit_ratio(obs, model, a=0.05)
np.float64(0.2857142857142857)
>>> hit_ratio(obs, model, a=0.1)
np.float64(0.8571428571428571)
>>> hit_ratio(obs, model, a=0.15)
np.float64(1.0)

Source code in modelskill/metrics.py

def hit_ratio(obs: np.ndarray, model: np.ndarray, a=0.1) -> Any:
    r"""Fraction within obs ± acceptable deviation

    $$
    HR = \frac{1}{n}\sum_{i=1}^n I_{|(model_i - obs_i)|} < a
    $$

    Range: [0, 1]; Best: 1

    Examples
    --------
    >>> obs = np.array([1.0, 1.1, 1.2, 1.3, 1.4, 1.4, 1.3])
    >>> model = np.array([1.02, 1.16, 1.3, 1.38, 1.49, 1.45, 1.32])
    >>> hit_ratio(obs, model, a=0.05)
    np.float64(0.2857142857142857)
    >>> hit_ratio(obs, model, a=0.1)
    np.float64(0.8571428571428571)
    >>> hit_ratio(obs, model, a=0.15)
    np.float64(1.0)
    """
    assert obs.size == model.size

    return np.mean(np.abs(obs - model) < a)

kge

kge(obs, model)

alias for kling_gupta_efficiency

Source code in modelskill/metrics.py

def kge(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for kling_gupta_efficiency"""
    return kling_gupta_efficiency(obs, model)

kling_gupta_efficiency

kling_gupta_efficiency(obs, model)

Kling-Gupta Efficiency (KGE)

\[ KGE = 1 - \sqrt{(r-1)^2 + \left(\frac{\sigma_{mod}}{\sigma_{obs}} - 1\right)^2 + \left(\frac{\mu_{mod}}{\mu_{obs}} - 1\right)^2 } \]

where \(r\) is the pearson correlation coefficient, \(\mu_{obs},\mu_{mod}\) and \(\sigma_{obs},\sigma_{mod}\) is the mean and standard deviation of observations and model.

Range: \((-\infty, 1]\); Best: 1

References

Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F., (2009), Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377(1-2), 80-91 https://doi.org/10.1016/j.jhydrol.2009.08.003

Knoben, W. J. M., Freer, J. E., and Woods, R. A. (2019) Technical note: Inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, Hydrol. Earth Syst. Sci., 23, 4323-4331 https://doi.org/10.5194/hess-23-4323-2019

Source code in modelskill/metrics.py

def kling_gupta_efficiency(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""
    Kling-Gupta Efficiency (KGE)

    $$
    KGE = 1 - \sqrt{(r-1)^2 + \left(\frac{\sigma_{mod}}{\sigma_{obs}} - 1\right)^2 +
                                \left(\frac{\mu_{mod}}{\mu_{obs}} - 1\right)^2 }
    $$

    where $r$ is the pearson correlation coefficient, $\mu_{obs},\mu_{mod}$ and $\sigma_{obs},\sigma_{mod}$ is the mean and standard deviation of observations and model.

    Range: $(-\infty, 1]$; Best: 1

    References
    ----------
    Gupta, H. V., Kling, H., Yilmaz, K. K. and Martinez, G. F., (2009), Decomposition of the mean squared error and NSE performance criteria: Implications for improving hydrological modelling, J. Hydrol., 377(1-2), 80-91 <https://doi.org/10.1016/j.jhydrol.2009.08.003>

    Knoben, W. J. M., Freer, J. E., and Woods, R. A. (2019) Technical note: Inherent benchmark or not? Comparing Nash–Sutcliffe and Kling–Gupta efficiency scores, Hydrol. Earth Syst. Sci., 23, 4323-4331 <https://doi.org/10.5194/hess-23-4323-2019>
    """
    assert obs.size == model.size

    if len(obs) == 0 or obs.std() == 0.0:
        return np.nan

    if model.std() > 1e-12:
        r = corrcoef(obs, model)
        if np.isnan(r):
            r = 0.0
    else:
        r = 0.0

    res = 1 - np.sqrt(
        (r - 1) ** 2
        + (model.std() / obs.std() - 1.0) ** 2
        + (model.mean() / obs.mean() - 1.0) ** 2
    )

    return res

lin_slope

lin_slope(obs, model, reg_method='ols')

Slope of the regression line.

\[ slope = \frac{\sum_{i=1}^n (model_i - \overline {model})(obs_i - \overline {obs})} {\sum_{i=1}^n (obs_i - \overline {obs})^2} \]

Range: \((-\infty, \infty )\); Best: 1

Source code in modelskill/metrics.py

def lin_slope(obs: np.ndarray, model: np.ndarray, reg_method="ols") -> Any:
    r"""Slope of the regression line.

    $$
    slope = \frac{\sum_{i=1}^n (model_i - \overline {model})(obs_i - \overline {obs})}
                    {\sum_{i=1}^n (obs_i - \overline {obs})^2}
    $$

    Range: $(-\infty, \infty )$; Best: 1
    """
    assert obs.size == model.size
    return _linear_regression(obs, model, reg_method)[0]

mae

mae(obs, model, weights=None)

alias for mean_absolute_error

Source code in modelskill/metrics.py

def mae(
    obs: np.ndarray, model: np.ndarray, weights: Optional[np.ndarray] = None
) -> Any:
    """alias for mean_absolute_error"""
    assert obs.size == model.size
    return mean_absolute_error(obs, model, weights)

mape

mape(obs, model)

alias for mean_absolute_percentage_error

Source code in modelskill/metrics.py

def mape(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for mean_absolute_percentage_error"""
    return mean_absolute_percentage_error(obs, model)

max_error

max_error(obs, model)

Max (absolute) error

\[ max_{error} = max(|model_i - obs_i|) \]

Range: \([0, \infty)\); Best: 0

Source code in modelskill/metrics.py

def max_error(obs, model) -> Any:
    r"""Max (absolute) error

    $$
    max_{error} = max(|model_i - obs_i|)
    $$

    Range: $[0, \infty)$; Best: 0
    """

    assert obs.size == model.size
    return np.max(np.abs(model - obs))

mean_absolute_error

mean_absolute_error(obs, model, weights=None)

Mean Absolute Error (MAE)

\[ MAE=\frac{1}{n}\sum_{i=1}^n|model_i - obs_i| \]

Range: \([0, \infty)\); Best: 0

Source code in modelskill/metrics.py

def mean_absolute_error(
    obs: np.ndarray, model: np.ndarray, weights: Optional[np.ndarray] = None
) -> Any:
    r"""Mean Absolute Error (MAE)

    $$
    MAE=\frac{1}{n}\sum_{i=1}^n|model_i - obs_i|
    $$

    Range: $[0, \infty)$; Best: 0
    """
    assert obs.size == model.size

    error = np.average(np.abs(model - obs), weights=weights)

    return error

mean_absolute_percentage_error

mean_absolute_percentage_error(obs, model)

Mean Absolute Percentage Error (MAPE)

\[ MAPE=\frac{1}{n}\sum_{i=1}^n\frac{|model_i - obs_i|}{obs_i}*100 \]

Range: \([0, \infty)\); Best: 0

Source code in modelskill/metrics.py

def mean_absolute_percentage_error(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Mean Absolute Percentage Error (MAPE)

    $$
    MAPE=\frac{1}{n}\sum_{i=1}^n\frac{|model_i - obs_i|}{obs_i}*100
    $$

    Range: $[0, \infty)$; Best: 0
    """

    assert obs.size == model.size

    if len(obs) == 0:
        return np.nan
    if np.any(obs == 0.0):
        warnings.warn("Observation is zero, consider to use another metric than MAPE")
        return np.nan  # TODO is it better to return a large value +inf than NaN?

    return np.mean(np.abs((obs - model) / obs)) * 100

mef

mef(obs, model)

alias for model_efficiency_factor

Source code in modelskill/metrics.py

def mef(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for model_efficiency_factor"""
    return model_efficiency_factor(obs, model)

metric_has_units

metric_has_units(metric)

Check if a metric has units (dimension).

Some metrics are dimensionless, others have the same dimension as the observations.

Parameters:

Name	Type	Description	Default
metric	str or callable	Metric name or function	required

Returns:

Type	Description
bool	True if metric has a dimension, False otherwise

Examples:

>>> metric_has_units("rmse")
True
>>> metric_has_units("kge")
False

Source code in modelskill/metrics.py

def metric_has_units(metric: Union[str, Callable]) -> bool:
    """Check if a metric has units (dimension).

    Some metrics are dimensionless, others have the same dimension as the observations.

    Parameters
    ----------
    metric : str or callable
        Metric name or function

    Returns
    -------
    bool
        True if metric has a dimension, False otherwise

    Examples
    --------
    >>> metric_has_units("rmse")
    True
    >>> metric_has_units("kge")
    False
    """
    if hasattr(metric, "__name__"):
        name = metric.__name__
    else:
        name = metric

    if name not in defined_metrics:
        raise ValueError(f"Metric {name} not defined. Choose from {defined_metrics}")

    return name in METRICS_WITH_DIMENSION

model_efficiency_factor

model_efficiency_factor(obs, model)

Model Efficiency Factor (MEF)

Scale independent RMSE, standardized by Stdev of observations

\[ MEF = \frac{RMSE}{STDEV}=\frac{\sqrt{\frac{1}{n} \sum_{i=1}^n(model_i - obs_i)^2}} {\sqrt{\frac{1}{n} \sum_{i=1}^n(obs_i - \overline{obs})^2}}=\sqrt{1-NSE} \]

Range: \([0, \infty)\); Best: 0

See Also

nash_sutcliffe_efficiency root_mean_squared_error

Source code in modelskill/metrics.py

def model_efficiency_factor(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Model Efficiency Factor (MEF)

    Scale independent RMSE, standardized by Stdev of observations

    $$
    MEF = \frac{RMSE}{STDEV}=\frac{\sqrt{\frac{1}{n} \sum_{i=1}^n(model_i - obs_i)^2}}
                                    {\sqrt{\frac{1}{n} \sum_{i=1}^n(obs_i - \overline{obs})^2}}=\sqrt{1-NSE}
    $$

    Range: $[0, \infty)$; Best: 0

    See Also
    --------
    nash_sutcliffe_efficiency
    root_mean_squared_error

    """
    assert obs.size == model.size

    return rmse(obs, model) / obs.std()

nash_sutcliffe_efficiency

nash_sutcliffe_efficiency(obs, model)

Nash-Sutcliffe Efficiency (NSE)

\[ NSE = 1 - \frac {\sum _{i=1}^{n}\left(model_{i} - obs_{i}\right)^{2}} {\sum_{i=1}^{n}\left(obs_{i} - {\overline{obs}}\right)^{2}} \]

Range: \((-\infty, 1]\); Best: 1

Note

r2 = nash_sutcliffe_efficiency(nse)

References

Nash, J. E.; Sutcliffe, J. V. (1970). "River flow forecasting through conceptual models part I — A discussion of principles". Journal of Hydrology. 10 (3): 282–290. https://doi.org/10.1016/0022-1694(70)90255-6

Source code in modelskill/metrics.py

def nash_sutcliffe_efficiency(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Nash-Sutcliffe Efficiency (NSE)

    $$
    NSE = 1 - \frac {\sum _{i=1}^{n}\left(model_{i} - obs_{i}\right)^{2}}
                    {\sum_{i=1}^{n}\left(obs_{i} - {\overline{obs}}\right)^{2}}
    $$

    Range: $(-\infty, 1]$; Best: 1

    Note
    ----
    r2 = nash_sutcliffe_efficiency(nse)

    References
    ----------
    Nash, J. E.; Sutcliffe, J. V. (1970). "River flow forecasting through conceptual models part I — A discussion of principles". Journal of Hydrology. 10 (3): 282–290. <https://doi.org/10.1016/0022-1694(70)90255-6>
    """
    assert obs.size == model.size

    if len(obs) == 0:
        return np.nan
    error = 1 - (np.sum((obs - model) ** 2) / np.sum((obs - np.mean(obs)) ** 2))  # type: ignore

    return error

nse

nse(obs, model)

alias for nash_sutcliffe_efficiency

Source code in modelskill/metrics.py

def nse(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for nash_sutcliffe_efficiency"""
    return nash_sutcliffe_efficiency(obs, model)

peak_ratio

peak_ratio(obs, model, inter_event_level=0.7, AAP=2, inter_event_time='36h')

Peak Ratio

PR is the mean of the largest-N individual ratios of identified peaks in the model / identified peaks in the measurements (N number of events defined by AAP). PR is calculated only for the joint-events, ie, events that ocurr simulateneously within a window +/- 0.5*inter_event_time.

Parameters:

Name	Type	Description	Default
inter_event_level	float	Inter-event level threshold (default: 0.7).	0.7
AAP	Union[int, float]	Average Annual Peaks (ie, Number of peaks per year, on average). (default: 2)	2
inter_event_time	str	Maximum time interval between peaks (default: 36 hours).	'36h'
Range			required

Source code in modelskill/metrics.py

def peak_ratio(
    obs: pd.Series,
    model: np.ndarray,
    inter_event_level: float = 0.7,
    AAP: Union[int, float] = 2,
    inter_event_time: str = "36h",
) -> Any:
    r"""Peak Ratio

    PR is the mean of the largest-N individual ratios of identified peaks in the
    model / identified peaks in the measurements (N number of events defined by AAP). PR is calculated only for the joint-events,
    ie, events that ocurr simulateneously within a window +/- 0.5*inter_event_time.

    Parameters
    ----------
    inter_event_level (float, optional)
        Inter-event level threshold (default: 0.7).
    AAP (int or float, optional)
        Average Annual Peaks (ie, Number of peaks per year, on average). (default: 2)
    inter_event_time (str, optional)
            Maximum time interval between peaks (default: 36 hours).

    $$
    \frac{\sum_{i=1}^{N_{joint-peaks}} (\frac{Peak_{model_i}}{Peak_{obs_i}} )}{N_{joint-peaks}}
    $$

    Range: $[0, \infty)$; Best: 1.0
    """

    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan
    assert isinstance(obs.index, pd.DatetimeIndex)
    time = obs.index

    # Calculate number of years
    dt_int = time[1:].values - time[0:-1].values
    dt_int_mode = float(stats.mode(dt_int, keepdims=False)[0]) / 1e9  # in seconds
    N_years = dt_int_mode / 24 / 3600 / 365.25 * len(time)
    peak_index, AAP_ = _partial_duration_series(
        time,
        obs,
        inter_event_level=inter_event_level,
        AAP=AAP,
        inter_event_time=inter_event_time,
    )
    peaks = obs[peak_index]
    found_peaks_obs = peaks.sort_values(ascending=False)

    peak_index, _ = _partial_duration_series(
        time,
        model,
        inter_event_level=inter_event_level,
        AAP=AAP,
        inter_event_time=inter_event_time,
    )
    peaks = model[peak_index]
    found_peaks_mod = peaks.sort_values(ascending=False)

    top_n_peaks = max(1, min(round(AAP_ * N_years), np.sum(peaks)))
    # Resample~ish, find peaks spread maximum Half the inter event time (if inter event =36, select data paired +/- 18h) (or inter_event) and then select
    indices_mod = (
        abs(found_peaks_obs.index.values[:, None] - found_peaks_mod.index.values)
        < pd.Timedelta(inter_event_time) / 2
    ).any(axis=0)
    indices_obs = (
        abs(found_peaks_mod.index.values[:, None] - found_peaks_obs.index.values)
        < pd.Timedelta(inter_event_time) / 2
    ).any(axis=0)
    # Find intersection (co-existing peaks, still a large number, O(1000s))
    obs_joint = found_peaks_obs.loc[indices_obs]
    mod_joint = found_peaks_mod.loc[indices_mod]
    # Now we forget about time index, as peaks have been paired already.
    df_filter = pd.DataFrame(
        data={
            "model": mod_joint.sort_index().values,
            "observation": obs_joint.sort_index().values,
        }
    )
    df_filter["Maximum"] = df_filter.max(axis=1)
    df_filter.sort_values(by="Maximum", ascending=False, inplace=True)
    # Finally we do the selection of the N- largest peaks from either model or measured
    df_filter = df_filter.iloc[0:top_n_peaks, :]
    # Rename to avoid further refactoring
    obs_joint = df_filter.loc[:, "observation"]
    mod_joint = df_filter.loc[:, "model"]

    if len(obs_joint) == 0 or len(mod_joint) == 0:
        return np.nan
    res = np.mean(mod_joint.values / obs_joint.values)
    return res

pr

pr(obs, model, inter_event_level=0.7, AAP=2, inter_event_time='36h')

alias for peak_ratio

Source code in modelskill/metrics.py

def pr(
    obs: pd.Series,
    model: np.ndarray,
    inter_event_level: float = 0.7,
    AAP: Union[int, float] = 2,
    inter_event_time: str = "36h",
) -> Any:
    """alias for peak_ratio"""
    assert obs.size == model.size
    return peak_ratio(obs, model, inter_event_level, AAP, inter_event_time)

r2

r2(obs, model)

Coefficient of determination (R2)

Pronounced 'R-squared'; the proportion of the variation in the dependent variable that is predictable from the independent variable(s), i.e. the proportion of explained variance.

\[ R^2 = 1 - \frac{\sum_{i=1}^n (model_i - obs_i)^2} {\sum_{i=1}^n (obs_i - \overline {obs})^2} \]

Range: \((-\infty, 1]\); Best: 1

Note

r2 = nash_sutcliffe_efficiency(nse)

Examples:

>>> obs = np.array([1.0,1.1,1.2,1.3,1.4])
>>> model = np.array([1.09, 1.16, 1.3 , 1.38, 1.49])
>>> r2(obs,model)
np.float64(0.6379999999999998)

Source code in modelskill/metrics.py

def r2(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Coefficient of determination (R2)

    Pronounced 'R-squared'; the proportion of the variation in the dependent variable that is predictable from the independent variable(s), i.e. the proportion of explained variance.

    $$
    R^2 = 1 - \frac{\sum_{i=1}^n (model_i - obs_i)^2}
                    {\sum_{i=1}^n (obs_i - \overline {obs})^2}
    $$

    Range: $(-\infty, 1]$; Best: 1

    Note
    ----
    r2 = nash_sutcliffe_efficiency(nse)

    Examples
    --------
    >>> obs = np.array([1.0,1.1,1.2,1.3,1.4])
    >>> model = np.array([1.09, 1.16, 1.3 , 1.38, 1.49])
    >>> r2(obs,model)
    np.float64(0.6379999999999998)
    """
    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan

    residual = model - obs
    SSr: Any = np.sum(residual**2)
    SSt: Any = np.sum((obs - obs.mean()) ** 2)

    return 1 - SSr / SSt

rho

rho(obs, model)

alias for spearmanr

Source code in modelskill/metrics.py

def rho(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for spearmanr"""
    return spearmanr(obs, model)

rmse

rmse(obs, model, weights=None, unbiased=False)

alias for root_mean_squared_error

Source code in modelskill/metrics.py

def rmse(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
    unbiased: bool = False,
) -> Any:
    """alias for root_mean_squared_error"""
    return root_mean_squared_error(obs, model, weights, unbiased)

root_mean_squared_error

root_mean_squared_error(obs, model, weights=None, unbiased=False)

Root Mean Squared Error (RMSE)

\[ res_i = model_i - obs_i \]

\[ RMSE=\sqrt{\frac{1}{n} \sum_{i=1}^n res_i^2} \]

Unbiased version:

\[ res_{u,i} = res_i - \overline {res} \]

\[ uRMSE=\sqrt{\frac{1}{n} \sum_{i=1}^n res_{u,i}^2} \]

Range: \([0, \infty)\); Best: 0

Source code in modelskill/metrics.py

def root_mean_squared_error(
    obs: np.ndarray,
    model: np.ndarray,
    weights: Optional[np.ndarray] = None,
    unbiased: bool = False,
) -> Any:
    r"""Root Mean Squared Error (RMSE)

    $$
    res_i = model_i - obs_i
    $$

    $$
    RMSE=\sqrt{\frac{1}{n} \sum_{i=1}^n res_i^2}
    $$

    Unbiased version:

    $$
    res_{u,i} = res_i - \overline {res}
    $$

    $$
    uRMSE=\sqrt{\frac{1}{n} \sum_{i=1}^n res_{u,i}^2}
    $$

    Range: $[0, \infty)$; Best: 0

    """
    assert obs.size == model.size

    residual = obs - model
    if unbiased:
        residual = residual - residual.mean()
    error = np.sqrt(np.average(residual**2, weights=weights))

    return error

scatter_index

scatter_index(obs, model)

Scatter index (SI)

Which is the same as the unbiased-RMSE normalized by the absolute mean of the observations.

\[ \frac{ \sqrt{ \frac{1}{n} \sum_{i=1}^n \left( (model_i - \overline {model}) - (obs_i - \overline {obs}) \right)^2} } {\frac{1}{n} \sum_{i=1}^n | obs_i | } \]

Range: \([0, \infty)\); Best: 0

Source code in modelskill/metrics.py

def scatter_index(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Scatter index (SI)

    Which is the same as the unbiased-RMSE normalized by the absolute mean of the observations.

    $$
    \frac{ \sqrt{ \frac{1}{n} \sum_{i=1}^n \left( (model_i - \overline {model}) - (obs_i - \overline {obs}) \right)^2} }
    {\frac{1}{n} \sum_{i=1}^n | obs_i | }
    $$

    Range: $[0, \infty)$; Best: 0
    """
    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan

    residual = obs - model
    residual = residual - residual.mean()  # unbiased
    return np.sqrt(np.mean(residual**2)) / np.mean(np.abs(obs))

scatter_index2

scatter_index2(obs, model)

Alternative formulation of the scatter index (SI)

\[ \sqrt {\frac{\sum_{i=1}^n \left( (model_i - \overline {model}) - (obs_i - \overline {obs}) \right)^2} {\sum_{i=1}^n obs_i^2}} \]

Range: [0, 100]; Best: 0

Source code in modelskill/metrics.py

def scatter_index2(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Alternative formulation of the scatter index (SI)

    $$
    \sqrt {\frac{\sum_{i=1}^n \left( (model_i - \overline {model}) - (obs_i - \overline {obs}) \right)^2}
    {\sum_{i=1}^n obs_i^2}}
    $$

    Range: [0, 100]; Best: 0
    """
    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan

    return np.sqrt(
        np.sum(((model - model.mean()) - (obs - obs.mean())) ** 2) / np.sum(obs**2)  # type: ignore
    )

si

si(obs, model)

alias for scatter_index

Source code in modelskill/metrics.py

def si(obs: np.ndarray, model: np.ndarray) -> Any:
    """alias for scatter_index"""
    return scatter_index(obs, model)

spearmanr

spearmanr(obs, model)

Spearman rank correlation coefficient

The rank correlation coefficient is similar to the Pearson correlation coefficient but applied to ranked quantities and is useful to quantify a monotonous relationship

\[ \rho = \frac{\sum_{i=1}^n (rmodel_i - \overline{rmodel})(robs_i - \overline{robs}) } {\sqrt{\sum_{i=1}^n (rmodel_i - \overline{rmodel})^2} \sqrt{\sum_{i=1}^n (robs_i - \overline{robs})^2} } \]

Range: [-1, 1]; Best: 1

Examples:

>>> obs = np.linspace(-20, 20, 100)
>>> mod = np.tanh(obs)
>>> rho(obs, mod)
np.float64(0.9999759973116955)
>>> spearmanr(obs, mod)
np.float64(0.9999759973116955)

See Also

corrcoef

Source code in modelskill/metrics.py

def spearmanr(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Spearman rank correlation coefficient

    The rank correlation coefficient is similar to the Pearson correlation coefficient but
    applied to ranked quantities and is useful to quantify a monotonous relationship

    $$
    \rho = \frac{\sum_{i=1}^n (rmodel_i - \overline{rmodel})(robs_i - \overline{robs}) }
                    {\sqrt{\sum_{i=1}^n (rmodel_i - \overline{rmodel})^2}
                    \sqrt{\sum_{i=1}^n (robs_i - \overline{robs})^2} }
    $$

    Range: [-1, 1]; Best: 1

    Examples
    --------
    >>> obs = np.linspace(-20, 20, 100)
    >>> mod = np.tanh(obs)
    >>> rho(obs, mod)
    np.float64(0.9999759973116955)
    >>> spearmanr(obs, mod)
    np.float64(0.9999759973116955)

    See Also
    --------
    corrcoef
    """
    import scipy.stats

    return scipy.stats.spearmanr(obs, model)[0]

urmse

urmse(obs, model, weights=None)

Unbiased Root Mean Squared Error (uRMSE)

\[ res_i = model_i - obs_i \]

\[ res_{u,i} = res_i - \overline {res} \]

\[ uRMSE = \sqrt{\frac{1}{n} \sum_{i=1}^n res_{u,i}^2} \]

Range: \([0, \infty)\); Best: 0

See Also

root_mean_squared_error

Source code in modelskill/metrics.py

def urmse(
    obs: np.ndarray, model: np.ndarray, weights: Optional[np.ndarray] = None
) -> Any:
    r"""Unbiased Root Mean Squared Error (uRMSE)

    $$
    res_i = model_i - obs_i
    $$

    $$
    res_{u,i} = res_i - \overline {res}
    $$

    $$
    uRMSE = \sqrt{\frac{1}{n} \sum_{i=1}^n res_{u,i}^2}
    $$

    Range: $[0, \infty)$; Best: 0

    See Also
    --------
    root_mean_squared_error
    """
    return root_mean_squared_error(obs, model, weights, unbiased=True)

willmott

willmott(obs, model)

Willmott's Index of Agreement

A scaled representation of the predictive accuracy of the model against observations. A value of 1 indicates a perfect match, and 0 indicates no agreement at all.

\[ willmott = 1 - \frac{\frac{1}{n} \sum_{i=1}^n(model_i - obs_i)^2} {\frac{1}{n} \sum_{i=1}^n(|model_i - \overline{obs}| + |obs_i - \overline{obs}|)^2} \]

Range: [0, 1]; Best: 1

Examples:

>>> obs = np.array([1.0, 1.1, 1.2, 1.3, 1.4, 1.4, 1.3])
>>> model = np.array([1.02, 1.16, 1.3, 1.38, 1.49, 1.45, 1.32])
>>> willmott(obs, model)
np.float64(0.9501403174479723)

References

Willmott, C. J. 1981. "On the validation of models". Physical Geography, 2, 184–194.

Source code in modelskill/metrics.py

def willmott(obs: np.ndarray, model: np.ndarray) -> Any:
    r"""Willmott's Index of Agreement

    A scaled representation of the predictive accuracy of the model against observations. A value of 1 indicates a perfect match, and 0 indicates no agreement at all.

    $$
    willmott = 1 - \frac{\frac{1}{n} \sum_{i=1}^n(model_i - obs_i)^2}
                        {\frac{1}{n} \sum_{i=1}^n(|model_i - \overline{obs}| + |obs_i - \overline{obs}|)^2}
    $$

    Range: [0, 1]; Best: 1

    Examples
    --------
    >>> obs = np.array([1.0, 1.1, 1.2, 1.3, 1.4, 1.4, 1.3])
    >>> model = np.array([1.02, 1.16, 1.3, 1.38, 1.49, 1.45, 1.32])
    >>> willmott(obs, model)
    np.float64(0.9501403174479723)

    References
    ----------
    Willmott, C. J. 1981. "On the validation of models". Physical Geography, 2, 184–194.
    """

    assert obs.size == model.size
    if len(obs) == 0:
        return np.nan

    residual = model - obs
    nominator: Any = np.sum(residual**2)
    denominator: Any = np.sum(
        (np.abs(model - obs.mean()) + np.abs(obs - obs.mean())) ** 2
    )

    return 1 - nominator / denominator