ノンパラメトリックな発想による確率降水量の標本分布のパラメトリックな近似

− 101 −

ノンパラメトリックな発想による確率降水量の標本分布のパラメトリックな近似 名古屋工業大学　北野 利一 １．はじめに 　河川整備の治水計画には，その流域に降る雨が河川に流出する基本高水流量を求める必要があり，これに より，対象となる河川の流域ならびに想定氾濫区域内の人口や資産などをもとに，氾濫した場合に想定され る最大被害額などを求め，その被害を防御もしくは減ずる対策を講じる1)．この時，計画の規模となるのが， 再現期間2)である．すなわち，再現期間となる R 年間に平均１回生じる降水量 µR を，確率降水量と定め るのである．総降水量に対応するような２日降水量や３日降水量を用いる場合もあるが，多くの場合は，日 降水量（あるいは，日界を変えて最大となる24 時間降水量）で代表させて，確率降水量（再現レベル）を 定めた後に，さまざまな降雨パターンを用いて，基本高水を決定している．自然外力の極値に関して，沿岸 域の防災に関わる高波・高潮潮位も再現期間を用いた類似した取扱いを行なっている3)． 　このような検討の際には，多くの不確実性が含まれることに留意する．下流域に大都市をかかえる主要な 河川に対しては，再現期間を100 年あるいは 200 年としており，これに対して，降雨の観測記録は限られる． 100 年を超える観測点もあるが，同一条件での観測の連続性などを考慮すると，観測記録から所与の再現期 間に対して推定される確率降水量には，無視できない推定誤差が含まれるため，従来から，その信頼区間に ついて議論されてきた．外挿となる確率降水量の信頼区間は，一般的に，右に歪む（点推定値から信頼区間 の上端までの長さが，下端までの長さよりも長い）と期待されるが，確率降水量の推定誤差分布を正規分布 で与えると，点推定値に対して対称な幅を与えてしまう．このような根本的な欠点を解消するために，プロ ファイル尤度を用いた信頼区間の算定法もある4, 5)．しかしながら，パラメトリックな方法は，統計学的な 推定理論に基づくため，当然，直感的に理解できるものではない．これは，降水や河川への流出の現象を扱 う気象学，水文学ならびに土木工学の立場からすると，統計理論が一種のブラックボックスに覆われたもの になり，不確実性の幅の大きさだけでなく，考え方そのものも不透明となる難点が生じている．そこで，本 研究では，直感的な思考につながりやすいノンパラメトリックな発想をもとに，パラメトリックな推定結果 を用いて，確率降水量の標本分布を検討するものである． ２．ノンパラメトリックな確率降水量の信頼区間 　本研究で例示に用いる対象データは，「地球温暖化対策に資するアンサンブル気候予測データベース （d4PDF）」 の成果6)を用いて，日本域 20km のダウンスケールによる出力のうち，名古屋を含む格子の過去 実験データ出力値（1951-2010 年の 60 年× 50 アンサンブル標本）を用いる．各アンサンブル標本は，独立 図 -1 POT による日降水量の順序統計量（3000 年分） 図 -2 生起数の変動（ポアソン分布） Count Number Probability 0 1 2 3 4 5 6 7 8 9 0. 1. 2 91.7 % Count Number Probability 0 5 10 15 20 25 30 0 .05 .1 90.9 % a) 平均 200 年に１回のイベントの生起数分布 b) 平均 1000 年に１回のイベントの生起数分布 50 100 200 500 1000 2000 250 300 350 400 450

Return Period [Year]

Precipitation [mm]

− 102 −

