hazardreturn.py

#
# Snippet to empirically determine credible intervals for return period plots
#
# aslak grinsted 2021
#


import numpy as np
from scipy.stats import beta

def returnperiod(N, ptiles=np.array([0.05, 0.17, 0.50, 0.83, 0.95]), dt=1):
    # PAPER: https://www.jstor.org/stable/2676784?seq=1#metadata_info_tab_contents
    # section 3.1.3 "jeffreys interval"
    X = np.arange(N,0,-1)
    RCI = np.empty(shape=(N,ptiles.shape[0]))
    a1 = 0.5
    a2 = 0.5
    for ii in range(ptiles.shape[0]):
        RCI[:,ii] = 1./beta.ppf(1-ptiles[ii],X+a1,N-X+a2);
    return RCI*dt


if __name__ == "__main__":
    #code for testing....
    import matplotlib.pyplot as plt

    #make some surrogate data
    from scipy.stats import genextreme
    dist = genextreme(c=-0.5, loc=1.5, scale=2.5)
    mag = dist.rvs(size=100)

    mag = np.sort(mag)
    RCI = returnperiod(mag.shape[0])

    plt.figure()
    plt.fill_betweenx(mag,RCI[:,0],RCI[:,-1],color='#BBBBFF')
    plt.fill_betweenx(mag,RCI[:,1],RCI[:,-2],color='#9999FF')
    plt.plot(RCI[:,2],mag,'k+')
    plt.gca().set_xscale('log')
    plt.ylabel('magnitude')
    plt.xlabel('Return period')

    Ctrue = dist.cdf(x=mag)
    Rtrue = 1./(1-Ctrue)

    plt.plot(Rtrue,mag,'r',label='R_true')