Vs vd

README.md

@@ -82,6 +82,8 @@ Returns a dictionary containing the calculated metrics and their corresponding v
   * `madp`: Mean Absolute Deviation of percentiles
   * `madc`: `mad + madp`
   * `kge`: Kling–Gupta Efficiency
+  * `vs`: Variance Similarity (Koh et al., 2012), unitless, ranges from 0 to 1 (1 = equal variances)
+  * `vd`: Variance Dissimilarity, complement of `vs` (`1 - vs`), ranges from -1 to 1 (-1 = noisy obs & model flat, 0 = equal variances, 1 = noisy model & obs flat)
 * [The storm metrics](#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).
   * `R1`: Difference between observed and modelled for the biggest storm
   * `R1_abs`: Absolute difference between observed and modelled for the biggest storm

pyproject.toml

@@ -1,6 +1,6 @@
 [tool.poetry]
 name = "seastats"
-version = "0.1.1"
+version = "0.1.2"
 description = "package for metocean statistics"
 authors = [
   "tomsail <saillour.thomas@gmail.com>",

seastats/__init__.py

@@ -23,6 +23,8 @@
 from seastats.stats import get_rmse
 from seastats.stats import get_slope_intercept
 from seastats.stats import get_slope_intercept_pp
+from seastats.stats import get_vd
+from seastats.stats import get_vs
 from seastats.storms import match_extremes
 
 __version__ = importlib.metadata.version(__name__)
@@ -46,6 +48,8 @@
     "get_slope_intercept",
     "get_slope_intercept_pp",
     "get_stats",
+    "get_vd",
+    "get_vs",
     "match_extremes",
     "STORM_METRICS",
     "STORM_METRICS_ALL",
@@ -76,6 +80,8 @@
     "madp",
     "madc",
     "kge",
+    "vs",
+    "vd",
 ]
 GENERAL_METRICS = ["bias", "rms", "rmse", "cr", "nse", "kge"]
 STORM_METRICS = ["R1", "R3", "error"]
@@ -137,6 +143,8 @@ def get_stats(  # noqa: C901
     - `madp`: The median absolute deviation of the simulated time series data from its median, calculated using the percentiles of the observed time series data.
     - `madc`: The median absolute deviation of the simulated time series data from its median, calculated by adding `mad` to `madp`
     - `kge`: The Kling-Gupta efficiency between the simulated and observed time series data.
+    - `vs`: Variance Similarity (Koh et al., 2012), unitless, ranges from 0 to 1 (1 = equal variances).
+    - `vd`: Variance Dissimilarity, the complement of `vs` (1 - vs), ranges from -1 to 1 (-1 = noisy obs & model flat, 0 = equal variances, 1 = noisy model & obs flat).
     - `R1`: Difference between observed and modelled for the biggest storm
     - `R1_abs`: Absolute R1 (R1 divided by observed value)
     - `R1_abs_norm`: Absolute normalized R1 (R1_abs divided by observed max peak)
@@ -209,6 +217,10 @@ def get_stats(  # noqa: C901
                 stats["madc"] = get_madc(sim, obs)
             case "kge":
                 stats["kge"] = get_kge(sim, obs)
+            case "vs":
+                stats["vs"] = get_vs(sim, obs)
+            case "vd":
+                stats["vd"] = get_vd(sim, obs)
             case "R1":
                 stats["R1"] = extreme_df["error"].iloc[0]
             case "R1_abs":

seastats/stats.py

@@ -47,6 +47,21 @@ def get_corr(sim: pd.Series[float], obs: pd.Series[float]) -> float:
     return float(sim.corr(obs))
 
 
+def get_vs(sim: pd.Series[float], obs: pd.Series[float]) -> float:
+    std_sim = sim.std()
+    std_obs = obs.std()
+    denominator = 0.5 * (std_sim**2 + std_obs**2)
+    if denominator == 0:
+        return float("nan")
+    return float((std_sim * std_obs) / denominator)
+
+
+def get_vd(sim: pd.Series[float], obs: pd.Series[float]) -> float:
+    std_sim = sim.std()
+    std_obs = obs.std()
+    return np.sign(std_sim - std_obs) * (1 - get_vs(sim, obs))
+
+
 def get_nse(sim: pd.Series[float], obs: pd.Series[float]) -> float:
     nse = 1 - np.nansum(np.subtract(obs, sim) ** 2) / np.nansum((obs - float(np.nanmean(obs))) ** 2)
     return float(nse)

tests/compute_stats_test.py

@@ -25,6 +25,8 @@
     "madp": 0.213,
     "madc": 0.266,
     "kge": 0.810,
+    "vs": 1.000,
+    "vd": 0.000,
 }
 EXPECTED_099 = {
     "R1": -0.292,
@@ -56,6 +58,8 @@
     "sim_std": 0.053,
     "slope": 0.043,
     "slope_pp": 0.626,
+    "vs": 0.884,
+    "vd": -0.116,
 }