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

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

<table>
<tr><th>Non-exceedance<br>probbility</th>
<th>Return period<br>(years)</th>
<th>Peak discharge<br>(m3/s)</th></tr>
<tr><td>0.50 </td><td>   2</td><td> 950</td></tr>
<tr><td>0.80 </td><td>   5</td><td>1650</td></tr>
<tr><td>0.90 </td><td>  10</td><td>2190</td></tr>
<tr><td>0.95 </td><td>  20</td><td>2780</td></tr>
<tr><td>0.98 </td><td>  50</td><td>3610</td></tr>
<tr><td>0.99 </td><td> 100</td><td>4330</td></tr>
<tr><td>0.995</td><td> 200</td><td>5090</td></tr>
<tr><td>0.998</td><td> 500</td><td>6190</td></tr>
<tr><td>0.999</td><td>1000</td><td>7090</td></tr>
</table>
Note: based on 2-parameter lognormal distribution.
<br>

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

+ 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)
```

#プログラム

```python: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)
```

![fig_flood_p.png](https://qiita-image-store.s3.amazonaws.com/0/129300/9ece9855-c9c2-34e0-e68c-81e2364e4cb5.png)



# 参考サイト
Scipyの統計関数
[http://kaisk.hatenadiary.com/entry/2015/02/17/192955](http://kaisk.hatenadiary.com/entry/2015/02/17/192955)

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

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