ノンパラメトリックな発想による確率降水量の標本分布のパラメトリックな近似 名古屋工業大学　北野 利一 １．はじめに 　河川整備の治水計画には，その流域に降る雨が河川に流出する基本高水流量を求める必要があり，これに より，対象となる河川の流域ならびに想定氾濫区域内の人口や資産などをもとに，氾濫した場合に想定され る最大被害額などを求め，その被害を防御もしくは減ずる対策を講じる1)．この時，計画の規模となるのが， 再現期間2)である．すなわち，再現期間となる R 年間に平均１回生じる降水量 µR を，確率降水量と定め るのである．総降水量に対応するような２日降水量や３日降水量を用いる場合もあるが，多くの場合は，日 降水量（あるいは，日界を変えて最大となる24 時間降水量）で代表させて，確率降水量（再現レベル）を 定めた後に，さまざまな降雨パターンを用いて，基本高水を決定している．自然外力の極値に関して，沿岸 域の防災に関わる高波・高潮潮位も再現期間を用いた類似した取扱いを行なっている3)． 　このような検討の際には，多くの不確実性が含まれることに留意する．下流域に大都市をかかえる主要な 河川に対しては，再現期間を100 年あるいは 200 年としており，これに対して，降雨の観測記録は限られる． 100 年を超える観測点もあるが，同一条件での観測の連続性などを考慮すると，観測記録から所与の再現期 間に対して推定される確率降水量には，無視できない推定誤差が含まれるため，従来から，その信頼区間に ついて議論されてきた．外挿となる確率降水量の信頼区間は，一般的に，右に歪む（点推定値から信頼区間 の上端までの長さが，下端までの長さよりも長い）と期待されるが，確率降水量の推定誤差分布を正規分布 で与えると，点推定値に対して対称な幅を与えてしまう．このような根本的な欠点を解消するために，プロ ファイル尤度を用いた信頼区間の算定法もある4, 5)．しかしながら，パラメトリックな方法は，統計学的な 推定理論に基づくため，当然，直感的に理解できるものではない．これは，降水や河川への流出の現象を扱 う気象学，水文学ならびに土木工学の立場からすると，統計理論が一種のブラックボックスに覆われたもの になり，不確実性の幅の大きさだけでなく，考え方そのものも不透明となる難点が生じている．そこで，本 研究では，直感的な思考につながりやすいノンパラメトリックな発想をもとに，パラメトリックな推定結果 を用いて，確率降水量の標本分布を検討するものである． ２．ノンパラメトリックな確率降水量の信頼区間 　本研究で例示に用いる対象データは，「地球温暖化対策に資するアンサンブル気候予測データベース （d4PDF）」 の成果6)を用いて，日本域 20km のダウンスケールによる出力のうち，名古屋を含む格子の過去 実験データ出力値（1951-2010 年の 60 年× 50 アンサンブル標本）を用いる．各アンサンブル標本は，独立 図 -1 POT による日降水量の順序統計量（3000 年分） 図 -2 生起数の変動（ポアソン分布） Count Number Probability 0 1 2 3 4 5 6 7 8 9 0. 1. 2 91.7 % Count Number Probability 0 5 10 15 20 25 30 0 .05 .1 90.9 % a) 平均 200 年に１回のイベントの生起数分布 b) 平均 1000 年に１回のイベントの生起数分布 50 100 200 500 1000 2000 250 300 350 400 450

Return Period [Year]

