2018.10.25 プログラムを母数推定部（計算のみ）と作図部の２部に書き換え
はじめに
2018年10月10日に海外での業務完了につき帰国した．現在は休暇中であり，のんびりしているが，全く何もやらないのもちょっと気が引けたので，仕事の練習として，Pythonを使った非線形回帰のプログラムを作ってみた．昔は同じことをFortranプログラムでやっていたが，scipyの統計関数や，非線形回帰機能を使うことにより，非常にシンプルにやりたいことを実現できることを確認した．
やっていることは，観測値を，仮定した確率分布関数に回帰し，母数（確率分布関数における係数）を推定するというものである．母数推定の方法は，これまでに積率法（モーメント法）などが研究されてきたが，ここでは Python：scipy の optimize.leastsq を用いる方法をとっている．
例として用いる観測値は，筆者の自宅に近い前橋地方気象台の各年雨量データ（前橋の気象データ）とした．
例題の中で扱った確率分布は，以下のものである．年最大日雨量データと下記の確率分布の関係を取扱った．
３変数対数正規分布 Log-Normal distribution with 3 parameters LN3 (ln3)
２変数対数正規分布 Log-Normal distribution with 2 parameters LN2 (ln2)
一般化極値分布 Generalized Extreme Value distribution GEV (gev)
グンベル分布 Gumbel distribution Gumbel (gum)
一般化パレート分布 Generalized Pareto distribution GPD (gpd)
指数分布 Exponential distribution Exponential (exp)
対数ピアソンIII型分布 Log-Pearson type III distribution LP3 (lp3)
３変数ワイブル分布 Weibull distribution with 3 parameters Weibull (wei)
なお，年最大日雨量データは，正規分布との適合性は良好ではないことが一般に知られている．このため，正規分布との適合性が比較的良い，各年の年間降雨量と正規分布との関係も扱ってみた．
当方の計算環境は以下の通り．
- MacBook Pro (Retina, 13-inch, Mid 2014)
- macOS Mojave
- Python 3.7.0
参考文献
データの取得
データとして，国土交通省気象庁のホームページ（前橋の気象データ）から，前述の前橋地方気象台における1987年から2017年までの観測データを取得した．ここでは難しいことはせず，ホームページよりテキストをクリップボードにコピーし，テキストエディタと簡単なプログラムにより整形して用いている．
作成したデータファイルの書式は以下の通り．
inp_maebashi.txt
year rfy rfd tave tmax tmin
1897 1384.6 67.5 12.5 17.7 8.3
1898 1099.1 107.6 13.3 18.5 8.9
1899 1097.0 66.3 13.0 18.5 8.7
.... ... ... ... ... ...
year : 観測年
rfy : 年間降雨量
rfd : 年最大日雨量
tave : 年平均気温
tmax : 年最高気温
tmin : 年最低気温
まずは，前段として，降雨量と気温の経時変化をみてみよう．降雨量については特に近年多雨化するような傾向は見られないが，気温については1950年頃から明らかに上昇傾向にある．地球温暖化の影響か？ この件については時間があったらということで，この投稿では，日最大降雨量（変数名：rfd）および年降雨量（変数名：rfy）について，確率分布に回帰する問題を考える．
プログラム中では，降雨量・気温を含めたファイルからデータを読み込んでいるが，興味ある読者が遊べるよう，分析に用いているデータを，１次元numpy配列として以下に示しておく．配列は1987年から2017年の値を順番に並べたものであり，ソートはかけていない．
# 前橋の各年の最大日降雨量（mm/day）1987〜2017年
import numpy as np
xx=np.array([67.5, 107.6, 66.3, 115.2, 123.7, 96.5, 93.2, 135.8, 68.6, 73.4, 71.5, 72.3,
52.4, 95.5, 164.3, 65.5, 110.1, 140.1, 104.2, 71.2, 59.9, 113.9, 67.8, 113.3,
94.2, 91.7, 92.3, 61.0, 98.0, 63.4, 99.7, 74.0, 119.8, 75.4, 76.5, 71.5,
50.8, 74.3, 95.7, 90.6, 141.7, 69.3, 84.7, 130.7, 163.1, 52.1, 112.6, 148.0,
128.5, 46.9, 357.4, 115.1, 75.6, 142.1, 67.6, 43.2, 69.6, 146.5, 262.4, 78.1,
50.0, 180.0, 122.8, 48.2, 106.2, 147.2, 48.7, 92.4, 70.4, 145.7, 61.0, 75.0,
89.0, 68.0, 100.5, 75.5, 73.0, 82.5, 94.5, 55.0, 98.0, 59.0, 72.0, 44.5,
111.5, 175.0, 115.5, 57.0, 71.0, 102.5, 78.5, 63.5, 80.0, 95.0, 89.0, 90.5,
85.0, 107.0, 56.5, 78.5, 157.5, 134.0, 201.0, 82.5, 103.5, 139.5, 53.5, 104.0,
100.5, 76.5, 82.0, 66.5, 64.5, 64.0, 178.5, 83.0, 84.5, 104.0, 70.5, 105.5,
91.0])
print(xx)
# 前橋の年降雨量（mm/year）1987〜2017年
import numpy as np
xx=np.array([1384.6, 1099.1, 1097.0, 927.9, 1235.6, 1574.0, 1213.6, 1308.1, 1304.4, 1272.4,
1292.0, 1231.5, 1284.4, 1733.0, 1396.4, 1426.5, 1212.3, 1230.0, 1534.1, 1556.5,
1327.5, 1078.8, 1321.8, 1538.0, 1526.4, 1161.2, 1459.7, 1008.8, 1339.6, 968.3,
1233.5, 1389.2, 1298.2, 1184.8, 1138.8, 1352.5, 911.3, 968.9, 1242.7, 1376.0,
1169.5, 1425.9, 1128.4, 1133.1, 1492.3, 868.2, 1168.5, 949.2, 1356.4, 1084.2,
1363.3, 1545.8, 1172.4, 1497.0, 1139.1, 1164.1, 1274.6, 1275.3, 1776.6, 1310.5,
1234.4, 1465.1, 1529.8, 879.5, 1310.2, 950.4, 799.0, 1291.5, 1148.2, 1245.5,
1067.2, 1391.0, 1050.5, 984.5, 1030.0, 1051.5, 917.5, 1122.0, 1047.0, 986.0,
1186.0, 847.5, 1100.0, 1060.0, 1283.0, 1275.0, 1163.0, 816.5, 1140.5, 1240.5,
1191.5, 1356.5, 1657.5, 1201.0, 1438.5, 1174.0, 1346.5, 970.5, 1075.5, 815.5,
1078.5, 1649.5, 1495.0, 1163.0, 1316.0, 1503.0, 1104.5, 1196.0, 1114.0, 1479.0,
1310.5, 1425.0, 986.0, 1490.5, 1340.0, 1074.0, 998.5, 1395.5, 1232.0, 1249.0,
1192.5])
print(xx)
プロッティング・ポジション公式
統計量を確率分布に回帰する場合に必要となるプロッティング・ポジション（非超過確率に相当）を定める公式を示す．プロッティング・ポジション公式を用いるため，標本値 $x$ は，昇順に並べ替えておくことに注意する．
\begin{equation}
F[x_{(i)}]=\frac{i-\alpha}{N+1-2\alpha}
\end{equation}
\begin{align}
&N & & \text{標本数} \\
&i & & \text{標本値を昇順に並べたときの小さいほうからの順位} \\
&x_{(i)} & & \text{$i$番目の順位標本値} \\
&F[x_{(i)}] & & \text{プロッティング・ポジション(非超過確率相当)} \\
&\alpha & & \text{$0 \sim 1$の定数}
\end{align}
定数$\alpha$の値としては以下のものが知られている．本稿ではCunnane式（$=0.4$）を用いる．
\begin{align}
&\text{Weibull} & & \alpha=0 \\
&\text{Blom} & & \alpha=0.375 \\
&\text{Cunnane} & & \alpha=0.40 \\
&\text{Gringorten} & & \alpha=0.44 \\
&\text{Hazen} & & \alpha=0.50 \\
\end{align}
確率分布の数式表現
統計量$x$が各種確率分布に従う場合の，確率密度関数，確率分布函数，クオンタイルの数式表現を以下に示す．
３変数対数正規分布（Log-Normal distribution with 3 parameters: LN3）
確率密度関数
\begin{equation}
f(x)=\frac{1}{(x-a)\sigma_y\sqrt{2\pi}}\cdot\exp\left\{-\frac{1}{2}\left[\frac{\ln(x-a)-\mu_y}{\sigma_y}\right]^2\right\}
\end{equation}
確率分布関数
\begin{equation}
F(x)=\Phi\left(\frac{\ln(x-a)-\mu_y}{\sigma_y}\right)
\qquad \Phi(z)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^z\exp\left(-\frac{1}{2}t^2\right)dt
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
z=\frac{\ln(x-a)-\mu_y}{\sigma_y} \qquad\rightarrow\qquad x=a+\exp(\mu_y+\sigma_y z)
\end{equation}
\begin{equation}
x_p=a+\exp(\mu_y+\sigma_y z_p) \qquad \text{$z_p$は$p=\Phi(z)$となる$z$の値}
\end{equation}
２変数対数正規分布（Log-Normal distribution with 2 parameters: LN2）
３変数対数正規分布において，$a=0$ としたもの．
確率密度関数
\begin{equation}
f(x)=\frac{1}{x\sigma_y\sqrt{2\pi}}\cdot\exp\left\{-\frac{1}{2}\left[\frac{\ln(x)-\mu_y}{\sigma_y}\right]^2\right\}
\end{equation}
確率分布関数
\begin{equation}
F(x)=\Phi\left(\frac{\ln(x)-\mu_y}{\sigma_y}\right)
\qquad \Phi(z)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^z\exp\left(-\frac{1}{2}t^2\right)dt
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
z=\frac{\ln(x)-\mu_y}{\sigma_y} \qquad\rightarrow\qquad x=\exp(\mu_y+\sigma_y z)
\end{equation}
\begin{equation}
x_p=\exp(\mu_y+\sigma_y z_p) \qquad \text{$z_p$は$p=\Phi(z)$となる$z$の値}
\end{equation}
一般化極値分布（Generalized Extreme Value distribution: GEV）
確率密度関数
\begin{equation}
f(x)=\frac{1}{a}\left(1-k\frac{x-c}{a}\right)^{1/k-1}\cdot\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right]
\qquad (k\neq 0)
\end{equation}
$c$は位置母数，$a$は尺度母数，$k$は形状母数である．
確率分布関数
\begin{equation}
F(x)=\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right]
\qquad (k\neq 0)
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
p=\exp\left[-\left(1-k\frac{x-c}{a}\right)^{1/k}\right] \qquad\rightarrow\qquad x=c+\frac{a}{k}\cdot\left\{1-[-\ln(p)]^k\right\}
\end{equation}
\begin{equation}
x_p=c+\frac{a}{k}\cdot\left\{1-[-\ln(p)]^k\right\}
\end{equation}
グンベル分布（Gumbel distribution）
グンベル分布は，一般化極値分布において，$k=0$の場合と一致する．
確率密度関数
\begin{equation}
f(x)=\frac{1}{a}\exp\left[-\frac{x-c}{a}-\exp\left(-\frac{x-c}{a}\right)\right] \qquad -\infty<x<\infty
\end{equation}
$c$は位置母数，$a$は尺度母数である．
確率分布関数
\begin{equation}
F(x)=\exp\left[-\exp\left(-\frac{x-c}{a}\right)\right]
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
p=\exp\left[-\exp\left(-\frac{x-c}{a}\right)\right] \qquad\rightarrow\qquad x=c-a\ln[-\ln(p)]
\end{equation}
\begin{equation}
x_p=c-a\ln[-\ln(p)]
\end{equation}
一般化パレート分布（Generalized Pareto distribution: GPD）
確率密度関数
\begin{equation}
f(x)=\frac{1}{a}\left(1-k\frac{x-c}{a}\right)^{1/k-1} \qquad (k\neq 0)
\end{equation}
$c$は位置母数，$a$は尺度母数，$k$は形状母数である．
確率分布関数
\begin{equation}
F(x)=1-\left(1-k\frac{x-c}{a}\right)^{1/k}
\qquad (k\neq 0)
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
p=1-\left(1-k\frac{x-c}{a}\right)^{1/k}\qquad\rightarrow\qquad x=c+\frac{a}{k}\left\{1-(1-p)^k\right\}
\end{equation}
\begin{equation}
x_p=c+\frac{a}{k}\cdot\left\{1-(1-p)^k\right\}
\end{equation}
指数分布（Exponential distribution）
指数分布は，一般化パレート分布において，$k=0$の場合と一致する．
確率密度関数
\begin{equation}
f(x)=\frac{1}{a}\exp\left(-\frac{x-c}{a}\right)
\end{equation}
$c$は位置母数，$a$は尺度母数である．
確率分布関数
\begin{equation}
F(x)=1-\exp\left(-\frac{x-c}{a}\right)
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
p=1-\exp\left(-\frac{x-c}{a}\right)\qquad\rightarrow\qquad x=c-a\ln(1-p)
\end{equation}
\begin{equation}
x_p=c-a\ln(1-p)
\end{equation}
対数ピアソンIII型分布（Log-Pearson type III distribution: LP3）
確率密度関数
\begin{equation}
f(x)=\frac{1}{|a|\cdot\Gamma(b)\cdot x}\left(\frac{\ln x-c}{a}\right)^{b-1}\cdot\exp\left(-\frac{\ln x-c}{a}\right)
\qquad a>0 : \exp(c)<x<\infty
\end{equation}
$c$は位置母数，$a$は尺度母数，$b$は形状母数である．
確率分布関数
\begin{equation}
F(x)=G\left(\frac{\ln x-c}{a}\right)
\qquad G(w)=\frac{1}{\Gamma(b)}\int_0^w t^{b-1}\exp(-t)dt \qquad (a>0)
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
w=\frac{\ln x-c}{a} \qquad\rightarrow\qquad x=\exp(c+a w)
\end{equation}
\begin{equation}
x_p=\exp(c+a w_p) \qquad \text{$w_p$は$p=G(w)$となる$w$の値}
\end{equation}
３変数ワイブル分布（Weibull distribution with 3 parameters）
確率密度関数
\begin{equation}
f(x)=\frac{k}{a}\left(\frac{x-c}{a}\right)^{k-1}\exp{\left[-\left(\frac{x-c}{a}\right)^k\right]}
\qquad (k\neq 0)
\end{equation}
$c$は位置母数，$a$は尺度母数，$k$は形状母数である．
確率分布関数
\begin{equation}
F(x)=1-\exp\left[-\left(\frac{x-c}{a}\right)^k\right]
\qquad (k\neq 0)
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
1-p=\exp\left[-\left(\frac{x-c}{a}\right)^k\right] \qquad\rightarrow\qquad x=c+a[-\ln{(1-p)}]^{1/k}
\end{equation}
\begin{equation}
x_p=c+a[-\ln{(1-p)}]^{1/k}
\end{equation}
各種確率分布と最大日雨量の関係
以下にプログラム事例を示す．プログラムは以下の２部となっている．
- 母数推定（計算のみ）
- 作図（Q-Qプロット，確率分布曲線，確率密度曲線）
結果グラフ
Q-Qプロット
確率分布曲線
確率密度曲線
母数推定プログラム
母数推定は，前述「非超過確率$p$に対する水文統計量$x_p$」に示した式を用い，左辺$x_p$に標本値（観測値）$x$を，右辺$p$に標本値$x$に対応する非超過確率を入れ，左辺と右辺の差が最小となるよう，右辺の未定係数（母数）を optimize.leastsq で求めるよう行なっている．
母数推定計算部の簡単な説明（対数ピアソンIII型：LP3 の例）
主要な変数・関数の簡単な説明を，対数ピアソンIIIを例に，以下に記す．
- xx : 標本値（各年の日最大雨量観測値）
- pp : プロッティング・ポジション公式で求めた非超過確率
- parameter0 : 未定係数の初期値（３変数回帰なので３個）
- func_lp3 : 最小化する関数．この場合ガンマ分布のパラメータ bb も未知数であることに注意．
- parameter : 求めるべき未定係数
- r : Q-Qプロットにおける相関係数
- 計算結果は，リスト res にて戻す．
対数ピアソンIII型分布の母数推定に関するコードは以下の通り．optimize.leastsq の利用により，自作で回帰するのと比べ，コードは非常にシンプルである．
def lp3(xx,pp):
def func_lp3(parameter,xx,pp):
bb=parameter[0]
cc=parameter[1]
aa=parameter[2]
ww=gamma.ppf(pp, bb, loc=0, scale=1)
est=np.exp(cc+aa*ww)
residual=xx-est
return residual
# parameter calculation
npara=3; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_lp3,parameter0,args=(xx,pp))
bb=result[0][0]
cc=result[0][1]
aa=result[0][2]
ww=gamma.ppf(pp, bb, loc=0, scale=1)
yy=np.exp(cc+aa*ww)
rr=np.corrcoef(xx, yy)[0][1]
res=['lp3(b,c,a,r)',bb,cc,aa,rr]
return res
ワイブル分布（関数 wei）での初期値
ワイブル分布におけるleastsqに与える初期値は，後述する別解法（Step-1）に示す方法のうち，$c=0$ とおいたときの，$k$ と $a$ の値を与えている．
母数推定プログラム全文
py_calc.py
#======================
# Parameter Calculation
#======================
import sys
import numpy as np
from scipy.stats import norm # normal distribution
from scipy.stats import gamma # Gamma distribution
from scipy import optimize # for use of leastsq
def wdata(fnameW,res,f):
print('{0:>15s} {1:15.7e} {2:15.7e} {3:15.7e} {4:15.7e}'.format(res[0],res[1],res[2],res[3],res[4]),file=f)
print('{0:>15s} {1:10.3f} {2:10.3f} {3:10.3f} {4:10.3f}'.format(res[0],res[1],res[2],res[3],res[4]))
def ln3(xx,pp):
def func_ln3(parameter,xx,zz):
aa=parameter[0]
mu=parameter[1]
sd=parameter[2]
est=aa+np.exp(mu+sd*zz)
residual=xx-est
return residual
# parameter calculation
zz=norm.ppf(pp, loc=0, scale=1)
npara=3; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_ln3,parameter0,args=(xx,zz))
aa=result[0][0]
mu=result[0][1]
sd=result[0][2]
yy=aa+np.exp(mu+sd*zz)
rr=np.corrcoef(xx, yy)[0][1]
res=['ln3(a,m,s,r)',aa,mu,sd,rr]
return res
def gev(xx,pp):
def func_gev(parameter,xx,pp):
kk=parameter[0]
cc=parameter[1]
aa=parameter[2]
est=cc+aa/kk*(1-(-np.log(pp))**kk)
residual=xx-est
return residual
# parameter calculation
npara=3; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_gev,parameter0,args=(xx,pp))
kk=result[0][0]
cc=result[0][1]
aa=result[0][2]
yy=cc+aa/kk*(1-(-np.log(pp))**kk)
rr=np.corrcoef(xx, yy)[0][1]
res=['gev(k,c,a,r)',kk,cc,aa,rr]
return res
def gpd(xx,pp):
def func_gpd(parameter,xx,pp):
kk=parameter[0]
cc=parameter[1]
aa=parameter[2]
est=cc+aa/kk*(1-(1-pp)**kk)
residual=xx-est
return residual
# parameter calculation
npara=3; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_gpd,parameter0,args=(xx,pp))
kk=result[0][0]
cc=result[0][1]
aa=result[0][2]
yy=cc+aa/kk*(1-(1-pp)**kk)
rr=np.corrcoef(xx, yy)[0][1]
res=['gpd(k,c,a,r)',kk,cc,aa,rr]
return res
def lp3(xx,pp):
def func_lp3(parameter,xx,pp):
bb=parameter[0]
cc=parameter[1]
aa=parameter[2]
ww=gamma.ppf(pp, bb, loc=0, scale=1)
est=np.exp(cc+aa*ww)
residual=xx-est
return residual
# parameter calculation
npara=3; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_lp3,parameter0,args=(xx,pp))
bb=result[0][0]
cc=result[0][1]
aa=result[0][2]
ww=gamma.ppf(pp, bb, loc=0, scale=1)
yy=np.exp(cc+aa*ww)
rr=np.corrcoef(xx, yy)[0][1]
res=['lp3(b,c,a,r)',bb,cc,aa,rr]
return res
def ln2(xx,pp):
def func_ln2(parameter,xx,zz):
mu=parameter[0]
sd=parameter[1]
est=np.exp(mu+sd*zz)
residual=xx-est
return residual
# parameter calculation
zz=norm.ppf(pp, loc=0, scale=1)
npara=2; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_ln2,parameter0,args=(xx,zz))
mu=result[0][0]
sd=result[0][1]
yy=np.exp(mu+sd*zz)
rr=np.corrcoef(xx, yy)[0][1]
res=['ln2(-,m,s,r)',0,mu,sd,rr]
return res
def gum(xx,pp):
def func_gum(parameter,xx,pp):
cc=parameter[0]
aa=parameter[1]
est=cc-aa*np.log(-np.log(pp))
residual=xx-est
return residual
# parameter calculation
npara=2; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_gum,parameter0,args=(xx,pp))
cc=result[0][0]
aa=result[0][1]
yy=cc-aa*np.log(-np.log(pp))
rr=np.corrcoef(xx, yy)[0][1]
res=['gum(-,c,a,r)',0,cc,aa,rr]
return res
def exp(xx,pp):
def func_exp(parameter,xx,pp):
cc=parameter[0]
aa=parameter[1]
est=cc-aa*np.log(1-pp)
residual=xx-est
return residual
# parameter calculation
npara=2; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_exp,parameter0,args=(xx,pp))
cc=result[0][0]
aa=result[0][1]
yy=cc-aa*np.log(1-pp)
rr=np.corrcoef(xx, yy)[0][1]
res=['exp(-,c,a,r)',0,cc,aa,rr]
return res
def wei(xx,pp):
def func_wei(parameter,xx,pp):
kk=parameter[0]
cc=parameter[1]
aa=parameter[2]
est=cc+aa*(-np.log(1-pp))**(1/kk)
residual=xx-est
return residual
# parameter calculation
n=len(xx)
cc=0.0
y=np.log(-np.log(1-pp))
x=np.log(xx-cc)
c1=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/(n*np.sum(x**2)-(np.sum(x))**2)
c2=(np.sum(x**2)*np.sum(y)-np.sum(x)*np.sum(x*y))/(n*np.sum(x**2)-(np.sum(x))**2)
kk=c1
aa=np.exp(-c2/c1)
parameter0 = np.array([kk,cc,aa])
result = optimize.leastsq(func_wei,parameter0,args=(xx,pp))
kk=result[0][0]
cc=result[0][1]
aa=result[0][2]
yy=cc+aa*(-np.log(1-pp))**(1/kk)
rr=np.corrcoef(xx, yy)[0][1]
res=['wei(k,c,a,r)',kk,cc,aa,rr]
return res
def main():
args = sys.argv
fnameL=args[1] # Location name
fnameR='inp_'+fnameL+'.txt' # input data file
fnameW='res_'+fnameL+'.txt' # output data file
# coefficient for plotting position formula
#alpha=0.0 # Weibull
alpha=0.4 # Cunnane
#alpha=0.5 # Hazen
# data input
data=np.loadtxt(fnameR, skiprows=1)
xx =data[:,2]
xx=np.sort(xx)
# setting protting position
n=len(xx)
ii=np.arange(n)+1
pp=(ii-alpha)/(n+1-2*alpha)
#parameter calculation
print('{0:>15s} {1:>10s} {2:>10s} {3:>10s} {4:>10s}'.format('Dist.','k or b','c or m','a or s','r'))
f=open(fnameW,'w')
res=ln3(xx,pp); wdata(fnameW,res,f)
res=gev(xx,pp); wdata(fnameW,res,f)
res=gpd(xx,pp); wdata(fnameW,res,f)
res=lp3(xx,pp); wdata(fnameW,res,f)
res=ln2(xx,pp); wdata(fnameW,res,f)
res=gum(xx,pp); wdata(fnameW,res,f)
res=exp(xx,pp); wdata(fnameW,res,f)
res=wei(xx,pp); wdata(fnameW,res,f)
f.close()
#==============
# Execution
#==============
if __name__ == '__main__': main()
出力ファイル
res_maebashi.txt
ln3(a,m,s,r) 4.6206741e+01 3.6279534e+00 7.7413782e-01 9.8997190e-01
gev(k,c,a,r) -2.6407481e-01 7.6352934e+01 2.2461378e+01 9.9360710e-01
gpd(k,c,a,r) -1.7110038e-01 5.9531466e+01 3.1083527e+01 9.8645487e-01
lp3(b,c,a,r) 3.7085592e+00 3.7885439e+00 1.9107401e-01 9.9187839e-01
ln2(-,m,s,r) 0.0000000e+00 4.4621660e+00 4.4827030e-01 9.7247983e-01
gum(-,c,a,r) 0.0000000e+00 7.7602703e+01 3.3029941e+01 9.6098207e-01
exp(-,c,a,r) 0.0000000e+00 5.3290528e+01 4.3488789e+01 9.7840226e-01
wei(k,c,a,r) 8.4712424e-01 5.9621867e+01 3.4147990e+01 9.8109111e-01
作図プログラム
- 基データファイルから観測値を読み込む
- 母数推定プログラムで計算した結果を書き込んだファイルから母数を推定結果を読み込む
- 観測値・母数推定結果からプロット用データを作成しファイル出力
- ファイルよりデータを読み込み作図．
- Q-Qプロット描画：観測値（observed）と推定結果（predicted）の関係をプロット．
- 確率分布関数描画：観測値ー非超過確率関係を描画し，その上にQ-Qプロット作図プログラムで出力した確率分布曲線を重ね書きする．
- 確率分布描画：観測値のヒストグラムを描画し，その上にQ-Qプロット作図プログラムで出力した確率密度曲線を重ね書きする．ヒストグラム作成時は，正規化しておく．
母数推定プログラムの出力ファイル（_res.txt）の読み込みは，pands により以下のように行う．空白区切テキストなので sep='\s+' とし，Header はないので header=None とする．
fnameR='res_maebashi.txt'
df=pd.read_csv(fnameR, sep='\s+', header=None)
print(df)
pandas で読み込んだデータ（print(df)の結果）は以下の通り．自動で行名と列名が追加される．
0 1 2 3 4
0 ln3(a,m,s,r) 46.206741 3.627953 0.774138 0.989972
1 gev(k,c,a,r) -0.264075 76.352934 22.461378 0.993607
2 gpd(k,c,a,r) -0.171100 59.531466 31.083527 0.986455
3 lp3(b,c,a,r) 3.708559 3.788544 0.191074 0.991878
4 ln2(-,m,s,r) 0.000000 4.462166 0.448270 0.972480
5 gum(-,c,a,r) 0.000000 77.602703 33.029941 0.960982
6 exp(-,c,a,r) 0.000000 53.290528 43.488789 0.978402
7 wei(k,c,a,r) 0.847124 59.621867 34.147990 0.981091
作図プログラム全文
py_plot.py
#======================
# Plot
#======================
import sys
import numpy as np
from scipy.stats import norm # normal distribution
from scipy.stats import gamma # Gamma distribution
import scipy.special as sp # for use of Gamma function
import pandas as pd
import matplotlib.pyplot as plt
# Declaration of global variables
xmin=0; xmax=600; dx=100
fsz=12 # font size
def drawqq(xx,yy,tstr,ss,j):
# Q-Q plot
ymin=xmin; ymax=xmax; dy=dx
nplot=240+j
plt.subplot(nplot)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Observed (mm/day)')
plt.ylabel('Predicted (mm/day)')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
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.scatter(xx,yy,s=20,c='#000080',marker='o')
plt.plot([xmin,xmax],[ymin,ymax],'-',color='#999999',lw=1)
xs=xmin+(xmax-xmin)*0.03
ys=ymin+(ymax-ymin)*0.97
plt.text(xs,ys,ss,ha='left',va='top',fontname='monospace')
def drawcum(xx,pp,yy,ff,tstr,j):
# CDF plot
ymin=0; ymax=1.0; dy=0.25
nplot=240+j
plt.subplot(nplot)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Non-exceedance Probability')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
plt.scatter(xx,pp,s=10,c='#000080',marker='o',linewidth=0.5,clip_on=False,label='Observed')
plt.plot(yy,ff,'-',color='#ff0000',lw=1,label='Predicted')
plt.legend(loc='lower right')
def drawhist(xx,xpdf,ypdf,tstr,j):
# PDF plot + Histogram
ymin=0; ymax=0.02; dy=0.005
nplot=240+j
plt.subplot(nplot)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Normalized frequency')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
edges = range(xmin,xmax,10)
plt.hist(xx,color='#00ccff',bins=edges,density=True,label='Observed')
plt.plot(xpdf,ypdf,'-',color='#ff0000',lw=2,label='Predicted')
plt.legend(loc='upper right')
def wdata(fnameW,vx,vy):
f=open(fnameW,'w')
for i in range(len(vx)):
print('{0:15.7e} {1:15.7e}'.format(vx[i],vy[i]),file=f)
f.close()
def data_ln3(xx,pp,aa,mu,sd,rr):
# plot data
zz=norm.ppf(pp, loc=0, scale=1)
yy=aa+np.exp(mu+sd*zz)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (LN3)'
ss='a={0:7.3f}\n'.format(aa)
ss=ss+'m={0:7.3f}\n'.format(mu)
ss=ss+'s={0:7.3f}\n'.format(sd)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
zz=norm.ppf(ff, loc=0, scale=1)
yc=aa+np.exp(mu+sd*zz)
xpdf=np.arange(aa+1,xmax,1)
ypdf=1/(xpdf-aa)/sd/np.sqrt(2*np.pi)*np.exp(-1/2*((np.log(xpdf-aa)-mu)/sd)**2)
fnameW='_q_ln3.txt'; wdata(fnameW,xx,yy)
fnameW='_c_ln3.txt'; wdata(fnameW,yc,ff)
fnameW='_d_ln3.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_gev(xx,pp,kk,cc,aa,rr):
# plot data
yy=cc+aa/kk*(1-(-np.log(pp))**kk)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (GEV)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa/kk*(1-(-np.log(ff))**kk)
xpdf=np.arange(xmin+1,xmax,1)
coef=1-kk*(xpdf-cc)/aa
ypdf=1/aa*coef**(1/kk-1)*np.exp(-coef**(1/kk))
fnameW='_q_gev.txt'; wdata(fnameW,xx,yy)
fnameW='_c_gev.txt'; wdata(fnameW,yc,ff)
fnameW='_d_gev.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_gpd(xx,pp,kk,cc,aa,rr):
# plot data
yy=cc+aa/kk*(1-(1-pp)**kk)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (GPD)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa/kk*(1-(1-ff)**kk)
xpdf=np.arange(xmin+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=1/aa*(1-kk*coef)**(1/kk-1)
fnameW='_q_gpd.txt'; wdata(fnameW,xx,yy)
fnameW='_c_gpd.txt'; wdata(fnameW,yc,ff)
fnameW='_d_gpd.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_lp3(xx,pp,bb,cc,aa,rr):
# plot data
ww=gamma.ppf(pp, bb, loc=0, scale=1)
yy=np.exp(cc+aa*ww)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (LP3)'
ss='b={0:7.3f}\n'.format(bb)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
ww=gamma.ppf(ff, bb, loc=0, scale=1)
yc=np.exp(cc+aa*ww)
xpdf=np.arange(np.exp(cc),xmax,1)
coef=(np.log(xpdf)-cc)/aa
ypdf=1/np.abs(aa)/sp.gamma(bb)/xpdf*coef**(bb-1)*np.exp(-coef)
fnameW='_q_lp3.txt'; wdata(fnameW,xx,yy)
fnameW='_c_lp3.txt'; wdata(fnameW,yc,ff)
fnameW='_d_lp3.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_ln2(xx,pp,mu,sd,rr):
# plot data
zz=norm.ppf(pp, loc=0, scale=1)
yy=np.exp(mu+sd*zz)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (LN2)'
ss='m={0:7.3f}\n'.format(mu)
ss=ss+'s={0:7.3f}\n'.format(sd)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
zz=norm.ppf(ff, loc=0, scale=1)
yc=np.exp(mu+sd*zz)
xpdf=np.arange(mu+1,xmax,1)
ypdf=1/xpdf/sd/np.sqrt(2*np.pi)*np.exp(-1/2*((np.log(xpdf)-mu)/sd)**2)
fnameW='_q_ln2.txt'; wdata(fnameW,xx,yy)
fnameW='_c_ln2.txt'; wdata(fnameW,yc,ff)
fnameW='_d_ln2.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_gum(xx,pp,cc,aa,rr):
# plot data
yy=cc-aa*np.log(-np.log(pp))
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (Gumbel)'
ss='c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
yc=cc-aa*np.log(-np.log(ff))
xpdf=np.arange(xmin+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=1/aa*np.exp(-coef-np.exp(-coef))
fnameW='_q_gum.txt'; wdata(fnameW,xx,yy)
fnameW='_c_gum.txt'; wdata(fnameW,yc,ff)
fnameW='_d_gum.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_exp(xx,pp,cc,aa,rr):
# plot data
yy=cc-aa*np.log(1-pp)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (Exponential)'
ss='c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
yc=cc-aa*np.log(1-ff)
xpdf=np.arange(xmin+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=1/aa*np.exp(-coef)
fnameW='_q_exp.txt'; wdata(fnameW,xx,yy)
fnameW='_c_exp.txt'; wdata(fnameW,yc,ff)
fnameW='_d_exp.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def data_wei(xx,pp,kk,cc,aa,rr):
# plot data
yy=cc+aa*(-np.log(1-pp))**(1/kk)
rr=np.corrcoef(xx, yy)[0][1]
tstr='Q-Q plot (Weibull)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr)
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa*(-np.log(1-ff))**(1/kk)
xpdf=np.arange(cc+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=kk/aa*coef**(kk-1)*np.exp(-coef**kk)
fnameW='_q_wei.txt'; wdata(fnameW,xx,yy)
fnameW='_c_wei.txt'; wdata(fnameW,yc,ff)
fnameW='_d_wei.txt'; wdata(fnameW,xpdf,ypdf)
return ss
def main():
args = sys.argv
fnameL=args[1] # location
# coefficient for plotting position formula
#alpha=0.0 # Weibull
alpha=0.4 # Cunnane
#alpha=0.5 # Hazen
# data input
fnameR='inp_'+fnameL+'.txt'
data=np.loadtxt(fnameR, skiprows=1)
xx =data[:,2]
xx=np.sort(xx)
# setting protting position
n=len(xx)
ii=np.arange(n)+1
pp=(ii-alpha)/(n+1-2*alpha)
# parameter input
fnameR='res_'+fnameL+'.txt'
df=pd.read_csv(fnameR, sep='\s+', header=None)
print(df)
# maing plot data
ss=[]
aa=df.iloc[0,1]
mu=df.iloc[0,2]
sd=df.iloc[0,3]
rr=df.iloc[0,4]
s=data_ln3(xx,pp,aa,mu,sd,rr)
ss=ss+[s]
kk=df.iloc[1,1]
cc=df.iloc[1,2]
aa=df.iloc[1,3]
rr=df.iloc[1,4]
s=data_gev(xx,pp,kk,cc,aa,rr)
ss=ss+[s]
kk=df.iloc[2,1]
cc=df.iloc[2,2]
aa=df.iloc[2,3]
rr=df.iloc[2,4]
s=data_gpd(xx,pp,kk,cc,aa,rr)
ss=ss+[s]
bb=df.iloc[3,1]
cc=df.iloc[3,2]
aa=df.iloc[3,3]
rr=df.iloc[3,4]
s=data_lp3(xx,pp,bb,cc,aa,rr)
ss=ss+[s]
dm=df.iloc[4,1] # dummy
mu=df.iloc[4,2]
sd=df.iloc[4,3]
rr=df.iloc[4,4]
s=data_ln2(xx,pp,mu,sd,rr)
ss=ss+[s]
dm=df.iloc[5,1] # dummy
cc=df.iloc[5,2]
aa=df.iloc[5,3]
rr=df.iloc[5,4]
s=data_gum(xx,pp,cc,aa,rr)
ss=ss+[s]
dm=df.iloc[6,1] # dummy
cc=df.iloc[6,2]
aa=df.iloc[6,3]
rr=df.iloc[6,4]
s=data_exp(xx,pp,cc,aa,rr)
ss=ss+[s]
kk=df.iloc[7,1]
cc=df.iloc[7,2]
aa=df.iloc[7,3]
rr=df.iloc[7,4]
s=data_wei(xx,pp,kk,cc,aa,rr)
ss=ss+[s]
slist=['LN3','GEV','GPD','LN3','LN2','Gumbel','Exponential','Weibull']
flist=['ln3','gev','gpd','lp3','ln2','gum','exp','wei']
# Q-Q plot
fig=plt.figure(figsize=(16,8),facecolor='w')
for j,fn in enumerate(flist):
fnameR='_q_'+fn+'.txt'
data=np.loadtxt(fnameR, skiprows=0)
xp=data[:,0]
yp=data[:,1]
tstr=slist[j]
drawqq(xp,yp,tstr,ss[j],j+1)
plt.tight_layout()
fnameF='fig_'+fnameL+'_qq.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
# CDF plot
fig=plt.figure(figsize=(16,8),facecolor='w')
for j,fn in enumerate(flist):
fnameR='_c_'+fn+'.txt'
data=np.loadtxt(fnameR, skiprows=0)
xp=data[:,0]
yp=data[:,1]
tstr=slist[j]
drawcum(xx,pp,xp,yp,tstr,j+1)
plt.tight_layout()
fnameF='fig_'+fnameL+'_cdf.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
# PDFplot + Histogram
fig=plt.figure(figsize=(16,8),facecolor='w')
for j,fn in enumerate(flist):
fnameR='_d_'+fn+'.txt'
data=np.loadtxt(fnameR, skiprows=0)
xp=data[:,0]
yp=data[:,1]
tstr=slist[j]
drawhist(xx,xp,yp,tstr,j+1)
plt.tight_layout()
fnameF='fig_'+fnameL+'_hist.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
#==============
# Execution
#==============
if __name__ == '__main__': main()
３変数ワイブル分布における別解法
３変数ワイブル分布で leastsq を用いる場合に，推定値が初期値に対して敏感であり，初期値の値によっては意図しない結果が得られてしまう．例えば，前橋のデータの場合，parameter0=np,array([1,1,1]) では非現実的な回帰推定値が計算されてしまうが，parameter0=np.array([10,10,10]) の場合は，現実的な回帰推定値が得られる．
このため，別解法による計算を試みた．
この問題では推定すべき母数が3個あるため，まず最初に位置母数 $c$ を推定し既知とした上で，形状母数 $k$ と尺度母数 $a$ を推定する方法をとる．
（Step-1）母数 c の推定
標本値 $x$ の非超過確率 $F(x)$ は，選定したプロッティング・ポジション公式により算定する．
確率分布関数を変形して，
\begin{equation}
1-F(x)=\exp\left[-\left(\frac{x-c}{a}\right)^k\right]
\end{equation}
上式の両辺の対数を2回とることにより
\begin{equation}
\ln\{-\ln[1-F(x)]\}=k \ln(x-c)-k \ln a \qquad\rightarrow\qquad Y=A\cdot X+B
\end{equation}
上記より，標本数を$N$として，
\begin{equation}
Y_i=\ln\{-\ln[1-F(x_i)]\} \qquad X_i=\ln(x_i-c)
\end{equation}
\begin{equation}
A=\cfrac{N \sum(X_i Y_i)-\sum X_i\cdot \sum Y_i}{N \sum X_i{}^2-(\sum X_i)^2} \qquad
B=\cfrac{\sum X_i{}^2\cdot \sum Y_i-\sum X_i\cdot \sum(X_i Y_i)}{N \sum X_i{}^2-(\sum X_i)^2}
\end{equation}
\begin{equation}
r=\cfrac{N \sum(X_i Y_i)-\sum X_i \cdot \sum Y_i}{\sqrt{[N \sum X_i{}^2-(\sum X_i)^2]\cdot[N \sum Y_i{}^2-(\sum Y_i)^2]}}
\end{equation}
\begin{equation}
k=-A \qquad a=\exp(-B/A)
\end{equation}
標本値を $Y_i=\ln{-\ln[1-F(x_i)]}$，$X_i=\ln(x_i-c)$ として直線回帰することにより，$k=A$，$a=\exp(-B/A)$ として求めることができる．
この際，$c$ の値を，１から徐々に増加させ，直線回帰による相関係数 $r$ が最大となる時の $c$ を求める $c$ とする．対数関数の引数は正の実数であるため，$c$を変化させる範囲は，標本値$x$の最小値未満である必要がある．
この段階で，一応3母数 $k$，$a$，$c$ は求まったことになるが，更に推定精度を上げるため，$c$ のみを固定し，次の段階で再度 $k$，$a$ を算定する．
（Step-2）k および a の推定
$c$ を既知として，次式を用いて，optimize.leastsq により $k$ および $a$ を求める．非線形回帰の初期値としては，Step-1で求めた $[k, a]$ を用いることができる．
\begin{equation}
x_p=c+a[-\ln{(1-p)}]^{1/k}
\end{equation}
結果グラフ
Step-1での結果
この方法では，確率密度関数と観測値ヒストグラムの関係は比較的良好であるが，大きな観測値の予測精度が他の方法に比べて良くない．
Step-1 + Step-2の結果
この方法では，Step-1までの結果に比べ，確率密度関数と観測値ヒストグラムの関係は悪くなるが，大きな観測値の予測精度がStep-1までの結果に比べ，よくなっている．
Step-1 + Step-2 で推定した母数を初期値として leastsq を利用
この方法では，確率密度関数と観測値ヒストグラムの関係は悪くなるが，大きな観測値の予測精度は，今回比較している３手法のうち最も良い．また，前述のプログラム py_calc.py での計算結果と同じ値に収束している．
プログラム
- 関数 wei1：Step-1まで
- 関数 wei2：Step-1 + Step-2
- 関数 wei3：関数 wei2 で推定した母数を初期値として leastsq を利用
py_wei.py
#=========================
# Weibull distribution
#=========================
import sys
import numpy as np
from scipy.stats import norm # normal distribution
from scipy import optimize
import matplotlib.pyplot as plt
# Declaration of global variables
fnameL='temp'
xmin=0; xmax=600; dx=100
# coefficient for plotting position formula
#alpha=0.0 # Weibull
alpha=0.4 # Cunnane
#alpha=0.5 # Hazen
def drawqq(xx,yy,tstr,ss):
ymin=xmin; ymax=xmax; dy=dx
fsz=12 # font size
plt.subplot(131)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Observed (mm/day)')
plt.ylabel('Predicted (mm/day)')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
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,yy,'o')
plt.plot([xmin,xmax],[ymin,ymax],'-',color='#999999',lw=1)
xs=xmin+(xmax-xmin)*0.03
ys=ymin+(ymax-ymin)*0.97
plt.text(xs,ys,ss,ha='left',va='top',fontname='monospace')
def drawcum(xx,pp,yy,ff,tstr):
ymin=0; ymax=1.0; dy=0.25
fsz=12 # font size
plt.subplot(132)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Non-exceedance Probability')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
plt.plot(xx,pp,'x',color='#0000ff',markersize=3,label='Observed')
plt.plot(yy,ff,'-',color='#ff0000',lw=1,label='Predicted')
plt.legend()
def drawhist(xx,xpdf,ypdf,tstr):
ymin=0; ymax=0.02; dy=0.005
fsz=12 # font size
plt.subplot(133)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/day)')
plt.ylabel('Normalized frequency')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
edges = range(xmin,xmax,10)
plt.hist(xx,color='#00ccff',bins=edges,density=True,label='Observed')
plt.plot(xpdf,ypdf,'-',color='#ff0000',lw=2,label='Predicted')
plt.legend()
def wei1(xx,pp):
# parameter calculation by Step-1
ac=np.arange(1,np.min(xx),0.1)
n=len(xx)
y=np.log(-np.log(1-pp))
rr=np.zeros(len(ac),dtype=np.float64)
for i,cc in enumerate(ac):
x=np.log(xx-cc)
rr[i]=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/np.sqrt(n*np.sum(x**2)-(np.sum(x))**2)/np.sqrt(n*np.sum(y**2)-(np.sum(y))**2)
cc=ac[np.argmax(rr)]
x=np.log(xx-cc)
c1=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/(n*np.sum(x**2)-(np.sum(x))**2)
c2=(np.sum(x**2)*np.sum(y)-np.sum(x)*np.sum(x*y))/(n*np.sum(x**2)-(np.sum(x))**2)
kk=c1
aa=np.exp(-c2/c1)
# plot data
yy=cc+aa*(-np.log(1-pp))**(1/kk)
rr=np.corrcoef(xx, yy)
tstr='Weibull (Step-1)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr[0][1])
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa*(-np.log(1-ff))**(1/kk)
xpdf=np.arange(cc+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=kk/aa*coef**(kk-1)*np.exp(-coef**kk)
fig=plt.figure(figsize=(12,4),facecolor='w')
drawqq(xx,yy,tstr,ss)
drawcum(xx,pp,yc,ff,tstr)
drawhist(xx,xpdf,ypdf,tstr)
plt.tight_layout()
fnameF='fig_'+fnameL+'_wei1.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
def wei2(xx,pp):
def func_wei(parameter,xx,pp,cc):
kk=parameter[0]
aa=parameter[1]
est=cc+aa*(-np.log(1-pp))**(1/kk)
residual=xx-est
return residual
# parameter calculation by Step-1
ac=np.arange(1,np.min(xx),0.1)
n=len(xx)
y=np.log(-np.log(1-pp))
rr=np.zeros(len(ac),dtype=np.float64)
for i,cc in enumerate(ac):
x=np.log(xx-cc)
rr[i]=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/np.sqrt(n*np.sum(x**2)-(np.sum(x))**2)/np.sqrt(n*np.sum(y**2)-(np.sum(y))**2)
cc=ac[np.argmax(rr)]
x=np.log(xx-cc)
c1=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/(n*np.sum(x**2)-(np.sum(x))**2)
c2=(np.sum(x**2)*np.sum(y)-np.sum(x)*np.sum(x*y))/(n*np.sum(x**2)-(np.sum(x))**2)
kk=c1
aa=np.exp(-c2/c1)
# parameter calculation by Step-2
parameter0 = [kk,aa]
result = optimize.leastsq(func_wei,parameter0,args=(xx,pp,cc))
kk=result[0][0]
aa=result[0][1]
# plot data
yy=cc+aa*(-np.log(1-pp))**(1/kk)
rr=np.corrcoef(xx, yy)
tstr='Weibull (Step-2)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr[0][1])
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa*(-np.log(1-ff))**(1/kk)
xpdf=np.arange(cc+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=kk/aa*coef**(kk-1)*np.exp(-coef**kk)
fig=plt.figure(figsize=(12,4),facecolor='w')
drawqq(xx,yy,tstr,ss)
drawcum(xx,pp,yc,ff,tstr)
drawhist(xx,xpdf,ypdf,tstr)
plt.tight_layout()
fnameF='fig_'+fnameL+'_wei2.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
def wei3(xx,pp):
def func_wei(parameter,xx,pp):
kk=parameter[0]
cc=parameter[1]
aa=parameter[2]
est=cc+aa*(-np.log(1-pp))**(1/kk)
residual=xx-est
return residual
# parameter calculation
ac=np.arange(1,np.min(xx),0.1)
n=len(xx)
y=np.log(-np.log(1-pp))
rr=np.zeros(len(ac),dtype=np.float64)
for i,cc in enumerate(ac):
x=np.log(xx-cc)
rr[i]=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/np.sqrt(n*np.sum(x**2)-(np.sum(x))**2)/np.sqrt(n*np.sum(y**2)-(np.sum(y))**2)
cc=ac[np.argmax(rr)]
x=np.log(xx-cc)
c1=(n*np.sum(x*y)-np.sum(x)*np.sum(y))/(n*np.sum(x**2)-(np.sum(x))**2)
c2=(np.sum(x**2)*np.sum(y)-np.sum(x)*np.sum(x*y))/(n*np.sum(x**2)-(np.sum(x))**2)
kk=c1
aa=np.exp(-c2/c1)
parameter0 = [kk,cc,aa]
result = optimize.leastsq(func_wei,parameter0,args=(xx,pp))
kk=result[0][0]
cc=result[0][1]
aa=result[0][2]
# plot data
yy=cc+aa*(-np.log(1-pp))**(1/kk)
rr=np.corrcoef(xx, yy)
tstr='Weibull (amended)'
ss='k={0:7.3f}\n'.format(kk)
ss=ss+'c={0:7.3f}\n'.format(cc)
ss=ss+'a={0:7.3f}\n'.format(aa)
ss=ss+'(r={0:.3f})'.format(rr[0][1])
ff=np.arange(0.01,1.00,0.01)
yc=cc+aa*(-np.log(1-ff))**(1/kk)
xpdf=np.arange(cc+1,xmax,1)
coef=(xpdf-cc)/aa
ypdf=kk/aa*coef**(kk-1)*np.exp(-coef**kk)
fig=plt.figure(figsize=(12,4),facecolor='w')
drawqq(xx,yy,tstr,ss)
drawcum(xx,pp,yc,ff,tstr)
drawhist(xx,xpdf,ypdf,tstr)
plt.tight_layout()
fnameF='fig_'+fnameL+'_wei3.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
def main():
global fnameL
args = sys.argv
fnameL=args[1] # location
# data input
fnameR='inp_'+fnameL+'.txt'
data=np.loadtxt(fnameR, skiprows=1)
xx =data[:,2]
xx=np.sort(xx)
# setting protting position
n=len(xx)
ii=np.arange(n)+1
pp=(ii-alpha)/(n+1-2*alpha)
# parameter calculation
wei1(xx,pp)
wei2(xx,pp)
wei3(xx,pp)
#==============
# Execution
#==============
if __name__ == '__main__': main()
正規分布と年降水量の関係
正規分布
統計量$x$が正規分布に従う場合の，確率密度関数，確率分布函数，クオンタイルの数式表現を以下に示す．
確率密度関数
\begin{equation}
f(x)=\frac{1}{\sigma_x\sqrt{2\pi}}\cdot\exp\left[-\frac{1}{2}\left(\frac{x-\mu_x}{\sigma_x}\right)^2\right]
\end{equation}
確率分布関数
\begin{equation}
F(x)=\Phi\left(\frac{x-\mu_x}{\sigma_x}\right)
\qquad \Phi(z)=\frac{1}{\sqrt{2\pi}}\int_{-\infty}^z\exp\left(-\frac{1}{2}t^2\right)dt
\end{equation}
非超過確率pに対する水文統計量xp
\begin{equation}
z=\frac{x-\mu_x}{\sigma_x} \qquad\rightarrow\qquad x=\mu_x+\sigma_x z
\end{equation}
\begin{equation}
x_p=\mu_x+\sigma_x z_p \qquad \text{$z_p$は$p=\Phi(z)$となる$z$の値}
\end{equation}
結果グラフ
母数推定プログラム
py_nor.py
#=========================
# Normal distribution
#=========================
import sys
import numpy as np
from scipy.stats import norm # normal distribution
from scipy import optimize
import matplotlib.pyplot as plt
# Declaration of global variables
fnameL='temp'
xmin=0; xmax=5000; dx=1000
# coefficient for plotting position formula
#alpha=0.0 # Weibull
alpha=0.4 # Cunnane
#alpha=0.5 # Hazen
def drawqq(xx,yy,tstr,ss):
ymin=xmin; ymax=xmax; dy=dx
fsz=12 # font size
plt.subplot(131)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Observed (mm/year)')
plt.ylabel('Predicted (mm/year)')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
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,yy,'o')
plt.plot([xmin,xmax],[ymin,ymax],'-',color='#999999',lw=1)
xs=xmin+(xmax-xmin)*0.03
ys=ymin+(ymax-ymin)*0.97
plt.text(xs,ys,ss,ha='left',va='top',fontname='monospace')
def drawcum(xx,pp,yy,ff,tstr):
ymin=0; ymax=1.0; dy=0.25
fsz=12 # font size
plt.subplot(132)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/year)')
plt.ylabel('Non-exceedance Probability')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
plt.plot(xx,pp,'x',color='#0000ff',markersize=3,label='Observed')
plt.plot(yy,ff,'-',color='#ff0000',lw=1,label='Predicted')
plt.legend()
def drawhist(xx,xpdf,ypdf,tstr):
ymin=0; ymax=0.003; dy=0.001
fsz=12 # font size
plt.subplot(133)
plt.rcParams['font.size']=fsz
plt.rcParams['font.family']='sans-serif'
plt.xlabel('Precipitation (mm/year)')
plt.ylabel('Normalized frequency')
plt.xlim([xmin,xmax])
plt.ylim([ymin,ymax])
plt.xticks(np.arange(xmin,xmax+dx,dx))
plt.yticks(np.arange(ymin,ymax+dy,dy))
plt.grid(color='#999999',linestyle='solid')
plt.title(tstr,loc='left',fontsize=fsz)
edges = range(xmin,xmax,100)
plt.hist(xx,color='#00ccff',bins=edges,density=True,label='Observed')
plt.plot(xpdf,ypdf,'-',color='#ff0000',lw=2,label='Predicted')
plt.legend()
def nor(xx,pp):
def func_nor(parameter,xx,zz):
mu=parameter[0]
sg=parameter[1]
est=mu+sg*zz
residual=xx-est
return residual
# parameter calculation
zz=norm.ppf(pp, loc=0, scale=1)
npara=2; parameter0 = np.ones(npara,dtype=np.float64)
result = optimize.leastsq(func_nor,parameter0,args=(xx,zz))
mu=result[0][0]
sg=result[0][1]
# plot data
yy=mu+sg*zz
rr=np.corrcoef(xx, yy)
tstr='Normal'
ss='m={0:7.3f}\n'.format(mu)
ss=ss+'s={0:7.3f}\n'.format(sg)
ss=ss+'(r={0:.3f})'.format(rr[0][1])
ff=np.arange(0.01,1.00,0.01)
zz=norm.ppf(ff, loc=0, scale=1)
yc=mu+sg*zz
xpdf=np.arange(xmin+1,xmax,1)
coef=(xpdf-mu)/sg
ypdf=1/sg/np.sqrt(2*np.pi)*np.exp(-1/2*(coef)**2)
fig=plt.figure(figsize=(12,4),facecolor='w')
drawqq(xx,yy,tstr,ss)
drawcum(xx,pp,yc,ff,tstr)
drawhist(xx,xpdf,ypdf,tstr)
plt.tight_layout()
fnameF='fig_'+fnameL+'_nor.jpg'
plt.savefig(fnameF, dpi=200, bbox_inches="tight", pad_inches=0.1)
plt.show()
def main():
global fnameL
args = sys.argv
fnameL=args[1] # location
# data input
fnameR='inp_'+fnameL+'.txt'
data=np.loadtxt(fnameR, skiprows=1)
xx =data[:,1]
xx=np.sort(xx)
# setting protting position
n=len(xx)
ii=np.arange(n)+1
pp=(ii-alpha)/(n+1-2*alpha)
# parameter calculation
nor(xx,pp)
#==============
# Execution
#==============
if __name__ == '__main__': main()
以 上