Add function to compute value at a given return period

Author: mbruggs

examples/extreme_response_full_sea_state_example.ipynb

@@ -464,7 +464,8 @@
     "This object provides common statistical functions (PDF, CDF, PPF, etc.) and metrics (expected value, median, etc). \n",
     "Here, we will look at the survival function and the 100-year return level. \n",
     "\n",
-    "The value of the survival function at a given return level (e.g. 100-years) (`s_t` in the code below) is $1/N$ where $N$ is the number of short-term periods in the return period. In this case $N$ is the number of 3-hour periods in a 100-year period, which gives `s_t`$\\approx 3e-6$. The corresponding response, i.e. the 100-year wave elevation in this case, is given as the inverse cumulative function (ppf) of $1-$`s_t`. This gives a 100-year wave of about 10.4 meters.  "
+    "The value of the survival function at a given return level (e.g. 100-years) \n",
+    "may be calculated using the `return_year_value` function. This gives a 100-year wave of about 11.1 meters.  "
    ]
   },
   {
@@ -483,8 +484,7 @@
    "source": [
     "t_st_hr = t_st/(60.0*60.0)\n",
     "t_return_yr = 100.0\n",
-    "s_t = 1.0/(365.25*24*t_return_yr/t_st_hr)\n",
-    "x_t = lte.ppf(1-s_t)\n",
+    "x_t = extreme.return_year_value(lte.ppf, t_return_yr, t_st_hr)\n",
     "\n",
     "print(f\"100-year elevation: {x_t} m\")"
    ]
@@ -493,7 +493,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Finally we plot the survival function and show the 100-year return level (dashed grey lines). The 100-year value is about 104 m (where the grey line intersects the x-axis)"
+    "Finally we plot the survival function and show the 100-year return level (dashed grey lines). The 100-year value is about 11.1 m (where the grey line intersects the x-axis)"
    ]
   },
   {

mhkit/loads/extreme.py

@@ -727,3 +727,31 @@ def mler_export_time_series(rao, mler, sim, k):
     mler_ts = mler_ts.rename(columns={0: 'WaveHeight', 1: 'LinearResponse'})
 
     return mler_ts
+
+
+def return_year_value(ppf, return_year, short_term_period_hr):
+    """
+    Calculate the value from a given distribution corresponding to a particular
+    return year.
+
+    Parameters
+    ----------
+    ppf: callable function of 1 argument
+        Percentage Point Function (inverse CDF) of short term distribution.
+    return_year: int
+        Return period in years.
+    short_term_period_hr: int, float
+        Short term period the distribution is created from in hours.
+
+    Returns
+    -------
+    value: float
+        The value corresponding to the return period from the distribution.
+    """
+    assert callable(ppf)
+    assert isinstance(return_year, int)
+    assert isinstance(short_term_period_hr, (float, int))
+
+    p = 1 / (return_year * 365.25 * 24 / short_term_period_hr)
+
+    return ppf(1 - p)

mhkit/tests/loads/test_loads.py

@@ -180,6 +180,20 @@ def test_mler_export_time_series(self):
 
         assert_frame_equal(self.mler_ts, mler_ts, atol=0.0001)
 
+    def test_return_year_value(self):
+        dist = stats.norm
+        return_year = 50
+        short_term_period = 1
+
+        val = loads.extreme.return_year_value(dist.ppf, return_year, short_term_period)
+
+        # This is a little hard to test as there is not an indepedent means by
+        # which we can do this calculation. The value below has been taken from
+        # an initial calculation and will indicate if the return_year_value is
+        # changed unexpectedly.
+        want = 4.5839339
+        self.assertAlmostEqual(want, val, 5)
+
     def test_longterm_extreme(self):
         ste_1 = stats.norm
         ste_2 = stats.norm