Precipitation [mm] 100 50 22 15 10Rank 6 3 1 と見なせるので，一括してまとめれば，3000 年分の日降水量の記録がある．この 3000 年分のうち，上位 100 番目までの日降水量（Peaks-over-threshold 法により抽出）に対して，順位 i の降水量に，その生起率を i /3000 で与え，横軸に再現期間（生起率の逆数）を対数でとり（補助軸として，上部の横軸には順位を付与）， 縦軸に降水量をプロットしたものを図 -1 に示す．一雨毎のピークとなる日降水量が大きな値をとる降雨イ ベントが，ポアソン過程で生起するとすれば，例えば，L = 200 年に平均１回生起するイベントが，3000 年間 の生起する回数 k に対する確率は，次式に示されるポアソン分布で与えられる（図 -2 a) を参照）． pL(k) = λkLe−λL/k! (1) したがって，生起率 λ1 = 1/200 yrs.−1(= λ3000/3000) のイベントは，確率 0.91 で，3000 年間に 10 回から 22 回の範囲で生起することから，200 年確率降水量の区間推定として，図 -1 の暖色の破線に示すとおり，上位 10 番目と22 番目の日降水量の順序統計量を区間の端点として求めることができる．データが 3000 年分ある強み から，生起率 1/1000 yrs.−1 のイベントも，図 -2 b) のポアソン分布に示すとおり，確率 0.92 で，3000 年間に 1 回から6 回の範囲で生起することから，1000 年確率降水量に対しても，同じようにノンパラメトリックな区間 推定が可能である（図 -1 での寒色の破線）．なお，この場合，3000 年間における期待生起数は 3 であり，この時， 生起数がゼロとなる確率がキリの良い数字 0.05 となることにちなんだ警句は，３の法則と知られる7, 8, 9)_． ３．レベル降水量を超過する生起率の推定 　50 アンサンブル標本の各々に対して，ある値の降水量（レベル降水量）を超える生起数を求め，これを 60 年間あたりの生起率の推定値とみなせば，生起率の標本値分布は，次式のガンマ分布で近似できる． fL(λ1) = LKλK−11 e−Lλ1/Γ(K) (2) 例えば，日降水量 200 mm を超える 60 年間の生起数の平均は 3.35 となり，K = 3.35, L = 60 yrs. を式 (2) に代 入して得られるガンマ分布は，図 -3a) に示すとおり，アンサンブル標本毎に得られる生起率のヒストグラムと 良好に適合する． 　期間最大値 y がしたがう漸近分布として知られる一般化極値分布（GEV）は，L 年間あたりの生起率を λL(y) =  1 + ξ y− µL σL −1/ξ (3) で与え，期間 L = 1 年間にレベル y を超える生起数をゼロ（k = 0）として，期間最大値 y の累積確率が導か れる．それにより得られる尤度に，各アンサンブル標本毎の年最大値データを代入して，母数 θ1 = {µ1, σ1, ξ} を推定し，式(3) を介することにより，レベル y = 200, 400 を超える生起率の推定値が標本サイズ50 で得ら れる．ここでは，母数推定には，無情報事前分布を用いたベイズ法（チューニングの必要のないRU 法10, 11)） で母数の乱数を生成し，式(3) により得られる生起率の平均を求めて点推定値 _λ˜_1(y) としている．このようにし て得られた生起率のヒストグラムを図 -3b) に示す．この場合はアンサンブル標本数がそれほどには大きくない 図 -3 アンサンブル標本毎に得られる生起率のヒストグラムとガンマ分布 　a) カウント数そのもの（200 mm）　 b) 生起率関数による生起率の推定（左：200 mm ／右：400 mm）

Occurrence rate [1/yr]

Density 0.00 0.05 0.10 0.15 0.20 05 10 15 u = 200 [mm] K = 6.57

Occurrence rate [1/yr]

Density 0.000 0.005 0.010 0.015 0.020 05 0 100 150 u = 400 [mm] K = 0.61

Occurrence rate [1/yr]

Density 0.00 0.05 0.10 0.15 0.20 05 10 15 u = 200 [mm] K = 3.35

− 103 −

が，汎用的な拡張を考える12)と，数値的な最適化を必要とする最尤法を避けた．また，図 -3b) で併せて示す

曲線は，次式で得られる K と L を用いて描かれるガンマ分布の密度関数である．図 -3a) のカウント数を期間 長で除して得られる生起率よりも，パラメトリックな生起率関数を介するため適合性は高まる．

