Python, matplotlibにより対数正規確率プロットを行う

はじめに

報告書に以下のようなデータが記載されている。

Non-exceedance probbility	Return period (years)	Peak discharge (m3/s)
0.50	2	950
0.80	5	1650
0.90	10	2190
0.95	20	2780
0.98	50	3610
0.99	100	4330
0.995	200	5090
0.998	500	6190
0.999	1000	7090

Note: based on 2-parameter lognormal distribution.

確率洪水流量は２変数対数正規分布で求めているようだ。
このデータを用いて、以下のことを行ってみる事にする。

• 10000年確率洪水流量を推定する
• ２変数対数正規分布に従っているか確認する

データに対して外挿することになるが、限られた情報の中で、勘で数字を出しましたというよりはマシであろう。まあ、プログラムの練習なので、統計的な厳密さには目をつぶることにして欲しい。

検討の流れ

回帰分析の準備

回帰分析の準備として，標準正規分布における非超過確率に対応するパーセント点，および流量の常用対数値を計算する．

ppp=norm.ppf(Pin, loc=0, scale=1)
qqq=np.log10(Qin)

回帰分析

次に目的変数を流量の対数値，説明変数を非超過確率のパーセント点として線形回帰を行ない、10000年確率に相当する洪水流量を推定する

scipy: optimize.leastsqを用いた線形回帰による回帰係数，相関係数の算出

parameter0 = [0.0,0.0]
result = optimize.leastsq(func,parameter0,args=(ppp,qqq))
aa=result[0][0]
bb=result[0][1]
rr=np.corrcoef(ppp,qqq)

optimize.leastsqのための関数（イコール０となる形式で表現する）

def func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    residual = y-(a*x+b)
    return residual

回帰結果を用いた1000年確率の推定

pp=0.9999
xq=norm.ppf(pp, loc=0, scale=1)
Qpp=10**(aa*xq+bb)

対数正規確率紙へのプロット

対数正規分布に従っているかどうかは、対数正規確率紙にプロットし、確認する。ここでは、グラフの横軸に洪水流量(対数値)、縦軸に確率値をプロットする。このため、デフォルトの軸目盛りは削除し、横軸、縦軸とも独自に設定する。

軸目盛りおよび目盛り線の削除

plt.tick_params(labelbottom='off')
plt.tick_params(labelleft='off')
plt.tick_params(which='both', width=0)

グラフはpng画像として保存するが、Qiitaなどに貼ることを考慮して、余白を削除しておく。

plt.savefig('fig_flood_p.png', bbox_inches="tight", pad_inches=0.2)

プログラム

py_qq.py

import numpy as np
from scipy.stats import norm # normal distribution
from scipy import optimize
import matplotlib.pyplot as plt


# function for optimize.leastsq
def func(parameter,x,y):
    a = parameter[0]
    b = parameter[1]
    residual = y-(a*x+b)
    return residual

#==============================
# data analysis
#==============================
# flood discharge data
Qin=np.array([950,1650,2190,2780,3610,4330,5090,6190,7090])
# non-exceedance probability data
Pin=np.array([0.50,0.80,0.90,0.95,0.98,0.99,0.995,0.998,0.999])

ppp=norm.ppf(Pin, loc=0, scale=1)
qqq=np.log10(Qin)

# least square method
parameter0 = [0.0,0.0]
result = optimize.leastsq(func,parameter0,args=(ppp,qqq))
aa=result[0][0]
bb=result[0][1]
rr=np.corrcoef(ppp,qqq)

# estimation of 10,000 years return period flood
pp=0.9999
xq=norm.ppf(pp, loc=0, scale=1)
Qpp=10**(aa*xq+bb)

print('log10(Q) = a * x + b')
print('a={aa:10.6f} b={bb:10.6f} r={rr[0][1]:10.6f}'.format(**locals()))
print('Q={Qpp:10.3f} (pp={pp:0.4f})'.format(**locals()))


#==============================
# plot by matplotlib
#==============================
fig = plt.figure(figsize=(7,8))
xmin=np.log10(100)
xmax=np.log10(20000)
ymin=norm.ppf(0.0001, loc=0, scale=1)
ymax=norm.ppf(0.9999, loc=0, scale=1)
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.tick_params(labelbottom='off')
plt.tick_params(labelleft='off')
plt.tick_params(which='both', width=0)

# data plot
plt.plot(qqq,ppp,'o')
# straight line by regression analysis
plt.plot([xmin, xmax], [(xmin-bb)/aa, (xmax-bb)/aa], color='k', linestyle='-', linewidth=1)

# setting of x-y axes
_dy=np.array([0.0001,0.001,0.01,0.1,0.5,0.9,0.99,0.999,0.9999])
dy=norm.ppf(_dy, loc=0, scale=1)
plt.hlines(dy, xmin, xmax, color='grey')
_dx=np.array([100,200,300,400,500,600,700,800,900,1000,2000,3000,4000,5000,6000,7000,8000,9000,10000,20000])
dx=np.log10(_dx)
plt.vlines(dx, ymin, ymax, color='grey')
fs=16
for i in range(2,5):
    plt.text(float(i), ymin-0.1, str(10**i), ha = 'center', va = 'top', fontsize=fs)
for i in range(0,9):
    plt.text(xmin-0.01, dy[i], str(_dy[i]), ha = 'right', va = 'center', fontsize=fs)
plt.text(0.5*(xmin+xmax), ymin-0.5, 'Flood discharge (m$^3$/s)', ha = 'center', va = 'center', fontsize=fs)
plt.text(xmin-0.25,0.5*(ymin+ymax),'Non-exceedance probability', ha = 'center', va = 'center', fontsize=fs, rotation=90)

# image saving and image showing
plt.savefig('fig_flood_p.png', bbox_inches="tight", pad_inches=0.2)
plt.show()

出力結果

$ python3 py_qq.py
log10(Q) = a * x + b
a=  0.282460 b=  2.978647 r=  0.999996
Q= 10693.521 (pp=0.9999)

参考サイト

Scipyの統計関数
http://kaisk.hatenadiary.com/entry/2015/02/17/192955

Scipyによる線形回帰事例
http://www2.kaiyodai.ac.jp/~kentaro/materials/new_HP/python/15fit_data3.html

matplotlibで画像の余白を削除する
https://mzmttks.blogspot.my/2012/01/pylab-2.html