設計 Pythonによる水文量頻度解析用自作関数とその利用

記事の最後に行く

自作関数

以下の確率分布に対するパラメータ推定と予測値計算の部分を自作関数としてプログラム本体から分離している。

含まれる関数（確率分布）

関数名	説明
ln3	３変数対数正規分布（LN3）
lp3	対数ピアソンIII型分布（LP3）
gev	一般化極値分布（GEV）
gum	グンベル分布（GUM）
sqe	平方根指数型最大値分布（SQRT-ET）

関数への入力値（引数）はいずれについても以下の３つ。

入力値（引数）

変数名	説明
xx	昇順に並び替えられた水文量観測値
pp	各観測値に相当する非超過確率（プロッティング・ポジション公式より算出）<.td>
pyear	予測値を求めたい確率年（昇順）

関数からの戻り値はいずれについても以下の４つ。

戻り値

変数名	説明
res	計算された確率分布関数のパラメータ（numpy.array）
ye	入力した非超過確率（pp）に対応する予測値（numpy.array）
yp	入力した確率年に対する予測値（numpy.array）
ss	確率分布名、算定されたパラメータなどの情報を記載した文字列

#im_hydrology.py
import numpy as np
from scipy.stats import norm  # normal distribution for LN3
from scipy.stats import gamma # Gamma distribution for LP3
import scipy.special as ssp   # Gamma function for GEV

