【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第2回 高瀬法(積率法)―
水文・水質・気象
VBA
シミュレーション
統計
確率密度関数
予測
1．高瀬法(積率法)による確率雨量の推定
確率雨量シリーズの第1回では，対数正規法で確率雨量を求めました． 本記事では2つ目の高瀬法(積率法)を紹介したいと思います． ・対数正規法 ・高瀬法(積率法) ・岩井法 ・ガンベル・チャウの方法
【確率雨量】シリーズ一覧 第1回(対数積法)で確率雨量などについて詳細に説明しているので，よかったらご覧ください．
【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第1回 対数正規法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第2回 高瀬法(積率法)― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第3回 岩井法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第4回 ガンベル・チャウの方法―
2．確率雨量について(復習)
前回の記事でも紹介した「確率雨量」についておさらいしたいと思います．
"〇〇年に1回の雨"という表現は，専門的には「確率雨量」と呼ばれたりしています．例えば，100年に1回起こるような雨のことを「100年確率雨量」と呼んだりします．また，〇〇年確率雨量の〇〇年のことを「再現期間」と呼びます．
単純に考えれば，100年に1回の雨が起きる確率(生起確率)は1/100=1%で，この確率のことを「超過確率」とも呼びます．
3．高瀬法(積率法)
高瀬法は確率雨量シリーズ第1回の対数正規法と基本的な考え方は大きく変わりません． ただし，雨量\(R\)と\( \xi \)の関係式が異なります．具体的に次式のようになります． なぜ，この式の形になるのかは勉強中なので，また後日記事にまとめたいと思います．
\[
\log _{10}R=\log _{10}R_{0}+\left( \dfrac {\sigma _{0}}{\sigma _{\xi}}\right) \xi
\]
\( \sigma \xi \)は，データ数\(N\)によって決まる値で，エクセルのJ列2行目〜K列16行目に示したような値をとります． 中間の値はラグランジュ補間法によって内挿します．
\( \xi \)は対数正規法と同じ方法で計算し，上記の式から確率雨量Rを求めます．
4．計算用シートとプログラム
計算に用いたデータおよび計算結果をエクセルシートにまとめています．
入力データはエクセルのC列2行目〜21行目(雨量データ)，H列2行目〜I列57行目(あらかじめ数値積分で求めた\(T\)と\( \xi \)およびJ列2行目〜K列16行目(\( \sigma \xi \)は，データ数\( N \))です． 計算結果はエクセルのP列〜R列に出力されます．
高瀬法(積率法)では，対数正規分布よりも確率雨量がやや大きく見積もられるようですね．
計算用エクセルファイルは下記URLからダウンロードできます． ご自由にお使いください． 1水文統計.xlsm
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 Takase_Method()
'最終更新日：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 SK As Double 'ラグランジュ補間法によって得られたσξの値
Dim SK_N As Integer 'σξとNのデータセットの数(p.12の表1.3.1)
Dim L1 As Double 'ラグランジュ補間法の途中式
Dim l As Integer 'ラグランジュ補間法でfor文のカウント用に使う
Dim PK As Double 'ラグランジュ補間法の多項式P(k)
Dim K1 As Double '式(1.3.1)の第二項
Dim K2 As Double '式(1.3.1)の第一項
Dim K3 As Double 'log10(x) 式(1.3.1)の右辺
'*********************定数などの定義*********************
m = 56 'ξデータの数
X5 = 0 '自然対数をとったデータの合計値
x9 = 0 '自然対数をとったデータの偏差二乗和
p = 3.14159 '円周率
PL = 2.30259 'ln10(10の自然対数)
SK_N = 15 'σξとNのデータセットの数(p.12の表1.3.1)
PK = 0 'ラグランジュ補間法の多項式P(k)
'*********************データ数および計算条件の入力*********************
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)
'*****************************************************************************************
'*********************ラグランジュ補間法によりNに対するσξの値を得る*********************
'*****************************************************************************************
'Nが70以上のときは線形補間する(予備的計算の結果)
If N > 70 Then
k = 1
Do Until N < Cells(k + 2, 10)
k = k + 1
Loop
SK = Cells(k + 1, 11) + (Cells(k + 2, 11) - Cells(k + 1, 11)) / (Cells(k + 2, 10) - Cells(k + 1, 10)) * (N - Cells(k + 1, 10))
Else
For k = 1 To SK_N
'N<10のときは出力して終了
If N < 10 Then
MsgBox "Nの値が小さすぎます．"
Exit Sub
End If
L1 = 1
For l = 1 To SK_N
If l <> k Then
L1 = L1 * (N - Cells(l + 1, 10)) / (Cells(k + 1, 10) - Cells(l + 1, 10))
End If
Next l
PK = PK + Cells(k + 1, 11) * L1
Next k
SK = PK
End If
Debug.Print "standard Deviation"
Debug.Print "σ0="; Cells(N + 7, 2)
Debug.Print "σξ="; SK
Debug.Print ""
'*********************計算結果用ラベル*********************
Cells(1, 16) = "高瀬法(積率法)"
Cells(2, 16) = "T(year)"
Cells(2, 17) = "ξ"
Cells(2, 18) = "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 / SK '(σ0/σξ)*ξ 式(1.3.1)の第二項
K2 = Log(X0) / PL 'log10(x0) 式(1.3.1)の第一項
K3 = K1 + K2 'log10(x) 式(1.3.1)の右辺
x = Exp(PL * K3) '式(1.3.1)の左辺のx，つまり，T年確率降水量
Cells(3 + cnt, 16) = t
Cells(3 + cnt, 17) = E2
Cells(3 + cnt, 18) = x
cnt = cnt + 1
Next t
Debug.Print Format(Date) & " " & Hour(Time) & "時" & Minute(Time) & "分" & Second(Time) & "秒(S.K)"; " 計算終了"
Application.ScreenUpdating = True
End Sub
5．おわりに
防災や減災，人命救助などに関心のある人が農業農村工学に興味を持つきっかけとして本記事が貢献できれば幸いです．
【確率雨量】シリーズ一覧
【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第1回 対数正規法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第2回 高瀬法(積率法)― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第3回 岩井法― 【確率雨量】100年に1度の大雨ってどうやって計算しているの？ ―第4回 ガンベル・チャウの方法―
本記事の執筆にあたり，書籍の内容について紹介することを快く承諾いただいた森北出版の皆様と，何より書籍の著者である河村三郎先生には心より感謝いたします．
確率雨量シリーズの第1回でも紹介しましたが，本記事の執筆にあたり，以下の書籍を参考にしました．この場を借りて厚く御礼申し上げます． 本記事で公開するプログラムは河村先生の書籍に書かれているもの(Basic)を元に，筆者がVBAで書いたものです． ミスには細心の注意を払っておりますが，完璧を保証するものではありませんので，以下の書籍をお手元にご用意いただくことをお勧めいたします．
河村三郎：土木工学プログラム集1 水文・水理1 実用プログラムとグラフィック表示，森北出版，1984 岩井重久，石黒政儀：応用水文統計学，森北出版，1970 石村園子：やさしく学べるラプラス変換・フーリエ解析 増補版，共立出版，2010
今後も農業農村工学(水文学，かんがい排水，土壌物理，水理学)を中心に記事を執筆していきたいと思います． リクエスト等も受け付けておりますので，ご遠慮なく連絡ください． Twitterアカウント：エビぐんかん@6LxAi9GCOmRigUI メール：nnCreatorCircle@gmail.com
検索用Key Words： 確率雨量，確率降水量，降雨，降水，水文事象，水文学，水文統計学，降雨事象，水利施設，設計，確率水文量，水文諸量，超過確率，非超過確率，リターンピリオド(return period)，T年確率値，事象値，計算方法，解き方，プログラム，マクロ，やり方