K = {λ1(y)}2V{ˆλ1(y)}, L = Kλ1(y) (4ab) これらは，ガンマ分布の期待値と分散の関係に基づくものであるが，真値は不明であるので，各アンサンブル 標本の母数乱数を平均して得られる推定値の平均 θ1 を代入して得られる生起率を λ1(y) とし，分散 V{ˆλ1(y)} は，Fisher 情報行列を θ1 で評価したものを用いて，デルタ法により求めた． ４．確率降水量と上位 K 番めの順序統計量の密度関数 　再現期間 R 年とする確率降水量は，年最大値分布の位置母数に一致し，式 (3) の生起率関数の逆関数として， λ1(µR) = 1/R → µR = µ1 + σ1Rξ− 1/ξ (5) として与えられる．ここで，各アンサンブル標本の年最大値から生成される母数 θ1 のベイズ乱数を上式に代入 し，得られる µR の平均を確率降水量の点推定値とし，その標本分布を図 -4 に示す（R = 20 年の場合）．なお， 図中の太線は，その理論分布として，式(2) のガンマ分布を式 (3) の生起率関数を介して， dFL(µR) =− fL(λ1)dλ1   λ1=λ1(µR)= dµR σLΓ(K)  1 + ξµR− µL σL −K/ξ−1 exp  −  1 + ξµR− µL σL −1/ξ (6) で得られるものであり，これは，L 年間における上位 K 番めの最大値の確率分布4, 5, 13, 14) の密度関数に一致 する．ここで，式(4b) ならびに (5) より，L = KR となる．R = 20 年の場合には，この場合，K = 4.70 と求 められるので，L = 94.0 年となる．したがって，順位の数字 K を本来の整数から実数値に拡張したものと理 解すれば，再現期間 R = 20 年の確率降水量の標本分布は，L = 94.0 年間における上位 K = 4.70 番めの最大 値分布に相当するものと解釈できる．ノンパラメトリックに考えれば，3 (= 60/20) 番めの順序統計量について， 標本によるヒストグラムとその理論分布を図 -5 に示す．当然であるが，パラメトリックな推定により，標本分 布の集中の度合いが高くなる（誤差のバラツキは狭くなる）ことがわかる．K は，北野ら15)で導入された経 験度とよばれる量で，対象とする再現期間の確率降水量を推定するのに用いた情報の量をデータの個数で表現 したものである．形状母数 ξ と観測年数（この場合，60 年）に依存し，図 -6 に示すとおり，再現期間が長く なるにつれて，経験度 K は低下する．再現期間 R = 200 年は，観測期間 60 年を超えるため，もはやノンパラ 図 -4 確率降水量の標本分布（再現期間 20 年）

Return level [mm] for 20 years retun period

Density

150 200 250 300

0.00

0.01

0.02 _{K = 4.70}

3rd. largest in 60 yrs. precipitation [mm]

Density 150 200 250 300 0.00 0.01 0.02 図 -5 期間 60 年における上位３番目の順序統計量の分布

Return Period [years]

Degree of Exper ience .1 1 10 100 1000 60 10 2 1 .5 R = 20 60 200 yrs. K = 4.7 1.8 0.75

Return level [mm] for 20 years retun period

Density 200 400 600 800 0.000 0.003 0.006 K = 0.75 図 -6 再現期間に応じた経験度の値 図 -7 確率降水量の標本分布（再現期間 200 年）

− 104 −

が，汎用的な拡張を考える12)と，数値的な最適化を必要とする最尤法を避けた．また，図 -3b) で併せて示す

曲線は，次式で得られる K と L を用いて描かれるガンマ分布の密度関数である．図 -3a) のカウント数を期間 長で除して得られる生起率よりも，パラメトリックな生起率関数を介するため適合性は高まる．