def ln3(xx,pp,pyear):
    def func_ln3(xx):
        n=len(xx)
        mx=np.sum(xx)/n
        s1=np.sqrt(np.sum((xx-mx)**2)/n)
        cx=np.sum(((xx-mx)/s1)**3)/n
        sx=np.sqrt(n/(n-1))*s1
        #rx=np.sqrt(n*(n-1))/(n-2)*cx
        a=1.01+7.01/n+14.66/n**2
        b=1.69/n+74.66/n**2
        rx=cx*(a+b*cx**3)
        beta=1+rx**2/2
        x=(beta+np.sqrt(beta**2-1))**(1/3)+(beta-np.sqrt(beta**2-1))**(1/3)-1
        sd=np.sqrt(np.log(x))
        mu=np.log(sx/np.sqrt(x*(x-1)))
        aa=mx-np.exp(mu)*np.exp(sd**2/2)
        return aa,mu,sd
    def est_ln3(aa,mu,sd,pp):
        zz=norm.ppf(pp, loc=0, scale=1)
        return aa+np.exp(mu+sd*zz)
    # parameter calculation
    aa,mu,sd=func_ln3(xx)
    # estimation
    ye=est_ln3(aa,mu,sd,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_ln3(aa,mu,sd,prob)
    res=np.array([aa,mu,sd,rr])
    ss='LN3\na={0:8.3f}\nm={1:8.3f}\ns={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def lp3(xx,pp,pyear):
    def func_lp3(xx):
        n=len(xx)
        yy=np.log(xx)
        my=np.sum(yy)/n
        sy=np.sqrt(np.sum((yy-my)**2)/n)
        cy=np.sum(((yy-my)/sy)**3)/n
        sy=np.sqrt(n/(n-1))*sy
        a=1+6.51/n+20.2/n**2
        b=1.48/n+6.77/n**2
        ry=cy*(a+b*cy**2)
        bb=4/ry**2
        aa=sy/np.sqrt(bb)
        if ry<0: aa=-aa
        cc=my-aa*bb
        return bb,cc,aa
    def est_lp3(bb,cc,aa,pp):
        ww=gamma.ppf(pp, bb, loc=0, scale=1)
        return np.exp(cc+aa*ww)
    # parameter calculation
    bb,cc,aa=func_lp3(xx)
    # estimation
    ye=est_lp3(bb,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_lp3(bb,cc,aa,prob)
    res=np.array([bb,cc,aa,rr])
    ss='LP3\nb={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gev(xx,pp,pyear):
    def func_gev(xx):
        n=len(xx)
        jj=np.arange(n)+1
        b0=np.sum(xx)/n
        b1=np.sum((jj-1)*xx)/n/(n-1)
        b2=np.sum((jj-1)*(jj-2)*xx)/n/(n-1)/(n-2)
        lam1=b0
        lam2=2*b1-b0
        lam3=6*b2-6*b1+b0
        d=2*lam2/(lam3+3*lam2)-np.log(2)/np.log(3)
        k=7.8590*d+2.9554*d**2
        a=k*lam2/(1-2**(-k))/ssp.gamma(1+k)
        c=lam1-a/k*(1-ssp.gamma(1+k))
        return k,c,a
    def est_gev(kk,cc,aa,pp):
        return cc+aa/kk*(1-(-np.log(pp))**kk)
    # parameter calculation
    kk,cc,aa=func_gev(xx)
    # estimation
    ye=est_gev(kk,cc,aa,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gev(kk,cc,aa,prob)
    res=np.array([kk,cc,aa,rr])
    ss='GEV\nk={0:8.3f}\nc={1:8.3f}\na={2:8.3f}\nr={3:8.3f}\nn={4:4.0f}'.format(res[0],res[1],res[2],res[3],len(xx))
    return res,ye,yp,ss


def gum(xx,pp,pyear):
    def func_gum(xx):
        n=len(xx)
        b0=np.sum(xx)/n
        j=np.arange(0,n)
        b1=np.sum(j*xx)/n/(n-1)
        lam1=b0
        lam2=2*b1-b0
        aa=lam2/np.log(2)
        cc=lam1-0.5772*aa
        return aa,cc
    def est_gum(aa,cc,pp):
        return cc-aa*np.log(-np.log(pp))
    aa,cc=func_gum(xx)
    # estimation
    ye=est_gum(aa,cc,pp)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_gum(aa,cc,prob)
    res=np.array([aa,cc,rr])
    ss='GUM\na={0:8.3f}\nc={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss


def sqe(xx,pp,pyear):
    def est_sqrtet(aa,bb,pp,x0):
        # estimation of quantile by Newton-Raphson Method
        def gx(x,aa,p):
            ggx=np.log(1+x)-x-np.log(-1/aa*np.log(p))
            gdx=1/(1+x)-1
            return ggx,gdx
        imax=100
        eps=1.0e-10
        tp=np.zeros(len(pp),dtype=np.float64)
        x1=x0
        for i,p in enumerate(pp):
            for j in range(imax):
                ggx,gdx=gx(x1,aa,p)
                x2=x1-ggx/gdx
                if np.abs(x2-x1)<eps: break
                x1=x2
            tp[i]=x2
        return tp**2/bb

    def bisection(xx):
        # estimation of coefficient bb by Bisection Method
        def func(bb,xx):
            # function for Bisection Method
            nn=len(xx)
            a1=(np.sum(np.sqrt(bb*xx))-2*nn)/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
            a2=nn/(np.sum(np.exp(-np.sqrt(bb*xx)))+np.sum(np.sqrt(bb*xx)*np.exp(-np.sqrt(bb*xx))))
            return a1-a2
        imax=100    # maximum number of iterations
        err=1.0e-10 # allowable error of the zero value
        x1=(2*len(xx)/np.sum(np.sqrt(xx)))**2
        x2=x1+0.5
        for i in range(0,imax):
            xi=0.5*(x1+x2)
            f1=func(x1,xx)
            f2=func(x2,xx)
            fi=func(xi,xx)
            if 0<f1*f2:
                x1=x2
                x2=x1+0.5
            elif abs(x2-x1) < err:
                break # converged
            else:
                if f1*fi < 0.0: x2=xi
                if fi*f2 < 0.0: x1=xi
        return xi

    bb=bisection(xx)
    aa=(np.sum(np.sqrt(bb*xx))-2*len(xx))/np.sum((bb*xx)*np.exp(-np.sqrt(bb*xx)))
    # estimation
    x0=np.sqrt(bb*xx[-1])
    ye=est_sqrtet(aa,bb,pp,x0)
    rr=np.corrcoef(xx, ye)[0][1]
    # probable value
    prob=1.0-1.0/pyear
    yp=est_sqrtet(aa,bb,prob,x0)
    res=np.array([aa,bb,rr])
    ss='SQRT-ET\na={0:8.3f}\nb={1:8.3f}\nr={2:8.3f}\nn={3:4.0f}'.format(res[0],res[1],res[2],len(xx))
    return res,ye,yp,ss

クオンタイルプロット

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology


fsz=12
xmin=0; xmax=5000; dx=1000
ymin=0; ymax=5000; dy=1000


def qqplot(fnameR,xx,ye,ss,n):
    tstr=fnameR
    nnn=230+n
    plt.subplot(nnn)
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xlabel('Observed (m$^3$/s)')
    plt.ylabel('Predicted (m$^3$/s)')
    plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(color='#999999',linestyle='solid')
    plt.gca().set_aspect('equal',adjustable='box')
    plt.title(tstr,loc='left',fontsize=fsz)
    plt.plot(xx,ye,'o',clip_on=False)
    plt.plot([xmin,xmax],[ymin,ymax],'-',color='#ff0000',lw=1)
    xs=xmin+0.05*(xmax-xmin)
    ys=ymin+0.95*(ymax-ymin)
    plt.text(xs,ys,ss,va='top',ha='left',fontsize=fsz,fontname='Ricty Diminished')


def calc(fnameR,alpha,pyear):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    res1,ye1,yp1,ss1=im_hydrology.ln3(xx,pp,pyear)    
    res2,ye2,yp2,ss2=im_hydrology.lp3(xx,pp,pyear)    
    res3,ye3,yp3,ss3=im_hydrology.gev(xx,pp,pyear)
    res4,ye4,yp4,ss4=im_hydrology.gum(xx,pp,pyear)
    res5,ye5,yp5,ss5=im_hydrology.sqe(xx,pp,pyear)
    
    plt.figure(figsize=(12,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    qqplot(fnameR,xx,ye1,ss1,1)
    qqplot(fnameR,xx,ye2,ss2,2)
    qqplot(fnameR,xx,ye3,ss3,3)
    qqplot(fnameR,xx,ye4,ss4,4)
    qqplot(fnameR,xx,ye5,ss5,5)
    plt.tight_layout()
    s=fnameR.replace('inp_','')
    s=s.replace('.txt','')
    fnameF='fig_qq_'+s+'.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()    

    s0='{0:8s}'.format('Year')
    s1='{0:8s}'.format('LN3')
    s2='{0:8s}'.format('LP3')
    s3='{0:8s}'.format('GEV')
    s4='{0:8s}'.format('GUM')
    s5='{0:8s}'.format('SQRT-ET')
    for i in range(len(pyear)):
        s0=s0+',{0:6.0f}'.format(pyear[i])
        s1=s1+',{0:6.0f}'.format(yp1[i])
        s2=s2+',{0:6.0f}'.format(yp2[i])
        s3=s3+',{0:6.0f}'.format(yp3[i])
        s4=s4+',{0:6.0f}'.format(yp4[i])
        s5=s5+',{0:6.0f}'.format(yp5[i])
    print(s0)
    print(s1)
    print(s2)
    print(s3)
    print(s4)
    print(s5)


def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        calc(fnameR,alpha,pyear)


#==============
# Execution
#==============
if __name__ == '__main__': main()

確率年と予測値

# Jackknife method
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import sys; sys.path.append('/Users/kk/pypro')
import im_hydrology

# Global variables
fsz=14
xmin=1
xmax=1000
dx=1
ymin=0
ymax=10000
dy=1000


def drawfig0(pyear,qqq,smt,ss):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qqq[0,:],'o--',ms=10,color=col,markerfacecolor=col,label=smt[0],lw=1,clip_on=False)
    plt.plot(pyear,qqq[1,:],'s--',ms=10,color=col,markerfacecolor=col,label=smt[1],lw=1,clip_on=False)
    plt.plot(pyear,qqq[2,:],'^--',ms=10,color=col,markerfacecolor=col,label=smt[2],lw=1,clip_on=False)
    plt.plot(pyear,qqq[3,:],'o-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[3],lw=1,clip_on=False)
    plt.plot(pyear,qqq[4,:],'s-',ms=10,color=col,markerfacecolor='#ffffff',label=smt[4],lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    fnameF='fig_comp.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()


def drawfig(pyear,qb,sd,ss,kkk):
    plt.figure(figsize=(10,8),facecolor='w')
    plt.rcParams['font.size']=fsz
    plt.rcParams['font.family']='sans-serif'
    tstr='Return Period vs Maximum Flood Peak Discharge at Tamor Intake'
    if kkk==1: tstr=tstr+' (LN3)'
    if kkk==2: tstr=tstr+' (LP3)'
    if kkk==3: tstr=tstr+' (GEV)'
    if kkk==4: tstr=tstr+' (GUM)'
    if kkk==5: tstr=tstr+' (SQRT-ET)'
    plt.xlim([xmin,xmax])
    plt.ylim([ymin,ymax])
    plt.xscale('log')
    plt.xlabel('Retuen Period (years)')
    plt.ylabel('Maximum Flood Peak Discharge (m$^3$/s)')
    #plt.xticks(np.arange(xmin,xmax+dx,dx))
    plt.xticks([1,10,100,1000],[1,10,100,1000])
    plt.yticks(np.arange(ymin,ymax+dy,dy))
    plt.grid(which='both',color='#999999',linestyle='solid')
    plt.title(tstr,loc='left',fontsize=fsz)

    col='#000080'
    plt.plot(pyear,qb   ,'o-'  ,ms=10,color=col,markerfacecolor=col,label='Average',lw=1,clip_on=False)
    plt.plot(pyear,qb+sd,'^--' ,ms=10,color=col,markerfacecolor=col,label='Ave.+ SD',lw=1,clip_on=False)
    plt.plot(pyear,qb-sd,'s--' ,ms=10,color=col,markerfacecolor=col,label='Ave.- SD',lw=1,clip_on=False)
    xs=1.2; ys=ymin+0.98*(ymax-ymin)
    bbox_props = dict(boxstyle='square,pad=0.1', fc='#ffffff', ec='#777777', lw=1)
    plt.text(xs,ys,ss,ha='left',va='top',fontsize=fsz,bbox=bbox_props,linespacing=1,fontfamily='monospace')
    plt.legend(loc='lower right',fontsize=fsz-2,handlelength=4)

    plt.tight_layout()
    if kkk==1: fnameF='fig_rp_ln3.jpg'
    if kkk==2: fnameF='fig_rp_lp3.jpg'
    if kkk==3: fnameF='fig_rp_gev.jpg'
    if kkk==4: fnameF='fig_rp_gum.jpg'
    if kkk==5: fnameF='fig_rp_sqe.jpg'
    plt.savefig(fnameF, dpi=100, bbox_inches="tight", pad_inches=0.1)
    plt.show()

    
def calc(fnameR,alpha,pyear,kkk):
    df = pd.read_csv(fnameR, index_col=False, delimiter='\s+')
    xx=np.array(df['Q2'])
    xx=np.sort(xx)
    # setting protting position
    n=len(xx)
    ii=np.arange(n)+1
    pp=(ii-alpha)/(n+1-2*alpha)
    if kkk==1: res,ye0,yp0,ss=im_hydrology.ln3(xx,pp,pyear)    
    if kkk==2: res,ye0,yp0,ss=im_hydrology.lp3(xx,pp,pyear)    
    if kkk==3: res,ye0,yp0,ss=im_hydrology.gev(xx,pp,pyear)    
    if kkk==4: res,ye0,yp0,ss=im_hydrology.gum(xx,pp,pyear)
    if kkk==5: res,ye0,yp0,ss=im_hydrology.sqe(xx,pp,pyear)
    print(ss)
    
    qq=np.zeros((n,len(pyear)),dtype=np.float64)
    qm=np.zeros(len(pyear),dtype=np.float64)
    qb=np.zeros(len(pyear),dtype=np.float64)
    sd=np.zeros(len(pyear),dtype=np.float64)
    for i in range(0,n):
        xi=np.delete(xx,i)
        m=len(xi)
        ii=np.arange(m)+1
        pp=(ii-alpha)/(m+1-2*alpha)
        if kkk==1: res,ye,yp,ss=im_hydrology.ln3(xi,pp,pyear)
        if kkk==2: res,ye,yp,ss=im_hydrology.lp3(xi,pp,pyear)
        if kkk==3: res,ye,yp,ss=im_hydrology.gev(xi,pp,pyear)
        if kkk==4: res,ye,yp,ss=im_hydrology.gum(xi,pp,pyear)
        if kkk==5: res,ye,yp,ss=im_hydrology.sqe(xi,pp,pyear)
        qq[i,:]=yp
    for j in range(0,len(pyear)):
        qm[j]=np.average(qq[:,j])
        qb[j]=n*yp0[j]-(n-1)*qm[j]
        sd[j]=np.sqrt((n-1)/n*np.sum((qq[:,j]-qm[j])**2))
    return qb,sd
    

def main():
    # coefficient for plotting position formula
    alpha=0.0 # Weibull
    #alpha=0.4 # Cunnane
    #alpha=0.5 # Hazen
    pyear=np.array([2,5,10,20,50,100,200,500,1000])
    qqq=np.zeros((5,len(pyear)),dtype=np.float64)
    lfile=['inp_flood.txt']
    for fnameR in lfile:
        for kkk in [1,2,3,4,5]:
            qb,sd=calc(fnameR,alpha,pyear,kkk)
            s0='{0:4s}'.format('Year')
            s1='{0:4s}'.format('Ave.')
            s2='{0:4s}'.format('SD')
            for j in range(len(pyear)):
                s0=s0+'{0:6.0f}'.format(pyear[j])
                s1=s1+'{0:6.0f}'.format(qb[j])
                s2=s2+'{0:6.0f}'.format(sd[j])
            ss='Retuen Period and Flood Discharge (m$^3$/s)'
            ss=ss+'\n'+s0+'\n'+s1+'\n'+s2
            ss=ss+'\n* Ave.: average, SD: standard deviation'
            print(ss)
            drawfig(pyear,qb,sd,ss,kkk)
            qqq[kkk-1,:]=qb[:]
        smt=['LN3','LP3','GEV','GUM','SQRT-ET']
        s0='{0:8s}'.format('Year')
        s1='{0:8s}'.format(smt[0])
        s2='{0:8s}'.format(smt[1])
        s3='{0:8s}'.format(smt[2])
        s4='{0:8s}'.format(smt[3])
        s5='{0:8s}'.format(smt[4])
        for j in range(len(pyear)):
            s0=s0+'{0:6.0f}'.format(pyear[j])
            s1=s1+'{0:6.0f}'.format(qqq[0,j])
            s2=s2+'{0:6.0f}'.format(qqq[1,j])
            s3=s3+'{0:6.0f}'.format(qqq[2,j])
            s4=s4+'{0:6.0f}'.format(qqq[3,j])
            s5=s5+'{0:6.0f}'.format(qqq[4,j])
        ss='Retuen Period and Flood Discharge (m$^3$/s)'
        ss=ss+'\n'+s0+'\n'+s1+'\n'+s2+'\n'+s3+'\n'+s4+'\n'+s5
        print(ss)
        drawfig0(pyear,qqq,smt,ss)


#==============
# Execution
#==============
if __name__ == '__main__': main()

Thank you.

記事の先頭に行く