1．確率雨量って？
最近テレビやネットニュースで「100年に1度の強さの豪雨・大雨」といった"〇〇年に1回の雨"というような表現を見ることが多くなりましたね．日本各地で豪雨による災害が頻発していると感じている方も少なくないでしょう．そこで，今回はそんな雨にまつわる話を水文学的な観点から紹介したいと思います．
さて，この"〇〇年に1回の雨"といった表現ですが，専門的には「確率雨量」と呼ばれたりしています．例えば，100年に1回起こるような雨のことを「100年確率雨量」と呼んだりします．また，〇〇年確率雨量の〇〇年のことを「再現期間」あるいは「リターンピリオド(return period)」と呼びます．
2．確率雨量の重要性
ところで，この"〇〇年に1回の雨"ってどうやって求めているんでしょう？ 当然ですが，100年前や500年前から降水量がずーっと測定されているなんてことはないですから，現時点で100年確率雨量や500年確率雨量を直接知ることは不可能です(50年に1回くらいなら可能かもしれないですが)．
ですが，農業農村工学分野では，なんとかしてこの〇〇年確率雨量を知る必要があるんです． というのもダムや頭首工，河川の堤防などの水利構造物の設計では，例えば，「結構強い雨が降っても決壊したり，水が溢れて大災害が起きないように，まあ，とりあえず100年に1回くらいの雨には耐えられるようにしよう」的な考え方に基づいて設計されているためです．
ここで，疑問持った人もいるのではないでしょうか？ 「じゃあ，500年に1回くらいの強い雨が降ったらどうするの？決壊して大災害になるんじゃない？市民を見捨てる気か！」
「はい．見捨てます!」 ．．．というわけでもないんですが，実際，全ての水害を防ぐことは難しいです． というのも，仮に500年確率雨量に耐えられる水利構造物を作ることが技術的に可能だとしても，建設費用が莫大になります． その費用には税金が充てられます． 完全防災のためには大きな構造物を作ればいいですが，お金がかかるというジレンマを抱えているわけです．加えて，水利構造物の材料のコンクリートの寿命は50〜100年とかだったりするので(コンクリートは専門じゃないでの間違ってたらごめんなさい)，500年に1回の雨が来る前に水利構造物の寿命が来ます．防災には水利構造物の必要性について国民の理解が必要になるわけですね．
3．確率雨量の求め方
さて，〇〇年に1回の雨をどうやって求めるのか？に話を戻します．
大まかな説明ですが，まず過去数十年分の雨量の観測データから，任意の強さの雨が起きる確率(雨量の確率密度関数)を求めます．
仮に，ある場所の雨量の確率密度関数が下図のようになっていたとします．
100年に1回(生起確率1/100=1%)の雨の強さを知りたいときは，その雨量の確率密度関数を使って上側1%点に相当する雨量の値を求めればいいわけです．なお，この生起確率のことは超過確率とも呼びます．
ちなみに渇水側を考える場合，例えば10年に1回の渇水の時の雨量はグラフの赤色の部分の面積が(1-1/10)つまり9/10になるときの横軸の値です．
4．確率雨量の推定手法
筆者が実装経験を持つ推定法は以下の4つになります．本記事では1つ目の対数正規法を紹介したいと思います．
・対数正規法 ・高瀬法(積率法) ・岩井法 ・ガンベル・チャウの方法
なお，本記事の執筆にあたり，以下の書籍を参考にしました．この場を借りて厚く御礼申し上げます． 本記事で公開するプログラムは河村先生の書籍に書かれているもの(Basic)を元に，筆者がVBAで書いたものです． ミスには細心の注意を払っておりますが，完璧を保証するものではありませんので，以下の書籍をお手元にご用意いただくことをお勧めいたします．
河村三郎：土木工学プログラム集1 水文・水理1 実用プログラムとグラフィック表示，森北出版，1984 岩井重久，石黒政儀：応用水文統計学，森北出版，1970 石村園子：やさしく学べるラプラス変換・フーリエ解析 増補版，共立出版，2010
5．対数正規法
対数正規法とは雨量の確率密度関数が対数正規分布に従うという前提に基づいて確率雨量を求める手法です． それでは，対数正規法の具体的な説明に入りたいと思います．
5.1 対数正規分布について
まず，対数正規分布の定義は以下の通りです． 「確率変数yが正規分布に従うとき，EXP(y)が従う分布を対数正規分布という．」 以下のサイトの説明がわかりやすいです．
高校数学の美しい物語〜定期試験から数学オリンピックまで800記事〜：対数正規分布の例と平均，分散 https://mathtrain.jp/lognormal
ここで，
とし，yの確率密度関数をg(y)，xの確率密度関数をf(x)とおくと，
となります． $$ $$
より
したがって，xの確率密度関数f(x)は以下のようになります．
5.2 雨量の確率密度関数
対数正規分布の考え方を雨量に適用してみましょう．
対数正規法では雨量R(確率変数)が対数正規分布に従うと仮定するので，
ちなみに，横軸の雨量Rの対数をとると確率密度関数は以下のような正規分布になります．
このとき，
となります．
5.3 ガウスの誤差関数
ある確率変数xに対するガウスの誤差関数erf(ξ)は次のように定義されます．
「なぜ急に誤差関数が出てくるの？」と思った人も多いでしょうが，これは非超過確率を求めるときに誤差関数を使うからです． 非超過確率とは，確率密度の総和から超過確率を引いたものです． 例えば，100年に1回の雨の生起確率(超過確率)は1/100=1%でしたね．なので，非超過確率とは100%(確率密度の総和)から1%(超過確率)を引いた残りの99%になります．
それでは，正規分布に従う確率変数log10Rに対してガウスの誤差関数の導出をしてみましょう．
確率変数log10Rは，
となりますが，計算上都合が悪いので基準化(標準化)します．
確率変数tの確率密度関数f(t)は，
ここで，誤差関数は以下のようにります．
ここで、
とおくと，
より
ここで、
とおくと，
5.4 超過確率と誤差関数
ここまでの話をまとめると， 雨量Rは対数正規分布に従う． ↓ log10R(雨量Rの常用対数)は正規分布に従う.
↓ log10R(雨量Rの常用対数)を基準化(標準化)すると標準正規分布に従う．
↓ このとき，誤差関数は次式のようになる．
ただし，
超過確率W(ξ)は，誤差関数を用いて次のよう表されます．
つまり，確率密度の総和(1)からある雨量ξに対する誤差関数を引けば，その雨量の再現期間を求めることができます． 逆に，ある再現期間の逆数1/Tと，1/2から誤差関数の1/2を引いた値が一致する雨量RがT年確率雨量になります．
本プログラムでは，計算時間短縮のために，あらかじめ数値積分により再現期間Tとそれに対応するξの値56組分が計算されており(下図のエクセルシートのH列とI列)，それらの間の値は線形補間法により推定します．なお，T<=10では，Tの変化に対するξの変化が大きいため，線形補間法によってξの第1次近似値を与え，試算法によってξを推定しています． 非超過確率の計算にはシンプソンの公式を用いています．
VBAで数値積分 シンプソンの公式(シンプソン法) https://nnhiroba.xsrv.jp/nnTechBlog/Article/Article.php?ID=23
再現期間Tに対するξを線形補間法(+試算法)で求めたら，上記の式で再現期間Tに対するT年確率雨量Rを求めます．
6．計算用シート
計算に用いたデータおよび計算結果をエクセルシートにまとめています．
入力データはエクセルのC列2行目〜21行目(雨量データ)およびH列2行目〜I列57行目(あらかじめ数値積分で求めたTとξ)です． 計算結果はエクセルのM列〜O列に出力されます．
例えば，今回の雨量データが年最大日雨量を20年分集めたものとしたら，100年確率日雨量は145.4 mmとなります．
つまり，水利施設などを作るときの計画基準を100年確率雨量としたとき，145.4 mmの日雨量に耐えうるように設計する必要があるということですね．
計算用エクセルファイルは下記URLからダウンロードできます． ご自由にお使いください． 1水文統計.xlsm
7．おわりに
防災や減災，人命救助などに関心のある人が農業農村工学に興味を持つきっかけとして本記事が貢献できれば幸いです． 高瀬法(積率法)，岩井法，ガンベル・チャウの方法に関しても記事にしていますので，よかったらご覧ください．
【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第1回 対数正規法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第2回 高瀬法(積率法)― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第3回 岩井法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第4回 ガンベル・チャウの方法―
本記事の執筆にあたり，書籍の内容について紹介することを快く承諾いただいた森北出版の皆様と，何より書籍の著者である河村三郎先生には心より感謝いたします．
検索用Key Words： 確率雨量，確率降水量，降雨，降水，水文事象，水文学，水文統計学，降雨事象，水利施設，設計，確率水文量，水文諸量，超過確率，非超過確率，再現期間，リターンピリオド(return period)，T年確率値，事象値，計算方法，解き方，プログラム，マクロ，やり方
8．プログラム
Option Explicit
Function FNE(x As Double)
FNE = Exp(-x ^ 2)
End Function
Function swap(ByRef val1 As Double, ByRef val2 As Double)
'2津の変数を交換する関数
Dim temp As Double
temp = val1
val1 = val2
val2 = temp
End Function
Function BubbleSort(ByRef Array1() As Double, ByVal N As Integer, ByVal SortDirection As String)
'変数の定義
Dim cnt As Integer, i As Integer, j As Integer, k As Integer
'変数の初期化
cnt = 1
'エラー防止
If SortDirection <> "descending" And SortDirection <> "ascending" Then
MsgBox ("引数に誤りがあります．")
End
End If
Select Case SortDirection
'昇順並べ替え
Case Is = "ascending"
k = N - 1
Do While k >= 0
j = -1
i = 1
Do
If Array1(i) > Array1(i + 1) Then
j = i - 1
Call swap(Array1(i), Array1(i + 1))
End If
i = i + 1
Loop While i <= k
k = j
Loop
'降順並べ替え
Case Is = "descending"
k = N - 1
Do While k >= 0
j = -1
i = 1
Do
If Array1(i) < Array1(i + 1) Then
j = i - 1
Call swap(Array1(i), Array1(i + 1))
End If
i = i + 1
Loop While i <= k
k = j
Loop
End Select
End Function
Sub Log_Normal_Distribution()
'最終更新日：2020/7/19
'川村三郎：土木工学プログラム集 水文・水理1，森北出版，1984に準拠
'*********************プログラム起動時のメッセージを出力*********************
Debug.Print "対数正規法"
Debug.Print Format(Date) & " " & Hour(Time) & "時" & Minute(Time) & "分" & Second(Time) & "秒(S.K)"; " 計算開始"
Worksheets("水文統計").Activate
Application.ScreenUpdating = False
'*********************変数の定義*********************
Dim N As Integer 'データ数Number
Dim TS As Double '再現期間の初期値Time Start(年)
Dim TE As Double '再現期間の最終値Time End(年)
Dim SY As Double '再現期間のキザミStep Year(年)
Dim m As Integer 'ξデータの数
Dim slct As Integer '入力データの種類(雨量or流出量)の判別用
Dim R() As Double '降水量データを格納する配列
'普遍定数など
Dim p As Double '円周率
Dim PL As Double 'ln10
'集計用変数
Dim X5 As Double '自然対数をとったデータの合計値
Dim x9 As Double '自然対数をとったデータの偏差二乗和
Dim k As Integer, i As Integer, cnt As Integer
Dim X0 As Double 'EXP(自然対数をとったデータの平均値)
'シンプソン公式および試算法などに関わる変数
Dim E2 As Double 'ξ
Dim SG As Double '自然対数をとったデータの標準偏差
Dim A1 As Double '区間[A1,B1]に関して数値積分する
Dim C1 As Double '補正係数
Dim B1 As Double '区間[A1,B1]に関して数値積分する
Dim ND As Integer '積分範囲の分割数
Dim W As Double '超過確率
Dim W1 As Double '超過確率の第一次近似値
Dim E3 As Double 'ξの補正値(試算法で使用)
Dim x1 As Double 'シンプソン公式第二項
Dim x2 As Double 'シンプソン公式第三項
Dim S As Double '(π)^0.5/2/erf(ξ)，シンプソン法を使用した数値積分で求める
Dim D As Double '数値積分における小区間の長さ
'確率雨量全般に関する変数
Dim x As Double '式(1.2.4)の左辺のx，つまり，T年確率降水量
Dim t As Double '再現期間
'対数正規法に固有の変数
Dim K1 As Double '式(1.2.4)の第二項
Dim K2 As Double '式(1.2.4)の第一項
Dim K3 As Double 'log10(x) 式(1.2.4)の右辺
'*********************定数などの定義*********************
m = 56 'ξデータの数
X5 = 0 '自然対数をとったデータの合計値
x9 = 0 '自然対数をとったデータの偏差二乗和
p = 3.14159 '円周率
PL = 2.30259 'ln10(10の自然対数)
'*********************データ数および計算条件の入力*********************
N = InputBox("降雨または流量のデータ数を入力してください．")
TS = InputBox("再現期間の初期値(年)を入力してください．")
TE = InputBox("再現期間の最終値(年)を入力してください．")
SY = InputBox("再現期間のキザミ(年)を入力してください．")
'*********************入力データの種類(雨量or流出量)の入力と計算条件の出力*********************
Do
slct = InputBox("雨量なら1を,流出量なら2を入力してください．")
If slct = 1 Then
Debug.Print "R：降雨(mm)"
ElseIf slct = 2 Then
Debug.Print "R:流出量(m3/s)"
End If
Loop Until slct = 1 Or slct = 2
Debug.Print "再現期間の初期値(年)TS="; TS; "(year)"
Debug.Print "再現期間の最終値(年)TE="; TE; "(year)"
Debug.Print "再現期間のキザミ(年)SY="; SY; "(year)"
Debug.Print ""
Debug.Print "ORIGINAL DATA"
Debug.Print "データ数N="; N
'*********************入力データを降順に並び替える*********************
Cells(1, 4) = "R (mm) 降順"
ReDim R(1 To N)
'降水量データを取得
For i = 1 To N
R(i) = Cells(i + 1, 3)
Next i
'バブルソート
Call BubbleSort(R, N, "descending")
'降順に並び替えたデータをワークシートに表示
For i = 1 To N
Cells(i + 1, 4) = R(i)
Next i
'*********************入力データを計算用に加工する*********************
'BasicのプログラムのLOG(R),SUM(log(R)),exp(SUM/N) 430〜480の部分
'データの自然対数をとる
Cells(1, 5) = "lnR R1"
Cells(N + 3, 1) = "合計ΣlnR x5"
Cells(N + 3, 2) = 0
For k = 1 To N
Cells(k + 1, 5) = Log(Cells(k + 1, 4))
Cells(N + 3, 2) = Cells(N + 3, 2) + Cells(k + 1, 5)
Next k
X5 = Cells(N + 3, 2) '自然対数をとったデータの合計値
X0 = Exp(Cells(N + 3, 2) / N) 'EXP(自然対数をとったデータの平均値)
Cells(N + 4, 1) = "EXP(Σ(logR)/N) x0"
Cells(N + 4, 2) = X0
Debug.Print "x0="; X0
'*********************log(R)の平均値および標準偏差を出力 500〜590*********************
Cells(N + 5, 1) = "ln(X0) / PL x7"
Cells(N + 5, 2) = Log(X0) / PL
Cells(N + 6, 1) = "Σx8 x9"
Cells(N + 6, 2) = 0
Cells(1, 6) = "log10R x6"
ActiveCell.Characters(start:=4, Length:=2).Font.Subscript = True
Cells(1, 7) = "(log10R-log10(x0))2 x8"
ActiveCell.Characters(start:=5, Length:=2).Font.Subscript = True
ActiveCell.Characters(start:=12, Length:=2).Font.Subscript = True
ActiveCell.Characters(start:=19, Length:=1).Font.Superscript = True
'*********************データの常用対数と偏差，標準偏差を求める*********************
Cells(N + 7, 1) = "σ0"
ActiveCell.Characters(start:=2, Length:=1).Font.Subscript = True
For k = 1 To N
Cells(k + 1, 6) = Cells(k + 1, 5) / PL
Cells(k + 1, 7) = (Cells(k + 1, 6) - Cells(N + 5, 2)) ^ 2
Cells(N + 6, 2) = Cells(N + 6, 2) + Cells(k + 1, 7)
Next k
Cells(N + 7, 2) = Sqr(Cells(N + 6, 2) / N)
SG = Cells(N + 7, 2)
Debug.Print "standard Deviation"
Debug.Print "sigma="; Cells(N + 7, 2)
Debug.Print ""
'*********************計算結果用ラベル*********************
Cells(1, 13) = "対数正規法"
Cells(2, 13) = "T(year)"
Cells(2, 14) = "ξ"
Cells(2, 15) = "R(mm)"
'*********************超過確率に対するξの値(E2)，T年確率雨量xを計算する*********************
'COMP. OF KUSI線形補間法
'T>10では，線形補間法によってξを求める．再現期間Tに対する超過確率W1はシンプソン公式による数値積分から求める．
'T<=10では，線形補間法によってξの第1次近似値を与え，シンプソン公式により超過確率の第一近似値W1を得る．
'(続き)第一次近似値ξ，W1を使って，真の超過確率Wとの超過確率の近似値W1の差が十分に小さくなるまで試算法によって計算し，ξの値を決定する．
'(続き)T<=10では，xの変化量に対してerf(x)の変化量つまり，再現期間および超過確率が大きく変化するので試算法を用いる
cnt = 0
For t = TS To TE Step SY
'********************************************************************************************
'*********************T>10のとき，シンプソン法による数値積分でξを求める*********************
'********************************************************************************************
'*********************線形補間によりξの値(E2)を求める*********************
k = 0
Do
k = k + 1
Loop Until t < Cells(k + 2, 8)
k = k - 1
E2 = (Cells(k + 3, 9) - Cells(k + 2, 9)) / (Cells(k + 3, 8) - Cells(k + 2, 8)) _
* (t - Cells(k + 2, 8)) + Cells(k + 2, 9)
'予めTに対するξの値が得られているときは，数値積分と試算法はスルー
'TS,TE,SYを実数型で定義しているため「if t=Cells(k+2,8) then」だと正しく処理できない
If Abs(t - Cells(k + 2, 8)) < 0.0000001 Then
GoTo PerfectMatchingKusi
End If
'*********************数値積分*********************
A1 = 0 '区間[A1,B1]に関して数値積分する
C1 = 0.5 '補正係数
B1 = E2 '区間[A1,B1]に関して数値積分する
ND = 50 '積分範囲の分割数
W = 1 / t '超過確率
D = (B1 - A1) / (2 * ND) '区間[0, E2]をND個に分割(分母の2はシンプソン公式の第二項のxの値にするためのもの)
S = FNE(A1) + FNE(B1) 'シンプソン公式のf(a),f(b)に相当するFNE(A1),FNE(B1)をまずSに入れる
x1 = D 'シンプソン公式第二項
x2 = D + x1 'シンプソン公式第三項
For k = 1 To ND - 1
S = S + 4 * FNE(x1) + 2 * FNE(x2)
x1 = x2 + D
x2 = x1 + D
Next k
S = S + 4 * FNE(x1)
S = (S * D / 6) * 2 'erf(ξ)の積分範囲は[-ξ,ξ]だが，数値積分では[0,ξ]で積分しており，確率変数xはN(0,1)に従うので
W1 = 0.5 - S / Sqr(4 * Atn(1)) '式(1.2.8)
'*************************************************************************************************************************************
'*********************T<=10のときは，シンプソン公式による数値積分に加えて試算法も行うことで，最適なξの値を求める*********************
'*************************************************************************************************************************************
If t <= 10 Then
'真の超過確率Wとの超過確率の近似値W1の差が十分に小さくなるまで繰り返す
Do Until Abs(W - W1) <= 0.00005
'*********************補正係数の計算*********************
E3 = E2 - C1 * (W - W1) / W '第一次近似値から相対誤差の50%を引く
E3 = (E3 + E2) / 2
E2 = E3
B1 = E3
'*********************数値積分*********************
D = (B1 - A1) / (2 * ND) '区間[0, E2]を2*50に分割
S = FNE(A1) + FNE(B1) 'シンプソン公式のf(a),f(b)に相当するFNE(A1),FNE(B1)をまずSに入れる
x1 = D
x2 = D + x1
For k = 1 To ND - 1
S = S + 4 * FNE(x1) + 2 * FNE(x2)
x1 = x2 + D
x2 = x1 + D
Next k
S = S + 4 * FNE(x1)
S = (S * D / 6) * 2
'第二次近似値を求める
W1 = 0.5 - S / Sqr(4 * Atn(1))
Loop
End If
'**************************************************************************
'*********************T年確率雨量xを求める(結果の出力)*********************
'**************************************************************************
PerfectMatchingKusi:
Debug.Print "T="; t; "年確率雨量のξ()は"; E2
K1 = SG * E2 * Sqr(2) 'σ0*√2*ξ 式(1.2.4)の第二項
K2 = Log(X0) / PL 'log10(x0) 式(1.2.4)の第一項
K3 = K1 + K2 'log10(x) 式(1.2.4)の右辺
x = Exp(PL * K3) '式(1.2.4)の左辺のx，つまり，T年確率降水量
Cells(3 + cnt, 13) = t
Cells(3 + cnt, 14) = E2
Cells(3 + cnt, 15) = x
cnt = cnt + 1
Next t
Debug.Print Format(Date) & " " & Hour(Time) & "時" & Minute(Time) & "分" & Second(Time) & "秒(S.K)"; " 計算終了"
Application.ScreenUpdating = True
End Sub