K = {λ1(y)}2V{ˆλ1(y)}, L = Kλ1(y) (4ab) これらは，ガンマ分布の期待値と分散の関係に基づくものであるが，真値は不明であるので，各アンサンブル 標本の母数乱数を平均して得られる推定値の平均 θ1 を代入して得られる生起率を λ1(y) とし，分散 V{ˆλ1(y)} は，Fisher 情報行列を θ1 で評価したものを用いて，デルタ法により求めた． ４．確率降水量と上位 K 番めの順序統計量の密度関数 　再現期間 R 年とする確率降水量は，年最大値分布の位置母数に一致し，式 (3) の生起率関数の逆関数として， λ1(µR) = 1/R → µR = µ1 + σ1Rξ− 1/ξ (5) として与えられる．ここで，各アンサンブル標本の年最大値から生成される母数 θ1 のベイズ乱数を上式に代入 し，得られる µR の平均を確率降水量の点推定値とし，その標本分布を図 -4 に示す（R = 20 年の場合）．なお， 図中の太線は，その理論分布として，式(2) のガンマ分布を式 (3) の生起率関数を介して， dFL(µR) =− fL(λ1)dλ1   λ1=λ1(µR)= dµR σLΓ(K)  1 + ξµR− µL σL −K/ξ−1 exp  −  1 + ξµR− µL σL −1/ξ (6) で得られるものであり，これは，L 年間における上位 K 番めの最大値の確率分布4, 5, 13, 14) の密度関数に一致 する．ここで，式(4b) ならびに (5) より，L = KR となる．R = 20 年の場合には，この場合，K = 4.70 と求 められるので，L = 94.0 年となる．したがって，順位の数字 K を本来の整数から実数値に拡張したものと理 解すれば，再現期間 R = 20 年の確率降水量の標本分布は，L = 94.0 年間における上位 K = 4.70 番めの最大 値分布に相当するものと解釈できる．ノンパラメトリックに考えれば，3 (= 60/20) 番めの順序統計量について， 標本によるヒストグラムとその理論分布を図 -5 に示す．当然であるが，パラメトリックな推定により，標本分 布の集中の度合いが高くなる（誤差のバラツキは狭くなる）ことがわかる．K は，北野ら15)で導入された経 験度とよばれる量で，対象とする再現期間の確率降水量を推定するのに用いた情報の量をデータの個数で表現 したものである．形状母数 ξ と観測年数（この場合，60 年）に依存し，図 -6 に示すとおり，再現期間が長く なるにつれて，経験度 K は低下する．再現期間 R = 200 年は，観測期間 60 年を超えるため，もはやノンパラ 図 -4 確率降水量の標本分布（再現期間 20 年）

Return level [mm] for 20 years retun period

Density

150 200 250 300

0.00

0.01

0.02 _{K = 4.70}

3rd. largest in 60 yrs. precipitation [mm]

Density 150 200 250 300 0.00 0.01 0.02 図 -5 期間 60 年における上位３番目の順序統計量の分布

Return Period [years]

Degree of Exper ience .1 1 10 100 1000 60 10 2 1 .5 R = 20 60 200 yrs. K = 4.7 1.8 0.75

Return level [mm] for 20 years retun period

