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()

以　上