#はじめに
報告書に以下のようなデータが記載されている。
| 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 |
確率洪水流量は２変数対数正規分布で求めているようだ。
このデータを用いて、以下のことを行ってみる事にする。
- 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