Density 200 400 600 800 0.000 0.003 0.006 K = 0.75 図 -6 再現期間に応じた経験度の値 図 -7 確率降水量の標本分布（再現期間 200 年） メトリックには推定できないが，K = 0.75 となり，パラメトリックには推定可能であり，図 -7 に示すとおり， 標本毎に推定される確率降水量のヒストグラムと，経験度を用いた上位 K 番めの最大値分布は，良好に一致 している．もう１つ興味深い例を付け加えると，図 -8 に示す再現期間 60 年の確率降水量の標本分布（左）と， 期間 60 年の最大値分布（右）の対比である．いずれも，アンサンブル標本によるヒストグラムと，その理論曲 線は良好な一致を確認できるだけでなく，その理論曲線は共通の数式で表されるものであり，一方は，標本か ら直接得られるもの（60 年最大値）という確率変数の分布であり，もう一方は，標本からパラメトリックな推 定で得られる（60 年最大値分布の位置母数の）推定値の標本分布である．ここで，60 年最大値分布の代表値 を位置母数とみなせば，図 -8 右は60 年最大値のみを用いて推定した場合に相当し，他方，図 -8 左は 60 年間 の年最大値 60 個を用いて推定しているので，後者の情報が大きく（K = 1.8 > 1），それゆえに，そのバラツキ の幅が小さいことが確認できる．また，これが，極値統計解析により外挿できる量的な根拠とも言える． ５．まとめ 　土木工学などの応用分野で，観測記録を用いて極値統計解析による外挿する場合に，どの程度の推定誤差が 含まれるのかを示すため，信頼区間を用いて提示してきた．しかしながら，信頼区間の端点を実際に計算しな いと見当もつなかいようでは，ブラックボックスと同じである．本研究で示すとおり，再現期間 R 年の確率降 水量の標本分布が，L (= KR) 年間における上位 K 番めの極値分布で近似できることを実例により確認した． なお，観測データの情報の量に応じて，経験度 K の値が与えられる．このような仕組みは，極値統計解析によ る推定結果を理解し，また，その不確実性を理解した上で推定結果を応用できるようになる点で有用である． 参考文献 1) 高橋　裕 : 河川工学，東京大学出版会，311p., 1990. 2) 池淵周一 ・椎葉充晴 ・宝　 馨 ・立川康人 : エース水文学，朝倉書店，216p., 2006. 3) 合田良実 : 耐波工学 − 港湾・海岸構造物の耐波設計，鹿島出版会，430p., 2008. 4) Coles, S.: An Introduction to statistical modeling of extreme values, Springer, 208p., 2001.

5) 高橋倫也・志村隆彰 : 極値統計学，ISM シリーズ：進化する数理統計 5，近代科学社，262p., 2016.

6) 地球温暖化対策に資するアンサンブル気候予測データベース（database for Policy Decision making for Future climate change, d4PDF），http://www.miroc-gcm.jp/~pub/d4PDF/index.html, 2016.

7) Hanley, J. A., and A. Lippman-Hand (1983): If nothing goes wrong, is everything all right? - Interpreting zero numerators, Jour. Amer. Medical Assoc., Vol. 249(13), pp.1743-5.

8) Jovanovic, B. D. and P. S. Levy (1997): A look at the rule of three. The American Statistician, Vol.51, pp.137-139. 9) 北野利一：伊勢湾台風級の高潮と確率潮位，第 49 回 水工学に関する夏期研修会講義集，B-3-1 ~ 20, 2013.

10) Wakefield, J.C., A. E. Gelfand & A. F. M. Smith: Efficient Generation of Random Variates via the Ratio-of-Uniforms Method, Statistics and Computing, Vol.1, pp.129–133, 1991.

11) Northrop, P. J.: rust: Ratio-of-Uniforms Simulation with Transformation. https://CRAN.R-project.org/package=rust, 2018. 12) 北野利一・上野玄太・森 聡紫・森 亮太 : 上位 K 番め最大値分布を用いた確率波高の標本分布の近似，土木学会論

文集B2（海岸工学），第 75 巻，（投稿中）

13) Weissman, I.: Estimation of parameters and quantile based on the k largest observations, Jour. of Amer. Stat. Assoc., Vol.73, pp.317-327, 1978.

14) Tawn, J. A.: An extreme value theory model for dependent observations, Jour. of Hydrology, Vol.101, pp.227-250.

15) 北野利一・森瀬喬士・喜岡　渉・高橋倫也 : 確率波高に対する推定の可否を決定づける新たな指標の提案，海岸工学 論文集，第55 巻，pp.141-145, 2008.

図 -8 再現期間 60 年の確率降水量の標本分布（右）と期間 60 年における最大値の分布（左）

Return level [mm] for return period 60 years

Density 200 250 300 350 400 450 0.000 0.006 0.012 _{K = 1.8} 60−year−max. precipitation [mm] Density 200 250 300 350 400 450 0.000 0.006 0.012