基于python的Poisson-Gumbel 泊松耿贝尔复合极值分布实现

1. 基本介绍

1973 年，马逢时等提出了海洋工程建筑中设计波高推算的一种新的方法。该方法在论证了一个离散分布和一个连续分布可构成“复合极值分布”这一理论问题的基础上，使用 Poisson-Gumbel 复合极值分布来推算台风影响下多年一遇的设计波高。在资料年限较长的情况下，使用该新方法的计算结果与 P-Ⅲ型曲
线、耿贝尔曲线大致接近，但在资料年限较短的情况下，该方法计算结果较为稳定，因此《港口与航道水文规范》(JTS145-2015)中规定，对于台风多发海域，某一波向一年中出现一个以上较大台风波高时，可按台风波高的最大值系列取样，采用泊松-耿贝尔复合极值分布确定不同重现期的设计波高。
泊松-耿贝尔分布的有效波高可表示为：

其中

H为波高平均值，S为波高方差，σn与yn均为需要从规范中查表的参数(已集成在程序中)，N为统计数据的年数，n为N年总共发生的台风数量，p为累积频率，即重现期的倒数。

2. 程序所需文件

需要的文件格式为csv文件，其中数据格式为从大到小排列的数据，如下图：

3. 实现程序

（变量命名十分混乱，并且使用了大量的拼音，懒得修改了-_-）

############################# by YGY from TJU
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
import math
from matplotlib.font_manager import FontProperties
import matplotlib

#######################原数据计算累积频率
def leijipinlv(n):
    p = []
    for i in range(n):
        p.append((i+1)/(n+1))
    return p

######################原数据计算分位数
def fenweishu(n,p):
    F = []   #分位数
    for i in range(n):
        F.append(norm.ppf(p[i]))  #返回标准正态分布p[i]处对应的累积频率积分
    return F

#####################原数据计算横坐标
def hengzuobiao(n,F):
    H = []
    for i in range(n):
        H.append(-norm.ppf(0.0001)+F[i])   #原数据横坐标
    return H
#####################原数据计算均值和标准差
def aveandvar(hei):
    ave = np.mean(hei)
    variance = np.std(hei)
    return ave,variance

#####################计算beita
def qbeita(n,t):  #n是年数，t是总数，下面的列表均为海港水文规范中提供的数据
    p2 = [0.01,0.05,0.1,0.2,0.5,1,2,4,5,10,25,50,55,60,70,80,90,95]
    beita = []
    jtstt = [8,9,10,11,12,13,14,15,16,17,18,19,20,22,24,26,28,30,35,40,45,50,60,70,80,90,100,200,500,1000]    #规范台风总数，yn，sigma对应
    ynt = [0.4843,0.4902,0.4952,0.4996,0.5035,0.5070,0.5100,0.5128,0.5157,0.5181,0.5202,0.5220,0.5236,0.5268,0.5296,
           0.5320,0.5343,0.5362,0.5403,0.5436,0.5463,0.5485,0.5221,0.5548,0.5569,0.5586,0.5600,0.5672,0.5724,0.5745]
    sgmt = [0.9043,0.9288,0.9497,0.9676,0.9833,0.9972,1.0095,1.0206,1.0316,1.0411,1.0493,1.0566,1.0268,1.0754,1.0864,
            1.0961,1.1047,1.1124,1.1285,1.1413,1.1519,1.1607,1.1747,1.1854,1.1938,1.2007,1.2065,1.2360,1.2588,1.2685]

    ########################通过台风总数直接或插值求解ynsigama
    if t in jtstt:
        print(jtstt.index(t))
        yn = ynt[jtstt.index(t)]
        sgm = sgmt[jtstt.index(t)]
    else:
        for i in range(len(jtstt)-1):
            if jtstt[i] < t < jtstt[i+1]:
                bizhong1 = (jtstt[i+1]-t)/(jtstt[i+1]-jtstt[i])
                bizhong2 = (t-jtstt[i])/(jtstt[i+1]-jtstt[i])
                yn = bizhong1*ynt[i]+bizhong2*ynt[i+1]
                sgm = bizhong1*sgmt[i]+bizhong2*sgmt[i+1]
    print(t,yn,sgm)
    ################################求beita
    for i in range(len(p2)):
        if 1+n*math.log(1-p2[i]/100)/t >0:
            if -(math.log(1+n*math.log(1-p2[i]/100)/t)) > 0:
                beita.append(-1/sgm*(yn+math.log(-(math.log(1+n*math.log(1-p2[i]/100)/t)))))
            else:
                break
    return beita

