はじめに
ここ最近、大雨のニュースを見るたびに「〇〇年に一度の大雨」という表現を目にします。これは「千年に一人の美少女」や「百年に一度の逸材」といったレトリックではなく、データに基づいて統計的に算出された根拠のある数字です。
気象庁や自治体では観測結果をもとにこのような量を算出し、発表することで防災に役立てています。このシリーズでは、算出方法を確認し、Pythonを使って、実装・計算してみたいと思います。
なお、ここで紹介する手法は、実際に公的機関が発表しているものとは異なる場合が多分にありますので、十分にご注意ください。
今回は、前回に引き続き、当てはめた分布関数から「〇〇年に一度の大雨」の値を算出してみたいと思います。
macOS Monterey 12.6.2, Python 3.9.15, SciPy 1.11.2
Pythonでやってみる気象データ解析の技術書を販売中
¥2,500 → ¥1,500 今なら40%OFF!!
気象の勉強をされている方必読！
気象予報士でもある著者がわかりやすく解説
プログラムを公開しているのでお手元で試せます！
「○○年に一度の大雨」の計算方法とは？
気象庁は観測結果をもとに、各地点について「〇〇年に一度の降水量」を算出し、発表しています。
ある期間内に1回起こると考えられる降水量のことを確率降水量といいます。気象庁の解説ページでは、確率降水量の推定方法が解説されていますので、これに基づいて「〇〇年に一度の大雨」の値を計算してみたいと思います。
大まかには次のような流れになっています。
- 年最大日降水量のヒストグラムを作成する
- 分布関数を当てはめる
- 分布関数の当てはまり具合を確認する
- 当てはめた分布関数から確率降水量を算出する
今回は、上記3. をとばして、4.の「当てはめた分布関数から確率降水量を算出する」を実際にやってみましょう！
再現年と確率降水量
気象庁の解説ページにもあるように、再現年 $T$ は
$$ T = \frac{1}{1 – F(x; \theta)} $$
で与えられます。ここで、$F(x; \theta)$ は確率分布関数で、当てはめた確率密度関数 $f(x; \theta)$ と
$$ F(x;\theta) = \int_{-\infty}^{x} f(x’;\theta)\mathrm{d}x’ $$
という関係があります。これより、再現年 $T$ の確率降水量 $x$ (例えば「 $T$ 年に1度の大雨 $x$ mm」) は
$$ x = F^{-1}\left(1-\frac{1}{T};\theta \right) $$
となります。ここで、$F^{-1}(y; \theta)$ は $F(x; \theta)$ の逆関数です。
3つの分布関数について、具体的に再現年から確率降水量を求める式を書くと次のようになります。
グンベル分布の場合
$$ x = \mu – \sigma\ln\left[-\ln \left(1-\frac{1}{T}\right)\right] $$
一般化極値分布の場合
$$ x = \mu + \frac{\sigma}{c}\left\{1-\left[-\ln \left(1-\frac{1}{T}\right)\right]^c\right\} $$
平方根指数型最大値分布の場合
グンベル分布や一般化極値分布のような初等関数の組み合わせで書くことはできません。以下求める手順を書きます。まず、
$$ x = \frac{t^2}{b} $$
として、新たに $t$ を定義します。代入して整理すると、
$$ \ln(1+t) – t – \ln\left[-\frac{1}{a}\ln\left(1-\frac{1}{T}\right)\right] = 0 $$
となりますので、これを $t$ について解いて、上の式に代入すると、確率降水量 $x$ が得られます。
実装と実行
プログラムの例
前回のプログラムを改良して、確率降水量を求める部分を追加しましょう。追加後のものを以下に示します。
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as ss
from scipy.optimize import minimize, root_scalar
#-----------------------------------------------------------------------------
# --- 平方根指数型最大値分布の関数
def sqet(x, a, b):
f = a*b*0.5*np.exp(-np.sqrt(b*x)-a*(1.0+np.sqrt(b*x))*np.exp(-np.sqrt(b*x)))
return f
def sqet_cdf(x, a, b):
F = np.exp(-a*(1.0+np.sqrt(b*x))*np.exp(-np.sqrt(b*x)))
return F
# SQETの対数尤度関数 (にマイナスをつけたもの)
def mL_sqet(b, data):
n = len(data)
sqbx = np.sqrt(b*data)
a = (np.sum(sqbx) -2.0*n) / np.sum(b*data*np.exp(-sqbx))
L = n*np.log(a) + n*np.log(b) -n*np.log(2) - np.sum(sqbx)
L = L - a*(np.sum(np.exp(-sqbx)) + np.sum(sqbx*np.exp(-sqbx)))
return -L
# bからaを求める式
def b2a_sqet(b, data):
n = len(data)
sqbx = np.sqrt(b*data)
a = (np.sum(sqbx) -2.0*n) / np.sum(b*data*np.exp(-sqbx))
return a
#------------------------------------------------------------------------------
# データの読み込み
year = []
annual_max_daily_prcp = []
i = 0
with open('data.csv', 'r') as fin:
for line in fin:
if i > 1:
data = line.split(',')
year.append(int(data[0]))
annual_max_daily_prcp.append(float(data[1].replace(']', '')))
i += 1
# ndarrayに変換
year = np.array(year)
annual_max_daily_prcp = np.array(annual_max_daily_prcp)
#------------------------------------------------------------------------------
# --- プロット
# 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'histogram-and-fitted-pdf-mod.png'
# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.80, 0.80, 0.15])
ax02 = fig.add_axes([0.10, 0.55, 0.80, 0.15])
ax03 = fig.add_axes([0.10, 0.35, 0.80, 0.15])
ax04 = fig.add_axes([0.10, 0.10, 0.80, 0.15])
# --- plot 1
# 年最大日降水量を西暦毎にバーグラフで表示
ax01.bar(year, annual_max_daily_prcp)
ax01.set_xlabel('year')
ax01.set_ylabel('annual maximum\\ndaily precipitation (mm)')
# グラフ上に最大値、最小値、平均値を表示
fig.text(0.01, 0.9, f'max: {annual_max_daily_prcp.max()} mm ({year[np.argmax(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.8, f'min: {annual_max_daily_prcp.min()} mm ({year[np.argmin(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.7, f'mean: {annual_max_daily_prcp.mean():.1f} mm', transform=ax01.transAxes)
# --- plot 2 と plot 3
# Gumbel分布を当てはめる
loc_gumbel, scale_gumbel = ss.gumbel_r.fit(annual_max_daily_prcp, method='MLE')
# 一般化極値分布を当てはめる
c_gev, loc_gev, scale_gev = ss.genextreme.fit(annual_max_daily_prcp, method='MLE')
# 平方根指数型最大値分布を当てはめる
res = minimize(mL_sqet, x0=[1.0], args=annual_max_daily_prcp,
bounds=[(0, None)], method='Nelder-Mead')
b_sqet = res.x[0]
a_sqet = b2a_sqet(b_sqet, annual_max_daily_prcp)
gbl = ss.gumbel_r(loc=loc_gumbel, scale=scale_gumbel)
gev = ss.genextreme(c=c_gev, loc=loc_gev, scale=scale_gev)
# bin, x のためのリストを作成する
bins = [10*(i+1) for i in range(int(annual_max_daily_prcp.max()/10) + 1)]
x = np.linspace(0.0, 10.0*(int(annual_max_daily_prcp.max()/10) + 1), 20*int(annual_max_daily_prcp.max()/10) + 1)
# ヒストグラムを表示
ax02.hist(annual_max_daily_prcp, bins=bins, density=True)
# 当てはめた分布をプロット
ax02.plot(x, gbl.pdf(x), label=f"scale={scale_gumbel:.2e}, loc={loc_gumbel:.2e}")
ax02.plot(x, gev.pdf(x), label=f"scale={scale_gev:.2e}, loc={loc_gev:.2e}, c={c_gev:.2e}")
ax02.plot(x, sqet(x, a_sqet, b_sqet), label=f"a={a_sqet:.2e}, b={b_sqet:.2e}")
ax02.legend()
# 累積ヒストグラムを表示
ax03.hist(annual_max_daily_prcp, bins=bins, cumulative=True, density=True)
ax03.plot(x, gbl.cdf(x))
ax03.plot(x, gev.cdf(x))
ax03.plot(x, sqet_cdf(x, a_sqet, b_sqet))
ax03.set_xlabel('annual maximum daily precipitation (mm)')
# --- plot 4
# 再現期間と確率降水量のプロット
ax04.plot(x, 1.0/(1.0-gbl.cdf(x)), color='red')
ax04.plot(x, 1.0/(1.0-gev.cdf(x)), color='blue')
ax04.plot(x, 1.0/(1.0-sqet_cdf(x, a_sqet, b_sqet)), color='green')
ax04.set_yscale('log')
ax04.set_ylim(1, 150)
ax04.set_xlabel('probable maximum precipitation (mm)')
ax04.set_ylabel('return period (years)')
# 確率降水量の値を求める式
# Gumbel分布
def pp_gumbel(scale, loc, T):
p = 1.0 - 1.0/T
x = loc - scale*np.log(-np.log(p))
return x
# 一般化極値分布
def pp_gev(c, scale, loc, T):
p = 1.0 - 1.0/T
x = loc + scale/c*(1.0-(-np.log(p))**c)
return x
# 平方根指数型最大値分布
def pp_sqet(a, b, T):
p = 1.0 - 1.0/T
def f(t):
return np.log(1.0 + t) - t - np.log(-1.0/a*np.log(p))
sol = root_scalar(f=f, method='newton', x0=100.0)
x = sol.root**2/b
return x
# 再現期間50年の確率降水量を算出し、線を引く
ax04.axhline(50.0, color='gray', linestyle='dashed')
ax04.axvline([pp_gumbel(scale_gumbel, loc_gumbel, 50.0)], color='salmon', linestyle='dashed')
ax04.axvline([pp_gev(c_gev, scale_gev, loc_gev, 50.0)], color='skyblue', linestyle='dashed')
ax04.axvline([pp_sqet(a_sqet, b_sqet, 50.0)], color='lightgreen', linestyle='dashed')
# グラフ上に再現期間50年の確率降水量の値を書く
fig.text(0.01, 0.9, f'Gumbel: {pp_gumbel(scale_gumbel, loc_gumbel, 50.0):.0f} mm', color='red', transform=ax04.transAxes)
fig.text(0.01, 0.8, f'GEV: {pp_gev(c_gev, scale_gev, loc_gev, 50.0):.0f} mm', color='blue', transform=ax04.transAxes)
fig.text(0.01, 0.7, f'SQET: {pp_sqet(a_sqet, b_sqet, 50.0):.0f} mm', color='green', transform=ax04.transAxes)
#------------------------------------------------------------------------------
# --- 図の保存
plt.savefig(out_fig_path, transparent=False)
このコードを fit-pdfs-mod.pyとして、data.csv というデータファイルのあるディレクトリ (ここでは、~/Desktop/labcode/python/probable-precipitationとします。) に保存します。
データファイルは、二行目以降が第一コラムに年、第二コラムが年最大日降水量となっているものを想定しています。作り方は、以前の記事で紹介しています。参考にしてください。
プログラムを実行する
ターミナルを開き、
$ cd ~/Desktop/labcode/python/probable-precipitationと入力し、ディレクトリを移動します（ $マークは無視してください）。次のようにして先ほどのコードを実行しましょう。
$ python fit-pdfs-mod.pyとします。histogram-and-fitted-pdf-mod.png というファイルができて次のようになっていればひとまず成功です！
四段目に当てはめた分布関数から求められる再現期間と確率降水量の関係がプロットされています。
グラフの左上には、再現年50年とした時の、それぞれの分布から計算される確率降水量が書かれています。グンベル分布では176 mm、一般化極値分布では184 mm、平方根指数型最大値分布では196mmとなっています。分布関数の選択によって確率降水量の推定値に1割程度の違いが見られます。
本来であれば、この後、関数の当てはまり具合 (適合度) や安定度 (極端なデータに引っ張られていないか) をチェックする処理が行われます。
ちなみに、気象庁が推定したものによると、大阪市の分布はGumbel分布が最も当てはまりが良く、再現年50年の値は173mmとなっていて、今回推定したものと近い値となっています。
値に若干の違いが見られた理由としては、使用したデータの期間が異なることがまず第一に考えられます。
コードの解説
上で紹介したコードの解説をします。前半部分は以前のコードと同じなので、追加した部分をメインに解説します。
import matplotlib.pyplot as plt
import numpy as np
import scipy.stats as ss
from scipy.optimize import minimize, root_scalar
まず、必要なモジュールをインポートします。 当てはめた平方根指数型最大値分布のパラメータから確率降水量を導くために使用する、root_scalar を追加でインポートしています。
#-----------------------------------------------------------------------------
# --- 平方根指数型最大値分布の関数
def sqet(x, a, b):
f = a*b*0.5*np.exp(-np.sqrt(b*x)-a*(1.0+np.sqrt(b*x))*np.exp(-np.sqrt(b*x)))
return f
def sqet_cdf(x, a, b):
F = np.exp(-a*(1.0+np.sqrt(b*x))*np.exp(-np.sqrt(b*x)))
return F
# SQETの対数尤度関数 (にマイナスをつけたもの)
def mL_sqet(b, data):
n = len(data)
sqbx = np.sqrt(b*data)
a = (np.sum(sqbx) -2.0*n) / np.sum(b*data*np.exp(-sqbx))
L = n*np.log(a) + n*np.log(b) -n*np.log(2) - np.sum(sqbx)
L = L - a*(np.sum(np.exp(-sqbx)) + np.sum(sqbx*np.exp(-sqbx)))
return -L
# bからaを求める式
def b2a_sqet(b, data):
n = len(data)
sqbx = np.sqrt(b*data)
a = (np.sum(sqbx) -2.0*n) / np.sum(b*data*np.exp(-sqbx))
return a
#------------------------------------------------------------------------------
# データの読み込み
year = []
annual_max_daily_prcp = []
i = 0
with open('data.csv', 'r') as fin:
for line in fin:
if i > 1:
data = line.split(',')
year.append(int(data[0]))
annual_max_daily_prcp.append(float(data[1].replace(']', '')))
i += 1
# ndarrayに変換
year = np.array(year)
annual_max_daily_prcp = np.array(annual_max_daily_prcp)
#------------------------------------------------------------------------------
# --- プロット
# 全般の設定
fig = plt.figure(figsize=(210/25.4, 294/25.4), dpi=100)
out_fig_path = 'histogram-and-fitted-pdf-mod.png'
# プロット枠の設定
ax01 = fig.add_axes([0.10, 0.80, 0.80, 0.15])
ax02 = fig.add_axes([0.10, 0.55, 0.80, 0.15])
ax03 = fig.add_axes([0.10, 0.35, 0.80, 0.15])
ax04 = fig.add_axes([0.10, 0.10, 0.80, 0.15])
# --- plot 1
# 年最大日降水量を西暦毎にバーグラフで表示
ax01.bar(year, annual_max_daily_prcp)
ax01.set_xlabel('year')
ax01.set_ylabel('annual maximum\\ndaily precipitation (mm)')
# グラフ上に最大値、最小値、平均値を表示
fig.text(0.01, 0.9, f'max: {annual_max_daily_prcp.max()} mm ({year[np.argmax(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.8, f'min: {annual_max_daily_prcp.min()} mm ({year[np.argmin(annual_max_daily_prcp)]})', transform=ax01.transAxes)
fig.text(0.01, 0.7, f'mean: {annual_max_daily_prcp.mean():.1f} mm', transform=ax01.transAxes)
# --- plot 2 と plot 3
# Gumbel分布を当てはめる
loc_gumbel, scale_gumbel = ss.gumbel_r.fit(annual_max_daily_prcp, method='MLE')
# 一般化極値分布を当てはめる
c_gev, loc_gev, scale_gev = ss.genextreme.fit(annual_max_daily_prcp, method='MLE')
# 平方根指数型最大値分布を当てはめる
res = minimize(mL_sqet, x0=[1.0], args=annual_max_daily_prcp,
bounds=[(0, None)], method='Nelder-Mead')
b_sqet = res.x[0]
a_sqet = b2a_sqet(b_sqet, annual_max_daily_prcp)
gbl = ss.gumbel_r(loc=loc_gumbel, scale=scale_gumbel)
gev = ss.genextreme(c=c_gev, loc=loc_gev, scale=scale_gev)
# bin, x のためのリストを作成する
bins = [10*(i+1) for i in range(int(annual_max_daily_prcp.max()/10) + 1)]
x = np.linspace(0.0, 10.0*(int(annual_max_daily_prcp.max()/10) + 1), 20*int(annual_max_daily_prcp.max()/10) + 1)
# ヒストグラムを表示
ax02.hist(annual_max_daily_prcp, bins=bins, density=True)
# 当てはめた分布をプロット
ax02.plot(x, gbl.pdf(x), label=f"scale={scale_gumbel:.2e}, loc={loc_gumbel:.2e}")
ax02.plot(x, gev.pdf(x), label=f"scale={scale_gev:.2e}, loc={loc_gev:.2e}, c={c_gev:.2e}")
ax02.plot(x, sqet(x, a_sqet, b_sqet), label=f"a={a_sqet:.2e}, b={b_sqet:.2e}")
ax02.legend()
# 累積ヒストグラムを表示
ax03.hist(annual_max_daily_prcp, bins=bins, cumulative=True, density=True)
ax03.plot(x, gbl.cdf(x))
ax03.plot(x, gev.cdf(x))
ax03.plot(x, sqet_cdf(x, a_sqet, b_sqet))
ax03.set_xlabel('annual maximum daily precipitation (mm)')
データの読み込み、関数の当てはめ、その後のプロットは前回と同じです。
ただし、グラフをもう一つ追加するために、プロット枠の設定の部分は変えています。
# --- plot 4
# 再現期間と確率降水量のプロット
ax04.plot(x, 1.0/(1.0-gbl.cdf(x)), color='red')
ax04.plot(x, 1.0/(1.0-gev.cdf(x)), color='blue')
ax04.plot(x, 1.0/(1.0-sqet_cdf(x, a_sqet, b_sqet)), color='green')
ax04.set_yscale('log')
ax04.set_ylim(1, 150)
ax04.set_xlabel('probable maximum precipitation (mm)')
ax04.set_ylabel('return period (years)')
ここからが、主に新たに追加した部分です。確率分布関数から再現期間を求める式を使って、それぞれの分布関数に対して降水量と再現期間の関係をプロットします。
# 確率降水量の値を求める式
# Gumbel分布
def pp_gumbel(scale, loc, T):
p = 1.0 - 1.0/T
x = loc - scale*np.log(-np.log(p))
return x
# 一般化極値分布
def pp_gev(c, scale, loc, T):
p = 1.0 - 1.0/T
x = loc + scale/c*(1.0-(-np.log(p))**c)
return x
# 平方根指数型最大値分布
def pp_sqet(a, b, T):
p = 1.0 - 1.0/T
def f(t):
return np.log(1.0 + t) - t - np.log(-1.0/a*np.log(p))
sol = root_scalar(f=f, method='newton', x0=100.0)
x = sol.root**2/b
return x
この部分は、与えられた再現期間に対して、確率降水量がいくつかを算出する式を定義しているところです。グンベル分布と一般化極値分布は解析的な式で書けますが、平方根指数型最大値分布は数値計算で求める必要があります。scipy.optimize の root_scalar を使って、Newton法で解を求めています。簡単に解説しておくと、まず、$f(t)$ という解を求めたい方程式を定義しておいて、root_scalarの第一引数に与えると、$f(t) = 0$ の解を求めてくれるというものです。
# 再現期間50年の確率降水量を算出し、線を引く
ax04.axhline(50.0, color='gray', linestyle='dashed')
ax04.axvline([pp_gumbel(scale_gumbel, loc_gumbel, 50.0)], color='salmon', linestyle='dashed')
ax04.axvline([pp_gev(c_gev, scale_gev, loc_gev, 50.0)], color='skyblue', linestyle='dashed')
ax04.axvline([pp_sqet(a_sqet, b_sqet, 50.0)], color='lightgreen', linestyle='dashed')
# グラフ上に再現期間50年の確率降水量の値を書く
fig.text(0.01, 0.9, f'Gumbel: {pp_gumbel(scale_gumbel, loc_gumbel, 50.0):.0f} mm', color='red', transform=ax04.transAxes)
fig.text(0.01, 0.8, f'GEV: {pp_gev(c_gev, scale_gev, loc_gev, 50.0):.0f} mm', color='blue', transform=ax04.transAxes)
fig.text(0.01, 0.7, f'SQET: {pp_sqet(a_sqet, b_sqet, 50.0):.0f} mm', color='green', transform=ax04.transAxes)
この部分は、再現期間を50年として確率降水量を実際に求めて、線をひいたり、グラフの左上に書いたりしている部分です。
上で定義した関数に、当てはめたパラメータと再現期間 50.0 を渡して確率降水量を求めています。
最後に
今回は、「〇〇年に一度の大雨」について理解するために、年最大日降水量のヒストグラムに分布関数を当てはめて、そのパラメータから実際に50年に一度の大雨の値を計算してみました。公的機関等では、実この後、適合度や安定度の評価をして、公表されているようです。
身の回りのデータセットについても、同様の統計処理が適用できる場合があります (「100年に一度の逸材」については、どのようなデータセットがあれば適用できるのでしょうか？) 。色々試してみてください。
参考文献
本記事を執筆するに当たり、次の資料・文献を参考にしました。
- 気象庁、確率降水量の推定方法：https://www.data.jma.go.jp/cpdinfo/riskmap/cal_qt.html
- 椎葉充晴、立川康人、市川温、『例題で学ぶ水文学』、森北出版
- M. Evans, N. Hastings, B. Peacock, “Statistical Distributions, Third Edition”, Wiley-Interscience.
- 星清、『水文統計解析』、開発土木研究月報 No.540、1998年5月：https://thesis.ceri.go.jp/db/files/0005005050.pdf
Pythonでやってみる気象データ解析の技術書を販売中
¥2,500 → ¥1,500 今なら40%OFF!!
気象の勉強をされている方必読！
気象予報士でもある著者がわかりやすく解説
プログラムを公開しているのでお手元で試せます！
