Add 'distribution_fitting.py'

slr/distribution_fitting.py

@@ -0,0 +1,80 @@
+"""Fit probability distributions to IPCC sea level rise forecasts.
+
+
+Reads:
+    'IPCC AR6.xlsx'
+
+Writes:
+    'png/*.png'
+
+D. Howe
+d.howe@wrl.unsw.edu.au
+2022-05-12
+"""
+import os
+import re
+import numpy as np
+import pandas as pd
+from scipy import stats, optimize
+import matplotlib.pyplot as plt
+
+# Read data
+df = pd.read_excel('IPCC AR6.xlsx', index_col=[0, 1, 2, 3, 4])
+df = df.sort_index()
+
+# Use all 'medium' confidence scenarios for intermediate quantiles
+scenarios = ['ssp119', 'ssp126', 'ssp245', 'ssp370', 'ssp585']
+dff = df.loc[838, 'total', 'medium', scenarios].groupby('quantile').mean()
+
+# Use ssp119/ssp585 for 5th and 95th quantiles
+dff.loc[5] = df.loc[838, 'total', 'medium', 'ssp119', 5]
+dff.loc[95] = df.loc[838, 'total', 'medium', 'ssp585', 95]
+dff = dff.T
+
+dff.index.name = 'year'
+percentiles = dff.columns.values / 100
+values = dff.loc[2150].values
+
+x_min = values.min() - 0.2
+x_max = values.max() + 0.2
+x = np.linspace(x_min, x_max, num=1000)
+
+# Get statistical distributions
+distributions = [
+    getattr(stats, d) for d in dir(stats)
+    if isinstance(getattr(stats, d), stats.rv_continuous)
+]
+
+for dist in distributions:
+
+    def cdf(x, loc, scale):
+        """Calculate cumulative density function"""
+        return dist(loc=loc, scale=scale).cdf(x)
+
+    try:
+        loc, scale = optimize.curve_fit(cdf, values, percentiles)[0]
+        p = {'loc': loc, 'scale': scale}
+
+    except TypeError:
+        continue
+
+    fig, ax = plt.subplots(1, 2, figsize=(6, 2))
+
+    ax[0].plot(x, 100 * dist.cdf(x, **p))
+    ax[0].plot(values, 100 * percentiles, '.', c='#444444')
+    ax[1].plot(x, 100 * dist.pdf(x, **p))
+
+    ax[0].axhline(y=100, c='#000000', lw=0.8, zorder=-1)
+    ax[0].set_ylim(0, 101)
+    ax[1].set_ylim(bottom=0)
+
+    ax[0].set_title(dist.name, x=-0.7, y=0.5, ha='left')
+
+    for a in ax.ravel():
+        a.spines['right'].set_visible(False)
+        a.spines['top'].set_visible(False)
+
+    ax[1].set_yticks([])
+
+    plt.savefig(f'png/{dist.name}.png', bbox_inches='tight', dpi=100)
+    plt.close()