Skip to content

Latest commit

 

History

History
142 lines (123 loc) · 7.69 KB

README.md

File metadata and controls

142 lines (123 loc) · 7.69 KB
Raw

SeaStats

seastats is a simple package to compare and analyse 2 time series. We use the following convention in this repo:

  • sim: modelled surge time series
  • mod: observed surge time series

The main function is:

def get_stats(
    sim: Series,
    obs: Series,
    metrics: Sequence[str] = SUGGESTED_METRICS,
    quantile: float = 0,
    cluster: int = 72,
    round: int = -1
) -> dict[str, float]

Calculates various statistical metrics between the simulated and observed time series data.

Parameters:

  • sim (pd.Series). The simulated time series data.
  • obs (pd.Series). The observed time series data.
  • metrics (list[str]). (Optional) The list of statistical metrics to calculate. If metrics = ["all"], all items in SUPPORTED_METRICS will be calculated. Default is all items in SUGGESTED_METRICS.
  • quantile (float). (Optional) Quantile used to calculate the metrics. Default is 0 (no selection)
  • cluster (int). (Optional) Cluster duration for grouping storm events. Default is 72 hours.
  • round (int). (Optional) Apply rounding to the results to. Default is no rounding (value is -1)

Returns a dictionary containing the calculated metrics and their corresponding values. With 2 types of metrics:

  • The "general" metrics: All the basic metrics needed for signal comparison (RMSE, RMS, Correlation etc..). See details below
    • bias: Bias
    • rmse: Root Mean Square Error
    • rms: Root Mean Square
    • rms_95: Root Mean Square for data points above 95th percentile
    • sim_mean: Mean of simulated values
    • obs_mean: Mean of observed values
    • sim_std: Standard deviation of simulated values
    • obs_std: Standard deviation of observed values
    • mae: Mean Absolute Error
    • mse: Mean Square Error
    • nse: Nash-Sutcliffe Efficiency
    • lamba: Lambda index
    • cr: Pearson Correlation coefficient
    • cr_95: Pearson Correlation coefficient for data points above 95th percentile
    • slope: Slope of Model/Obs correlation
    • intercept: Intercept of Model/Obs correlation
    • slope_pp: Slope of Model/Obs correlation of percentiles
    • intercept_pp: Intercept of Model/Obs correlation of percentiles
    • mad: Mean Absolute Deviation
    • madp: Mean Absolute Deviation of percentiles
    • madc: mad + madp
    • kge: Kling–Gupta Efficiency
  • The storm metrics: a PoT selection is done on the observed signal (using the match_extremes() function). Function returns the decreasing extreme event peak values for observed and modeled signals (and time lag between events). See details below.
    • R1: Difference between observed and modelled for the biggest storm
    • R1_norm: Normalized R1 (R1 divided by observed value)
    • R3: Average difference between observed and modelled for the three biggest storms
    • R3_norm: Normalized R3 (R3 divided by observed value)
    • error: Average difference between observed and modelled for all storms
    • error_norm: Normalized error (error divided by observed value)

General metrics

A. Dimensional Statistics:

Mean Error (or Bias)

$$\langle x_c - x_m \rangle = \langle x_c \rangle - \langle x_m \rangle$$

RMSE (Root Mean Squared Error)

$$\sqrt{\langle(x_c - x_m)^2\rangle}$$

Mean-Absolute Error (MAE):

$$\langle |x_c - x_m| \rangle$$

B. Dimentionless Statistics (best closer to 1)

Performance Scores (PS) or Nash-Sutcliffe Eff (NSE): $$1 - \frac{\langle (x_c - x_m)^2 \rangle}{\langle (x_m - x_R)^2 \rangle}$$

Correlation Coefficient (R):

$$\frac {\langle x_{m}x_{c}\rangle -\langle x_{m}\rangle \langle x_{c}\rangle }{{\sqrt {\langle x_{m}^{2}\rangle -\langle x_{m}\rangle ^{2}}}{\sqrt {\langle x_{c}^{2}\rangle -\langle x_{c}\rangle ^{2}}}}$$

Kling–Gupta Efficiency (KGE):

$$1 - \sqrt{(r-1)^2 + b^2 + (g-1)^2}$$ with :

  • r the correlation
  • b the modified bias term (see ref) $$\frac{\langle x_c \rangle - \langle x_m \rangle}{\sigma_m}$$
  • g the std dev term $$\frac{\sigma_c}{\sigma_m}$$

Lambda index ($\lambda$), values closer to 1 indicate better agreement:

$$\lambda = 1 - \frac{\sum{(x_c - x_m)^2}}{\sum{(x_m - \overline{x}_m)^2} + \sum{(x_c - \overline{x}_c)^2} + n(\overline{x}_m - \overline{x}_c)^2 + \kappa}$$

  • with kappa $$2 \cdot \left| \sum{((x_m - \overline{x}_m) \cdot (x_c - \overline{x}_c))} \right|$$

Storm metrics

The functions uses the match_extremes() function (detailed below) and returns:

  • R1: the error for the biggest storm
  • R3: the mean error for the 3 biggest storms
  • error: the mean error for all the storms above the threshold.
  • R1_norm/R3_norm/error: Same methodology, but values are in normalised (in %) relatively to the observed peaks.

case of NaNs

The storm_metrics() might return:

{'R1': np.nan,
 'R1_norm': np.nan,
 'R3': np.nan,
 'R3_norm': np.nan,
 'error': np.nan,
 'error_norm': np.nan}

Extreme events

Example of implementation:

from seastats.storms import match_extremes
extremes_df = match_extremes(sim, obs, 0.99, cluster = 72)
extremes_df

The modeled peaks are matched with the observed peaks. Function returns a pd.DataFrame of the decreasing observed storm peaks as follows:

time observed observed time observed model time model diff error error_norm tdiff
2022-01-29 19:30:00 0.803 2022-01-29 19:30:00 0.565 2022-01-29 17:00:00 -0.237 0.237 0.296 -2.5
2022-02-20 20:30:00 0.639 2022-02-20 20:30:00 0.577 2022-02-20 20:00:00 -0.062 0.062 0.0963 -0.5
...
2022-11-27 15:30:00 0.386 2022-11-27 15:30:00 0.400 2022-11-27 17:00:00 0.014 0.014 0.036 1.5

with:

  • diff the difference between modeled and observed peaks
  • error the absolute difference between modeled and observed peaks
  • tdiff the time difference between modeled and observed peaks

NB: the function uses pyextremes in the background, with PoT method, using the quantile value of the observed signal as physical threshold and passes the cluster_duration argument.

this happens when the function storms/match_extremes.py couldn't finc concomitent storms for the observed and modeled time series.

Usage

see notebook for details

get all metrics in a 3 liner:

from seastats import get_stats, GENERAL_METRICS_ALL, STORM_METRICS_ALL
general = get_stats(sim, obs, metrics = GENERAL_METRICS)
storm = get_stats(sim, obs, quantile = 0.99, metrics = STORM_METRICS) # we use a different quantile for PoT selection
pd.DataFrame(dict(general, **storm), index=['abed'])
bias rmse rms rms_95 sim_mean obs_mean sim_std obs_std nse lamba cr cr_95 slope intercept slope_pp intercept_pp mad madp madc kge R1 R1_norm R3 R3_norm error error_norm
abed -0.007 0.086 0.086 0.088 -0 0.007 0.142 0.144 0.677 0.929 0.817 0.542 0.718 -0.005 1.401 -0.028 0.052 0.213 0.265 0.81 0.237364 0.295719 0.147163 0.207019 0.0938142 0.177533