##################################计算曲线纵坐标
def qjieguo(ave,variance,beita):
    jieguo = []
    for i in range(len(beita)):
        jieguo.append(ave+beita[i]*variance)
    return jieguo

################################计算曲线横坐标
def quxianhengzuobiao(JG):
    JH = []
    p3 = [0.01, 0.05, 0.1, 0.2, 0.5, 1, 2, 4, 5, 10, 25, 50, 55, 60, 70, 80, 90, 95]
    for i in range(len(JG)):
        JH.append(-norm.ppf(0.0001)+norm.ppf(p3[i]/100))
    return JH

##############################出图
def jieguoplot(H,hei,jieguohengzuobiao,jieguo,A,V):
    legend_font = {"family": "Times New Roman"}
    legend_font2 = {"family": "Simsun"}
    plt.rcParams['font.sans-serif'] = ['SimHei']   #解决中文不显示问题
    #matplotlib.rcParams['font.family'] = 'SimHei'
    A2 = round(A,2)   #保留两位小数
    V2 = round(V,2)
    duishuc = [0.01,0.1,0.2,0.5,1,2,5,10,20,30,40,50,60,70,80,90,95]
    duishuh =[]
    h = []
    for i in range(len(duishuc)):
        duishuh.append(-norm.ppf(0.0001)+norm.ppf(duishuc[i]/100))
        h.append(10)
    #print('duishu:',duishuh)
    plt.grid()  #显示网格
    plt.scatter(H, hei,color='k',s=15,label=u'经验累积频率')
    plt.plot(jieguohengzuobiao,jieguo,color='r',label=u'泊松耿贝尔理论分布累积频率曲线')
    plt.xlim(0.00,-norm.ppf(0.0001)+norm.ppf(duishuc[len(duishuc)-1]/100))  #x轴范围
    plt.xlabel(u'频率(%)',fontproperties='SimHei')   #x轴名称
    plt.ylabel(u'有效波高(m)',fontproperties='SimHei')   #y轴名称
    #plt.text(duishuh[11],jieguo[1],'H='+str(A2)+'S='+str(V2))
    ax = plt.gca()
    ax.set_xticks(duishuh)   #x轴刻度
    plt.xticks(duishuh,duishuc,fontproperties = 'Times New Roman')   #替换x轴坐标显示内容
    plt.yticks(fontproperties = 'Times New Roman')
    plt.legend(ncol=1)   #图上显示图例
    plt.savefig('示例.jpg')
    plt.show()

##########################主程序
total = int(input('请输入台风总数：'))
nian = int(input('请输入总年数：'))
height = np.loadtxt('波高.csv',usecols=(0))   #原始数据
resultp= leijipinlv(total)    #累积频率
resultF = fenweishu(total,resultp)   #分位数
resultH = hengzuobiao(total,resultF) #横坐标
resultA,resultV = aveandvar(height)   #均值，标准差
print(resultA,resultV)
resultB = qbeita(nian,total)   #beita
resultJ = qjieguo(resultA,resultV,resultB)   #曲线纵坐标
resultJH = quxianhengzuobiao(resultJ)        #曲线横坐标
print('resultH',len(resultH),resultH)
print('height',len(height),height)
print('resultJH',len(resultJH),resultJH)
print('resultJ',len(resultJ),resultJ)
print('100年一遇：',resultJ[5],'50年一遇：',resultJ[6],'20年一遇：',resultJ[8],'10年一遇：',resultJ[9],'4年一遇：',resultJ[10])   #重现期
jieguoplot(resultH,height,resultJH,resultJ,resultA,resultV)   #绘图

3. 程序运行及计算效果

程序运行需要输入两个参数，为台风场次（即第1步中文件的行数和统计的台风发生总年数）：

在此我们假设20年间共发生31场台风。
程序绘图如图：

注意，这里的原始数据是我们随便给的。
程序可输出不同重现期的计算结果如下：
100年一遇： 65.02338903947609 50年一遇： 57.57381724808603 20年一遇： 47.57075647818699 10年一遇： 39.73328705552583 4年一遇： 28.366663823007947