極値解析が用いられる事例紹介(2)
豪雨災害（河川工学）

京都大学 防災研究所 水資源環境研究センター
田中茂信

数学協働プログラム
甚大災害の外力想定に必要となる極値統計解析法の背景と活用
20141208

本日の内容

1)治水計画の歴史と確率を取入れた河川計画とその背景

2)水文時系列資料の極値解析

Mann - Kendall検定によるトレンドの評価
	AMS資料の確率評価
	POT資料の確率評価
	AMSとPOTの比較

治水計画の歴史と確率を取り入れた河川計画とその背景

明治初期  欧米の技術がオランダ人技術者をはじめとする欧米の技術者の来日
	国　：港湾の整備と一体となった河口部処理及び舟運を主目的とする低水工事
	地方：高水工事（氾濫が一地域に限る地先的なものとの認識）
M18年，20年代の大水害をうけ明治29年旧河川法制定，直轄工事が始まる
M43年(1910)　全国的な大水害をうけ，勅令を以て内閣に臨時治水調査会設置，第一期治水計画検討
T10 年(1921)　第二期治水計画策定
1927 M. Fréshet，1928 R.A. Fisher and L.H.C. Tippett，1943 Gnedenko　極値理論の研究
	極値分布には3つの型があって，かつ，それ以外にない．
S8年(1933)　　臨時治水調査会に代わり土木会議により第三期治水計画策定
S9年(1934)　　水害防止協議会
氾濫区域内の道路,鉄道などの路線のうち重要なものの高さは大体20年に1回起こる程度の洪水の水位(最 大洪水位)より30cm以上を標準とするとしている
富永は計画の洪水流量としては,既往の最高流量を原則とするが,周囲の状況または財政の関係より10年に1回起こる程度の洪水流量で我慢することもあるとしている (山本,1999)
改修途中でその計画規模を超える洪水が発生した利根川等では,その度に計画の見直し
1945年，47年，48年と1910年に匹敵する大災害が相次いで起こり,その対応として従来の河川改修計画を拡大した他に,できるだけ下流河川の拡幅を避けるため,貯水池,遊水地などが上流部に計画されている.また,多目的ダムの計画・着工が行われた
S25年(1950)　中安米蔵「治水計画における洪水流量についてー千代川を中心としてー」発表　防災利益率の概念導入　岩井（1947）を使用
S27年(1952)　米田正文「淀川計畫高水論」発表　確率主義は,治水計画を理論的にし,安全かつ信頼性のある計画と主張するが,
		既往最大洪水にも重きを置いている　確率紙に平分線を引く方法を使用
1955 Jenkinson 　一般極値分布(generalized extreme- value distribution)
1956 角屋　睦「極値分布とその一解法」　対数極値分布A型および対数極値分布B型
1958 E.Gumbel “Statistics of Extremes”
1958「建設省河川砂防技術基準」制定
「確率計算の方法としてはFuller, Hazen, Foster, Goodrich, Slade, Kimball, 岩井,Gumbel等の方法,非母数法と漸次改良され新しい提案がなされてきた.しかしいず れの方法とも全国一律の基準と決定するに至らないが,現在のところ岩井法を統計の要素の下限値が 正又は零,上限値が無限大になるという条件のもとに採用するか,Gumbel法によるか,あるいは確 率紙にHazenのプロッティングポジションを用いてプロットし,平分線を引く簡略法を用いる場合が 多い.ただし,簡略法は取り扱いに便利であるが,統計的健全性に欠け,客観性に乏しい欠点があ る」
1958山本三郎「河川工学」
Fullerの方法,正規分布,対数正規分布,岩井法,Hazen紙(対数確率紙)による方法およびGumbelの方法 が16頁にわたって説明されている.　なお,Gumbelの方法については極値理論に基づくものと説明
1961 角屋　睦　年最大値(AMS; Annual Maximum Series)資料と非毎年値(POT)資料の理論的関係や応用例を示した
1963 和訳「極値統計学」出版
1964 角屋　睦「水文統計論」
1970 岩井重久，石黒政儀「応用水文統計学」
1974,1975 Balkema and de Haan(1974)とPikands(1975)　一般パレート分布の導出
1979 GreenwoodらPWM(確率加重積率)法を提案
1988 宝・高棹　分布モデルの評価基準として適合度のみでなく確率水文量の安定性も考慮すべき
	SLSC（Standardized Least Square Criterion)
1990 宝・高棹　3母数対数正規分布，Gumbel分布一般極値分布の母数推定法：小標本の場合L積率法がよい

Mann-Kendall検定によるトレンドの評価

a special case of Kendall’ s Correlation Coefficient with independent variable of t(time)

Mann - Kendall検定はトレンドが線形か非線形かを問わずに時系列資料のトレンドを検定するノンパラメトリックな方法である．
本検定の帰無仮説(Subscript[H, 0])と対立仮説(Subscript[H, 1])は以下に示すとおりである．
	(Subscript[H, 0]) : n個のデータ {(Subscript[x, 1]), (Subscript[x, 2]), ... , (Subscript[x, n])} が独立で同一の確率分布にしたがう．
	(Subscript[H, 1]) : n個のデータ{(Subscript[x, 1]), (Subscript[x, 2]), ... , (Subscript[x, n])} が同一の確率分布にしたがわない．
Mann - Kendall検定においては，次式で与えられる統計量Zを定義する．

S=(∑_(k=1)^(n-1))(∑_(j=k+1)^n)sign((Subscript[x, j])-(Subscript[x, k]))                                                            (1)

sign(θ)=({{\[Piecewise], ({{1, θ>0}, {0, θ=0}, {-1, θ<0}})}})                                                                    (2)

S は漸近的に正規分布し，平均がゼロで分散が式(3)で表される．

Var(S)=((n(n-1)(2n+5)-(∑_(i=1)^m)(Subscript[e, i])((Subscript[e, i])-1)(2(Subscript[e, i])+5))/18)       (3)

Z=({{\[Piecewise], ({{((S-1)/(Sqrt[Var(S)])), S>0}, {0, S=0}, {((S+1)/(Sqrt[Var(S)])), S<0}})}})                                                                  (4)

ここに，(Subscript[e, i])はデータを昇順に並べたとき，同じ値が連続して出現する個数を表し，mはその組数を表す．式 (4) で表されるZは n が小さい場合でも近似的に標準正規変量であることが示されている．

RunLengthEncode::usage="RunLengthEncode[dataseries]は,一連のリストの中で連続して同じ値がいくつ続くかを順に{値,続く個数}のリストで返す.";
RunLengthEncode[x_List]:={First[#],Length[#]}&/@Split[x];
MannKendall::usage="MannKendall[dataseries]はZ値を返す.MannKendall[dataseries,True]はZ値,S,S/D,分散を返す.\[LeftBracketingBar]Z\[RightBracketingBar]<1.96\!\(\*SubscriptBox[\(のときは帰無仮説H\), \(0\)]\)は棄却されない.
帰無仮説\!\(TraditionalForm\`\*SubscriptBox[\(H\), \(0\)]\\  : \\ n個のデータ\\ {\*SubscriptBox[\(x\), \(1\)], \\ \*SubscriptBox[\(x\), \(2\)], \\  ... \\ , \\ \*SubscriptBox[\(x\), \(n\)]}\\ が独立で同一の確率分布にしたがう．\*SubscriptBox[\(対立仮説H\), \(1\)]\\  : \\ n個のデータ {\*SubscriptBox[\(x\), \(1\)], \\ \*SubscriptBox[\(x\), \(2\)], \\  ... \\ , \\ \*SubscriptBox[\(x\), \(n\)]}\\ が同一の確率分布にしたがわない．\)";
MannKendall[x_List,ansfull_:False]:=Module[{n=Length[x],SRunLen,Vs2,S,Z2,D,tt2,ttti},SRunLen=SortBy[RunLengthEncode[Sort[x]],Last][[All,2]];ttti=#(#-1)(2#+5)&/@SRunLen//Total;tt2=#(#-1)&/@SRunLen//Total;
Vs2=(1/18)(n(n-1)(2n+5)-ttti);S=(∑_(k=1)^(n-1))(∑_(j=k+1)^n)Sign[x[[j]]-x[[k]]];Z2=((S-Sign[S])/(Sqrt[Vs2]));D=(Sqrt[(1/2)n(n-1)-(1/2)tt2])(Sqrt[(1/2)n(n-1)]);If[ansfull,{Z2,S,(S/D),Vs2}//N,N[Z2]]];

トレンド評価の事例

気象官署の資料

気象官署の年降水量や年最大値は定常かどうかを調べる 
資料は気象庁のWebsiteより取り込む

高知地点資料=Import["http://www.data.jma.go.jp/obd/stats/etrn/view/annually_s.php?prec_no=74&block_no=47893&year=&month=&day=&view=","Data"];

年を括弧でくくる[data_]:={{data[[1]]},data[[2]]};

Needs["PlotLegends`"]

年を括弧でくくる/@高知地点資料[[2,2,2,5;;,{1,5}]]

DateListPlot[年を括弧でくくる/@高知地点資料[[2,2,2,5;;,{1,5}]],PlotRange->All,FrameLabel->{"year","日降水量 [mm]"},PlotLabel->"高知 年最大日降水量"]

MannKendall[高知地点資料[[2,2,2,5;;-2,5]]]

これは有意水準5 % で帰無仮説が棄却されないことになる.すなわち, 年最大日雨量は定常な変動の範囲である.

ついでに，平均気温についてMannKendall検定を行う

MannKendall/@{高知地点資料[[2,2,2,5;;-2,8]]}

MannKendall[高知地点資料[[2,2,2,5;;-2,#]]]&/@{8,9,10}

これは有意水準5 % で帰無仮説が棄却されることになる.すなわち, 気温は有意に増加していると解釈できる.

DateListPlot[年を括弧でくくる/@高知地点資料[[2,2,2,5;;-2,{1,8}]],GridLines->{Automatic,{15,16,17,18}}]

DateListPlot[Map[年を括弧でくくる,{高知地点資料[[2,2,2,5+75;;-2,{1,11}]],高知地点資料[[2,2,2,5;;-2,{1,9}]],高知地点資料[[2,2,2,5;;-2,{1,8}]],高知地点資料[[2,2,2,5;;-2,{1,10}]]},{2}],GridLines->{Automatic,Automatic},PlotRange->{10,40},PlotStyle->{Red,Orange,Black,Blue},PlotLegend->{"年最高気温","日最高の年平均","日平均の年平均","日最低の年平均"},LegendPosition->{1.1,-0.4},LegendSize->{0.8,0.9},Joined->{True,True,True}]

極値分布による頻度解析

極値

連続資料や継続的に観測された資料のうち，局所的な最大値または最小値のことをそれぞれ極大値，極小値といい、これらを併せて極値とよぶ．一定区間，または一定期間の極大値を区分最大値という．期間が一年間の最大値は年最大値(Annual Maximum Series, AMS)と呼ばれ，区分最大値の一つである．AMSは取り扱いが容易なことからよく使われる．
年最大日平均気温や年最大時間降水量は年最大値であるが，年降水量は極値ではない．

日吉dam流域平均2日雨量の時系列とAMS

AMS日吉2day

AMS日吉2day={{{1952,7,10,0,0,0.`},170.3`},{{1953,9,24,0,0,0.`},245.5`},{{1954,9,17,0,0,0.`},116.3`},{{1955,8,21,0,0,0.`},80.80000000000001`},{{1956,9,25,0,0,0.`},105.6`},{{1957,6,26,0,0,0.`},127.5`},{{1958,8,24,0,0,0.`},88.9`},{{1959,8,12,0,0,0.`},305.8`},{{1960,8,29,0,0,0.`},327.5`},{{1961,10,26,0,0,0.`},224.60000000000002`},{{1962,6,9,0,0,0.`},97.19999999999999`},{{1963,6,13,0,0,0.`},62.4`},{{1964,9,24,0,0,0.`},100.69999999999999`},{{1965,9,16,0,0,0.`},279.`},{{1966,9,17,0,0,0.`},129.2`},{{1967,7,8,0,0,0.`},155.`},{{1968,8,28,0,0,0.`},72.7`},{{1969,6,25,0,0,0.`},125.1`},{{1970,6,14,0,0,0.`},158.3`},{{1971,8,30,0,0,0.`},177.6`},{{1972,9,15,0,0,0.`},202.5`},{{1973,10,12,0,0,0.`},70.9`},{{1974,9,8,0,0,0.`},88.1`},{{1975,8,22,0,0,0.`},160.4`},{{1976,9,8,0,0,0.`},127.60000000000001`},{{1977,9,28,0,0,0.`},92.6`},{{1978,6,21,0,0,0.`},95.4`},{{1979,6,28,0,0,0.`},140.1`},{{1980,11,21,0,0,0.`},69.5`},{{1981,6,25,0,0,0.`},72.4`},{{1982,7,31,0,0,0.`},179.3`},{{1983,9,27,0,0,0.`},204.9`},{{1984,6,25,0,0,0.`},71.`},{{1985,6,24,0,0,0.`},108.9`},{{1986,7,9,0,0,0.`},76.8`},{{1987,7,18,0,0,0.`},110.8`},{{1988,7,14,0,0,0.`},144.60000000000002`},{{1989,9,2,0,0,0.`},203.5`},{{1990,9,18,0,0,0.`},187.6`},{{1991,9,30,0,0,0.`},80.69999999999999`},{{1992,8,18,0,0,0.`},175.9`},{{1993,6,29,0,0,0.`},107.9`},{{1994,9,15,0,0,0.`},97.9`},{{1995,5,11,0,0,0.`},218.5`},{{1996,8,27,0,0,0.`},223.8`},{{1997,7,26,0,0,0.`},149.9`},{{1998,10,16,0,0,0.`},169.8`},{{1999,6,29,0,0,0.`},106.89999999999999`},{{2000,9,11,0,0,0.`},138.3`},{{2001,8,21,0,0,0.`},133.79999999999998`},{{2002,9,6,0,0,0.`},63.5`},{{2003,8,8,0,0,0.`},93.10000000000001`},{{2004,10,19,0,0,0.`},217.7`},{{2005,7,3,0,0,0.`},114.4`},{{2006,7,17,0,0,0.`},193.7`},{{2007,7,12,0,0,0.`},76.5`},{{2008,6,19,0,0,0.`},85.10000000000001`},{{2009,10,7,0,0,0.`},92.9769`},{{2010,7,13,0,0,0.`},169.1573`},{{2011,9,20,0,0,0.`},189.62099999999998`},{{2012,6,16,0,0,0.`},94.4328`},{{2013,9,15,0,0,0.`},340.9286`}}

年毎POTcr62

年毎POTcr62={{170.3`,141.2`,72.5`,130.1`,62.`,67.6`},{245.5`,120.8`,104.30000000000001`,65.2`,87.7`,86.8`,69.2`},{116.3`,108.30000000000001`,95.2`,70.1`,69.`,68.7`},{80.80000000000001`,68.`,74.`},{105.6`,88.3`},{127.5`,73.6`,104.9`,66.5`},{88.9`,75.7`},{305.8`,81.1`,166.2`,66.1`},{327.5`,117.30000000000001`,74.60000000000001`,67.6`,62.2`},{224.60000000000002`,111.5`,111.`,97.8`,87.`,71.`},{97.19999999999999`,78.8`,62.6`},{62.4`,62.4`},{100.69999999999999`},{279.`,132.1`,92.8`,129.4`,121.2`,129.`,76.69999999999999`,103.5`,80.39999999999999`},{129.2`,84.5`,64.5`,66.`,65.7`},{155.`,124.1`},{72.7`,70.9`},{125.1`,64.7`,114.3`,102.6`,84.3`,89.5`,70.2`},{158.3`,94.80000000000001`,78.`,70.`,65.7`},{177.6`,70.2`,130.70000000000002`},{202.5`,166.70000000000002`,95.9`,92.6`,72.8`,73.6`},{70.9`,62.4`},{88.1`,88.`,76.1`,66.5`,62.5`},{160.4`,74.6`,87.9`},{127.60000000000001`,71.9`,75.4`},{92.6`,63.5`,90.2`},{95.4`,75.69999999999999`,82.6`,74.6`},{140.1`,86.5`,82.`},{69.5`,65.6`,64.`},{72.4`,62.9`},{179.3`,79.7`},{204.9`,131.9`,116.69999999999999`,84.2`,87.9`},{71.`,63.4`},{108.9`,80.4`,88.7`,75.9`,76.2`,62.800000000000004`},{76.8`,71.4`,70.5`,70.9`,62.8`},{110.8`,92.69999999999999`,70.4`},{144.60000000000002`,132.6`,102.`,79.7`},{203.5`,72.7`,86.2`,65.6`,76.`,75.`},{187.6`,66.1`,69.9`},{80.69999999999999`,79.8`},{175.9`,82.3`},{107.9`,63.7`,98.5`,72.2`,64.7`,68.39999999999999`,62.6`},{97.9`,66.`,81.`},{218.5`,108.8`,78.7`,68.19999999999999`},{223.8`,76.5`},{149.9`,112.30000000000001`,130.9`,102.6`,74.4`,74.`,62.599999999999994`},{169.8`,146.20000000000002`},{106.89999999999999`,69.7`,94.5`,88.30000000000001`,77.3`},{138.3`,99.5`},{133.79999999999998`,100.39999999999999`},{63.5`},{93.10000000000001`,74.9`},{217.7`,109.`,65.1`,98.3`,75.5`,73.9`},{114.4`,65.3`},{193.7`},{76.5`},{85.10000000000001`,66.2`},{92.9769`,70.52590000000001`,69.9134`},{169.1573`,79.1596`,69.2817`,67.8896`},{189.62099999999998`,183.0625`,75.2091`,173.21409999999997`,177.7353`,156.96669999999997`},{94.4328`,91.65360000000003`,78.3237`},{340.9286`,90.73899999999999`,119.286`,80.2496`,62.960499999999996`}}

AMS高知

極値分布

極値水文資料を対象にする水文頻度解析において，多くの確率分布モデルが用いられてきた経緯がある．複数のモデルについて具体的な解説が行われている資料として，初期には角屋(1964)による「水文統計論」，岩井・石黒(1970)の「応用水文統計学」があげられる．近年では，水文・水資源ハンドブック(1997)の「7.3 水文頻度解析」がかなり詳しく紹介している．海外では，Handbook of Hydrology (1993）の第18章にStedingerらによる ”Frequency Analysis of Extreme Events” があり，水文極値解析に関する幅広い解説が述べられている．これらの文献で扱われている分布モデルと母数推定法について次表に示す．

一方，極値理論の研究者の間では，極値の解析には極値理論を用いることが当然のごとく王道と考えられている．計画外力の再現期間は既存の観測の蓄積期間に比べて大きい場合が一般的であるので，外挿のための理論的根拠を有する極値分布を極値解析に用いることは至極当然である．

【年最大値資料 (AMS) に用いる確率分布】

極値資料は3つの型の極値分布で表されることが証明されており，これらを一つの式で表したものが一般極値分布 (GEV : generalized extreme value) である．一般極値分布の形状母数が 0 の場合が Gumbel 分布であり，xを変量とするとき確率密度関数 f (x)，確率分布関数F(x) はそれぞれ次のように表される.

a) Gumbel 分布

f (x)=(1/α)exp {-((x-ξ)/α)}exp [-exp {-((x-ξ)/α)}]                            (5)

F(x)=exp [-exp {-((x-ξ)/α)}]                                                                (6)

変数xの分布範囲 :   -∞<x<∞

b) 一般極値分布 (GEV)

f (x)=(1/α)(({1-((k (x-ξ))/α)})^((1/k)-1))exp[-(({1-((k (x-ξ))/α)})^(1/k))]          (7)

F(x)=exp(-(((1-((k (x-ξ))/α)))^(1/k)))                                                         (8)

変数xの分布範囲 :  (Piecewise[{{-∞<x≤ξ+α/k ,    (k>0)}, {-∞<x<∞  ,    (k=0)}, { ξ+α/k≤x<∞,    (k<0)}}])   Gumbel分布

ここに ξ : 位置母数，α : 尺度母数，k : 形状母数である．

【閾値超過資料 (POT) に用いる確率分布】

水文時系列資料の頻度解析の方法として，AMS以外に閾値超過資料 (POT : peaks over threshold) がある．POTは，閾値を超過する独立なピーク値を全て取り出した資料である．AMSが他の年の年最大値より大きな年間第 2 位や 3 位などのデータを使わず，また，年によっては豪雨や洪水とみなせないような事象を含むのに対し，適切に閾値を選べばこのように特性の異なる水文資料を含むことを回避できる利点がある．閾値の選定によく用いられる手法としては，標本平均超過関数を用いる手法がある．POTに対して用いられる指数分布および一般Pareto分布 (GP : generalized Pareto distribution) の確率密度関数 f (x)，確率分布関数F (x) はそれぞれ次のように表される．

c) 指数分布 (Exp)

f (x)=(1/α)exp {-((x-ξ)/α)}                                                                         (11)

F(x)=1-exp {-((x-ξ)/α)}                                                                       (12)

変数xの分布範囲 :  ξ<x<∞

d) 一般パレート分布 (GPまたはGPD)

f (x)=(1/α) (({1-k((x-ξ)/α)})^(1/k-1))                                                               (13)

F(x)=1-(((1-k((x-ξ)/α)))^(1/k))                                                                  (14)

変数xの分布範囲:ξ<x≤ξ+α/k      (k > 0

ここに ξ : 位置母数（閾値），α : 尺度母数，k : 形状母数であり，k＝0の場合が指数分布である．

母数推定法(1)　積率法

表1より，古くは積率法が主に用いられていたが，「水文頻度解析」およびHandbook of HydrologyにおいてPWM法およびL積率法が加わっていることがわかる．この背景にはその当時の計算機資源などの環境が影響している．PWM(Probability Weighted Moments)法はGreenwoodら(1979)により提案された方法であり，L積率法(L-Moments)はそれを書き換えたものである．また，Regional Frequency Analysis(Hosking and Wallis, 1997)はL積率に関する詳しい解説と各分布のL積率による母数推定法を示している．実務では，積率法がアウトライヤー（はずれ値）の影響を受けやすいという欠点を改善するロバストな方法としてPWM(probability weighted moment:確率加重モーメント)法やL積率法が最近注目されている．このような状況に対し，極値理論の専門家の間では最尤法（さいゆうほう，Maximum likelihood estimation：MLE）がもっぱら用いられている．（「尤」は「もっともらしい」，likelihoodはlikely+hoodで起こる可能性の高いまたはもっともらしい＋状態という意味である．）しかし，最尤法は簡単な式で表せないことが多く，その場合，数値的に解かねばならないが，多くのソフトで実装されている．

【積率法】

確率分布モデルの原点及び平均値まわりの積率，すなわち，平均，分散，歪みをそれぞれ標本から得られる標本平均，不偏分散，不偏歪みに等しいとおいて分布モデルの母数を求める方法である．確率分布モデルの積率は母数の関数で表されているので，2母数の確率分布モデルは平均と分散に関する方程式を連立させて解く．3母数の場合はさらに歪みが不偏歪みに等しいとおく式を加えて連立方程式を解く．なお，3 母数の場合，種々の歪みの補正が提案され ている． 
確率分布モデルの確率密度関数をf(x)とするとき，この平均，分散及び歪みはそれぞれ式(15)，(16)および(17)で表される．

μ=(∫_(-∞)^∞)x f(x) \[DifferentialD]x                                                                                              (15)

(σ^2)=(∫_(-∞)^∞)(((x-μ))^2)f(x) \[DifferentialD]x                                                                             (16)

γ=(1/(σ^3))(∫_(-∞)^∞)(((x-μ))^3)f(x) \[DifferentialD]x                                                                        (17)

一方，標本 {(Subscript[X, 1]), (Subscript[X, 2]), ⋯, (Subscript[X, n])} から得られる標本平均，不偏分散，不偏歪みはそれぞれ式 (18)，(19) および (20) で表される．

(Subscript[(OverHat[μ]), X])=(OverBar[X])=(((∑_(i=1)^n)(Subscript[X, i]))/n)                                                                                             (18)

(Subscript[(OverHat[σ]), X]^2)=(S^2)=(((∑_(i=1)^n)((((Subscript[X, i])-(OverBar[X])))^2))/(n-1))                                                                            (19)

(Subscript[(OverHat[γ]), X])=((n(∑_(i=1)^n)((((Subscript[X, i])-(OverBar[X])))^3))/((n-1)(n-2)(S^3)))                                                                                 (20)

事例（Gumbel分布）

式 (15)の確率密度関数f(x)に式(5)を代入して積分する．

(∂_x) (E^(-(E^(-((x-ξ)/α)))))

fgum[x]:=(1/α)Exp [-((x-ξ)/α)]Exp [-Exp[-((x-ξ)/α)]]

model平均=Integrate[x fgum[x],{x,-∞,∞},Assumptions->{{x,ξ}∈Reals,α>0}](*21*)

EulerGammaは普通 γ で表され，約0.5772である．

EulerGamma//N

model分散=Integrate[(((x-model平均))^2) fgum[x],{x,-∞,∞},Assumptions->{{x,ξ}∈Reals,α>0}](*22*)

標本として日吉ダムの年最大2日雨量を用いる．式(18)と式(19)は

標本平均=Mean[AMS日吉2day[[All,2]]]

標本分散=Variance[AMS日吉2day[[All,2]]]

モデルの統計量と標本の統計量が1次と2次でそれぞれ等しいとおく連立方程式から, 母数 {α, ξ} を求める.

prmans=Solve[model平均==標本平均&&model分散==標本分散,{α,ξ}]

αが正のprmansの2番目の答えを使う．標本のヒストグラムと得られた母数で確率密度関数を描き比較する. このため，Histogramコマンドはオプション”PDF”を使用する．標本とモデルをShowコマンドで重ねて比較する.この時,最初の要素のフォーマットが使われる．描画範囲を指定し直すことも可能．

ヒストグラムと積率法の確率密度分布=Show[Histogram[AMS日吉2day[[All,2]],{40},"PDF"],Plot[Evaluate[fgum[x]/.prmans[[2]]],{x,0,400},PlotRange->All],PlotRange->{{0,400},All},Frame->True]

母数推定法(2)　L積率法

【L積率法】

大きく外れたデータが含まれ歪んだ水文資料に対処するために考案された手法がPWM(Probability-Weighted Moment, 確率加重積率)やL積率を用いる手法である．L積率は順序統計量の線形和で表される特徴を持つ(L は linear combinations に由来する)．
PWM は式(23)で定義される．通常の積率が変量xについて考えるのに対し，PWMでは確率Fあるいはuで考える．次数 r が大きくなっても常に変数の1乗しか含んでいないので，通常の積率に比べてアウトライヤーの影響を受けにくい．

(Subscript[β, r])=(∫_0^1)x(u)(u^r)\[DifferentialD]u                                                                                            (23)

(Subscript[X, j])(j=1, 2, ⋯, n)を標本から得られた順序統計量((Subscript[X, 1]) ≥ (Subscript[X, 2]) ≥ (Subscript[X, 3]), ⋯ , ≥ (Subscript[X, n]))とするとき，標本を用いたPWM推定量は式(24)~式(27)で表される．

(OverHat[(Subscript[β, 0])])=(OverBar[X])                                                                                                                    (24)

(OverHat[(Subscript[β, 1])])=(∑_(j=1)^(n-1))(((n-j)(Subscript[X, j]))/(n(n-1)))                                                                                            (25)

(OverHat[(Subscript[β, 2])])=(∑_(j=1)^(n-2))(((n-j)(n-j-1)(Subscript[X, j]))/(n(n-1)(n-2)))                                                                     (26)

(OverHat[(Subscript[β, 3])])=(∑_(j=1)^(n-3))(((n-j)(n-j-1)(n-j-2)(Subscript[X, j]))/(n(n-1)(n-2)(n-3)))                                              (27)

L積率は式 (28) 〜式 (31) のようにPWMを用いて表される．

(Subscript[λ, 1])=(Subscript[β, 0])                                                                                                                    (28)

(Subscript[λ, 2])=2(Subscript[β, 1])-(Subscript[β, 0])                                                                                                     (29)

(Subscript[λ, 3])=6(Subscript[β, 2])-6(Subscript[β, 1])+(Subscript[β, 0])                                                                                      (30)

(Subscript[λ, 4])=20(Subscript[β, 3])-30(Subscript[β, 2])+12(Subscript[β, 1])+(Subscript[β, 0])                                                                (31)

式 (32) および式 (33) で表されるL積率比はそれぞれ歪みおよび尖鋭度を表すものでモデルの当てはまり具合を判断する材料にもなる．

(Subscript[τ, 3])=(Subscript[λ, 3])/(Subscript[λ, 2])                                                                                                           (32)

(Subscript[τ, 4])=(Subscript[λ, 4])/(Subscript[λ, 2])                                                                                                           (33)

従来，わが国ではL積率法 （PWM法） の積率法に対するバイアスや誤差分散が小さいことが評価されてきているが，さらなる利点として，積率法では標本歪み係数のとりうる限界値が標本サイズによって変化するので標本サイズに関係せずに種々の分布モデルを比較
することができないのに対し，L積率法では下図 のように種々の分布モデルのL - skewnessとLkurtosisが一枚の図に表されるので，標本のL - skewness τ3 とL - kurtosis τ4 を知れば，概ねどの分布に近いかどうかが判断できる利点が有る．

2 母数の分布は点で，3 母数の分布は線で示されている．着色部分は3次および4次のL積率の
		存在する範囲である．各記号は，それぞれ，E：指数分布，G：Gumbel分布，L：Logistic
		分布，N：正規分布，GLO：generalized logistic分布，GEV：一般極値分布，GPA：一般
		パレート分布，LN3：3 母数対数正規分布，Pearson 3：Pearson type III 分布を示す

事例（Gumbel分布）

Clear[α,ξ]

まず，モデルのL積率を求める．Gumbel分布は2母数なので2次までのL積率が必要である．このため，(Subscript[β, 0])および(Subscript[β, 1])を求める．式(23)でr=0より

(Subscript[β, 0])=(∫_0^1)x(u)(u^0)du =(∫_0^1)x(u)du

F=Exp [-Exp [-((x-ξ)/α)]]より逆関数x(u)を求めると次式となる．

x (u) = ξ - α Log[-Log[F]]

λ1=β0=(∫_0^1)(ξ- α Log[-Log[u]])\[DifferentialD]u

なお，累積分布関数をu = F[x] とおくとdu/dx=f(x)となり，du=f(x)dx．結局，積分領域[0,1]が[-∞,∞]に変わって次式となり，積率法の平均と同じものが得られることになる．

(∫_(-∞)^∞)x (∂_x) Fgum[x] dx  = (∫_(-∞)^∞)x fgum[x] dx

Fgum[x]:=Exp [-Exp [-((x-ξ)/α)]]

次に (Subscript[β, 1])について求める. r=1として，

(Subscript[β, 1])=(∫_0^1)x(u)(u^1)du =(∫_0^1)x(u)udu

β1=(∫_0^1)(ξ- α Log[-Log[u]])u\[DifferentialD]u

(Subscript[λ, 2])=2(Subscript[β, 1])-(Subscript[β, 0]) より

λ2=2β1-β0

次に標本のL積率推定量を求める.まず，標本をSortして順序統計量を求め，次に 式(24),(25)を計算し，最後に(Subscript[λ, 1]), (Subscript[λ, 2])を求める．

いろいろな標本のL積率推定量を求めることになるので式 (23)~式 (33) をまとめて求めるモジュ - ル Lmoment を作成する.

Lmoment::usage="Lmoment[data_List]は,資料dataのL積率推定量を順に{λ1,λ2,λ3,λ4,λ2/λ1,λ3/λ2,λ4/λ2}のリストで返す.";

Lmoment[data_List]:=Module[{X,n=Length[data],B0,B1,B2,B3,λ1,λ2,λ3,λ4},
(* L moment method from HandBook of Hydrology Chap.18 FREQUENCY ANALYSIS OF EXTREME EVENTS *)
(* unbiased PWM estimators *)
	X=Sort[data,Greater];
B0:=Mean[X];
B1:=(∑_(j=1)^(n-1))(((n-j)X[[j]])/(n(n-1)));
B2:=(∑_(j=1)^(n-2))(((n-j)(n-j-1)X[[j]])/(n(n-1)(n-2)));                                         (* 18.1.14 *)
B3:=(∑_(j=1)^(n-3))(((n-j)(n-j-1)(n-j-2)X[[j]])/(n(n-1)(n-2)(n-3)));
(* L-moment estimators *)
λ1=B0;
λ2=2B1-B0;
λ3=6B2-6B1+B0;
λ4=20B3-30B2+12B1-B0;                                        (* 18.1.16 *)
{λ1,λ2,λ3,λ4,λ2/λ1,λ3/λ2,λ4/λ2}];

{標本λ1,標本λ2}=Lmoment[AMS日吉2day[[All,2]]][[;;2]]

prmLmom=Solve[λ1==標本λ1&&λ2==標本λ2,{α,ξ}]

L積率法gumPDF=Plot[Evaluate[fgum[x]/.prmLmom],{x,0,700},PlotRange->All,PlotStyle->Orange]

Show[ヒストグラムと積率法の確率密度分布,L積率法gumPDF,PlotRange->{{0,700},All}]

母数推定法(3)　L積率法　GEV

事例（GEV分布）

確率密度関数と確率密度関数と累積分布関数を再掲する.

f (x)=(1/α)(({1-((k (x-ξ))/α)})^((1/k)-1))exp[-(({1-((k (x-ξ))/α)})^(1/k))]          (7)

F(x)=exp(-(((1-((k (x-ξ))/α)))^(1/k)))                                                         (8)

Fgev[x]:=(E^(-(((1-((k (x-ξ))/α)))^(1/k))))

fgev[x]:=(∂_x) Fgev[x]

fgev[x]

GEVは3母数なので,3次までのL積率が必要になり,式(26)(Subscript[よりβ, 2])まで必要.　上記Gumbel分布の時の値をクリアしてから行う．

Clear[λ1,β0,model平均,λ2,β1,model分散]

まず，(Subscript[β, 0])=(∫_0^1)x(u)(u^0)du =(∫_0^1)x(u)du　であるので，式(8)から逆関数x(u)を求める

x(u)=ξ+(α/k)(1-(((-Log[u]))^k))

β0=(∫_0^1)(ξ+(α/k)(1-(((-Log[u]))^k)))\[DifferentialD]u

λ1=β0=Integrate[(ξ+(α/k)(1-(((-Log[u]))^k))),{u,0,1},Assumptions -> {{x, ξ} ∈ Reals, α > 0, k >-1}]

(∫_0^1)((α+k ξ-α (((-Log[u]))^k))/k)u\[DifferentialD]u

β1=Integrate[(ξ+(α/k)(1-(((-Log[u]))^k)))u, {u,0, 1}, Assumptions -> {{x, ξ} ∈ Reals, α > 0, k >-1}]

λ2=2β1-β0

(∫_0^1)((α+k ξ-α (((-Log[u]))^k))/k)(u^2)\[DifferentialD]u

β2=Integrate[(ξ+(α/k)(1-(((-Log[u]))^k)))(u^2),{u,0, 1}, Assumptions -> {{x, ξ} ∈ Reals, α > 0, k >-1}]

λ3=6β2-6β1+β0

Simplify[%]

%/λ2

Factor[%]

((分子分母に2)^(-k))を掛けると分子は2 (3^(-k))-3 (2^(-k))+1, 2 (3^(-k))-3 (2^(-k))+3-2,-2(1-(3^(-k))) +3(1-(2^(-k)))より

τ3=2(1-(3^(-k)))/(1-(2^(-k)))-3

λ1, λ2, τ3からα, ξ, kを求めることになるが, ちょっと難しい.よく使われている近似解法を表3に示した．

母数推定法(4)　POTに用いる分布

POTに用いる確率分布（再掲）

c) 指数分布 (Exp)

f (x)=(1/α)exp {-((x-ξ)/α)}                                                                         (11)

F(x)=1-exp {-((x-ξ)/α)}                                                                       (12)

変数xの分布範囲 :  ξ<x<∞

d) 一般パレート分布 (GPまたはGPD)

f (x)=(1/α) (({1-k((x-ξ)/α)})^(1/k-1))                                                               (13)

F(x)=1-(((1-k((x-ξ)/α)))^(1/k))                                                                  (14)

変数xの分布範囲:ξ<x≤ξ+α/k      (k > 0

ここに ξ : 位置母数（閾値），α : 尺度母数，k : 形状母数であり，k＝0の場合が指数分布である．

AMSとPOTそれぞれに用いる分布のL積率解

AMSに用いる分布とPOTに用いる分布のL積率解

具体的事例

対象資料　日吉ダム上流流域平均2日雨量

Gumbel確率紙へのプロット

点でプロット

GumbelPaperplot::usage="GumbelPaperplot[X_List]は区間最大値のリストXをCannaneのPlotting Positionを用いて,グンベル確率紙にプロットする.オプション(default値)はPlotting Positionを決めるα(0.4),横軸の名称(Rainfall [mm]),プロットの色(Red),縦軸のラベルを回転するかどうかのスィッチ(False)である.縦軸のラベルを回転すると右側のリタ-ンピリオドは示されず,縦にスリムな図となる.";GumbelPaperplot[x_List,α_:0.4,vUnit_:"Rainfall [mm]",col_:Black,Rotatelbl_:False]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,Tlb,yp},
ynorm[data_]:=-Log[-Log[data]];
		orderx=Sort[x];ytick={99.5,99,98,95,90,80,70,50,30,20,10,5,1};
lab1=Transpose[{ynorm[ytick/100],ytick}];
Tnen={200,100,50,20,10,5};
lab2=If[Rotatelbl,None,Transpose[{Take[ynorm[ytick/100],Length[Tnen]],Tnen}]];Tlb=If[Rotatelbl,None,"T"];
yp=Table[N[(i-α)/(n+1-2α)],{i,n}];      grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotRange->{{Min[x]/2,1.6Max[x]},{-1.54,6}},PlotStyle->col,AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick/100]},Frame->True,Axes->None,FrameTicks->{Automatic,lab1,Automatic,lab2},FrameLabel-> {vUnit, "F[%]",None,Tlb},Prolog->AbsolutePointSize[4],RotateLabel->Rotatelbl,PlotLabel->"Gumbel Probability Paper"]
	];

点ではなくマ - クでプロット

gumbellistplotmarknew[x_List,α_,vUnit_,mark_:{m1,0.04},xrange_:{0,0},Rotatelbl_:False]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,yp,Tlb},
ynorm[data_]:=-Log[-Log[data]];
		orderx=Sort[x];ytick=0.01{99.5,99,98,95,90,80,70,50,30,20,10};
lab1=Transpose[{ynorm[ytick],ytick}];
Tnen={200,100,50,20,10,5};
lab2=If[Rotatelbl,Transpose[{Take[ynorm[ytick],Length[Tnen]],Tnen}],None];Tlb=If[Rotatelbl,"T",None];
yp=Table[N[(i-α)/(n+1-2α)],{i,n}];      grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotMarkers->mark,PlotRange->{If[xrange=={0,0},{Min[x]/2,1.6Max[x]},xrange],{-1.54,6}},AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick]},Frame->True,Axes->None,FrameTicks->{Automatic,lab1,Automatic,lab2},FrameLabel-> {vUnit, "F",None,Tlb},Prolog->AbsolutePointSize[4],RotateLabel->Rotatelbl,PlotLabel->"Gumbel Probability Paper"]
	];

PlottingPosition[i_,n_,α_]:=((i-α)/(n+1-2 α));(*Weibull:0, Cunnane:0.4,Hazen:0.5*)

指数確率紙へのプロット

(*データをマークで指数確率紙にプロット*)POTExpPlotMark[x_List,title_:"",α_:0.4,vUnit_:"Rainfall [mm]",mark_:{m1,0.03},xrange_:{0,0},sz_:3.5 72]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,Tlb,yp},
ynorm[data_]:=-Log[1-data];
		orderx=Sort[x];ytick=0.01{99.98,99.95,99.9,99.8,99.5,99,98,95,90,80,70,50};
lab1=Transpose[{ynorm[ytick],ytick}];
Tnen={5000,2000,1000,500,200,100,50,20,10,5};
lab2=If[False,Transpose[{Take[ynorm[ytick],Length[Tnen]],Tnen}],None];Tlb=If[False,"T",None];
yp= Table[N[(i-α)/(n+1-2α)],{i,n}];      grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotMarkers->mark,PlotRange->{If[xrange=={0,0},{Max[0,Min[x]-0.1(Max[x]-Min[x])],Max[x]+0.2(Max[x]-Min[x])},xrange],{0.001,Max[ynorm[yp]]+1.0}},AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick]},Axes->None,Frame->True,FrameTicks->{Automatic,lab1,Automatic,lab2},FrameLabel-> {vUnit, "G",None,Tlb},ImageSize->sz,RotateLabel->False,PlotLabel->"Exponential Probability Paper",PlotRangePadding->None,Epilog->{Black,Text[Style[title,FontSize->12],ImageScaled[{0.19,0.86}],{Left,Top},Background->White]}]
	];

square={{-0.5,-0.5},{0.5,-0.5},{0.5,0.5},{-0.5,0.5},{-0.5,-0.5}};triangl={{0,1},{Cos[7Pi/6],Sin[7Pi/6]},{Cos[11Pi/6],Sin[11Pi/6]},{0,1}};

(({{m1, m2, m3, m4, m5, m6}, {m1r, m2r, m3r, m4r, m5r, m6r}, {m1b, m2b, m3b, m4b, m5b, m6b}, {m1g, m2g, m3g, m4g, m5g, m6g}, {m1p, m2p, m3p, m4p, m5p, m6p}, {m1c, m2c, m3c, m4c, m5c, m6c}, {m1m, m2m, m3m, m4m, m5m, m6m}, {m1br, m2br, m3br, m4br, m5br, m6br}, {m1o, m2o, m3o, m4o, m5o, m6o}}))=Graphics/@{{#,Circle[{0,0},1]},{#,Disk[{0,0},1]},{#,Line[square]},{#,Polygon[square]},{#,Line[triangl]},{#,Polygon[triangl]}}&/@{Black,Red,Blue,Green,Pink,Cyan,Magenta,Brown,Orange};

Gumbel確率紙=gumbellistplotmarknew[AMS日吉2day[[All,2]],0.4,"2日降水量 (mm/2day)"]

母数推定，確率評価，確率水文量

表3を用いて標本から母数を推定する．AMSの場合，AMS_Listは標本の区分最大値のリストである．POTの場合，POT_Listは標本の閾値を超える独立な極大値のリストである．閾値の設定方法としては種々提案されているが，ここでは標本平均超過関数(SMEF)を使う．
ここまで，すべての確率分布において非超過確率をFで表してきたが，AMSの確率評価では1-Fが年超過確率を表すが，POTの確率評価の場合は観測資料全体について閾値を超える事象の確率を示しているだけで時間の概念がないので直接比較することができない．これ以降，このことを明確に区別するため，POTの閾値を超過する事象の非超過確率をG(x)で表す

AMS  グンベル分布

Fgum[x]:=Exp [-Exp [-((x-ξ)/α)]];

Gumbelprm::usage="Gumbelprm[AMS]は区間最大値のリストAMSにフィットするグンベル分布の母数{α,ξ}を返す.";Gumbelprm[AMS_List]:=Module[{α,ξ,λ},
				 λ=Lmoment[AMS];
          α=λ[[2]]/Log[2.](* (Subscript[λ, 2])/Log2 *);
          ξ=λ[[1]]-0.5772α (* (Subscript[λ, 1])-γα *);
     	{α,ξ}];

GumbelF::usage="GumbelF[AMS,x]は区間最大値AMSに対して母数推定し,レベルxに対応する非超過確率を返す";
GumbelF[AMS_List,x_]:=Module[{α,ξ},{α,ξ}=Gumbelprm[AMS];
Exp [-Exp [-((x-ξ)/α)]]];

GumbelXq::usage="GumbelXq[AMS,ReturnPeriod]は区間最大値AMSに対して母数推定し,ReturnPeriodに対応する非超過確率を返す.ReturnPeriodは単独でもリストでもかまわない.";
GumbelXq[AMS_List,ReturnPeriod_]:=Module[{α,ξ},{α,ξ}=Gumbelprm[AMS];
	ξ-α  Log[-Log[1-(1/ReturnPeriod)]]] ;

AMS  一般極値分布(GEV)

Fgev[x]:=(E^(-(((1-((k (x-ξ))/α)))^(1/k))))

GEVprm::usage="GEVprm[AMS_List]は区間最大値のリストAMSにフィットするGEV分布の母数{α,κ,ξ}を返す.";GEVprm[AMS_List]:=Module[{κ,α,ξ,c,λ},
			 λ=Lmoment[AMS];
   		c=((2λ[[2]])/(λ[[3]]+3λ[[2]]))-((Log[2])/(Log[3]));
        κ=7.859c + 2.9554c^2;                                                   (* 18.2.22a *)
        α=((κ  λ[[2]])/(Gamma[1+κ](1-(2^(-κ)))));
        ξ=λ[[1]]+((α( Gamma[1+κ]-1))/κ);
		{α,κ,ξ}];

GEVF::usage="GEVF[AMS,x]は区間最大値AMSに対して母数推定し,レベルxに対応する非超過確率を返す";GEVF[AMS_List,x_]:=Module[{α,κ,ξ},
		{α,κ,ξ}=GEVprm[AMS];Exp[-(((1-κ((x-ξ)/α)))^(1/κ))]];

GEVXq::usage="GEVXq[AMS,ReturnPeriod]は区間最大値AMSに対して母数推定し,ReturnPeriodに対応する確率水文量を返す.ReturnPeriodは単独でもリストでもかまわない.";
GEVXq[AMS_List,ReturnPeriod_]:=Module[{α,κ,ξ},
		{α,κ,ξ}=GEVprm[AMS];
		ξ+(α/κ)(1-(((-Log[1-(1/ReturnPeriod)]))^κ))] ;

POT  指数分布 (Exp)

Expprm::usage="Expprm[POT_List]は閾値超過資料のリストPOTに,閾値ξとしてMin[POT]を用い,Exp分布の母数{α,ξ}を返す.Expprm[POT_List,threshold]の形で閾値thresholdが指定された場合はそれをξとして,Exp分布の母数{α,ξ}を返す.";
Expprm[POT_List,threshold_:-9999]:=Module[{α,ξ},
  ξ=If[threshold==-9999, Min[POT],threshold]; α=Mean[POT]-ξ;
		{α,ξ}];

ExpG::usage="ExpG[POT,x,threshold]は閾値超過資料のリストPOTに対してthresholdを踏まえて母数推定し,レベルxに対応する非超過確率を返す";
ExpG[POT_List,x_,threshold_:-9999]:=Module[{α,ξ},{α,ξ}=Expprm[POT,threshold];1-Exp[-((x-ξ)/α)]];

ExpXq::usage="ExpXq[POT,G,threshold]は閾値超過資料のリストPOTに対してthresholdを踏まえて母数推定し,非超過確率Gに対応する確率水文量を返す.";
ExpXq[POT_List,G_,threshold_:-9999]:=Module[{α,ξ},
		{α,ξ}=Expprm[POT,threshold];
		ξ-α Log[1-G]] ;

POT  一般パレ - ト分布(GP)

GPprm::usage="GPprm[POT_List]は閾値超過資料のリストPOTに,閾値ξとしてMin[POT]を用い,GP分布の母数{α,κ,ξ}を返す.GPprm[POT_List,threshold]の形で閾値thresholdが指定された場合はそれをξとして,GP分布の母数{α,κ,ξ}を返す.";
GPprm[POT_List,threshold_:-9999]:=Module[{κ,α,ξ,λ},  ξ=If[threshold==-9999, Min[POT],threshold];
			 λ=Lmoment[POT];
   		κ=((λ[[1]]-ξ)/(λ[[2]]))-2.;
        α=(λ[[1]]-ξ)(1+κ);
		{α,κ,ξ}];

GPG::usage="GPG[POT,x,threshold]は閾値超過資料のリストPOTに対してthresholdを踏まえて母数推定し,レベルxに対応する非超過確率を返す";
GPG[POT_List,x_,threshold_:-9999]:=Module[{α,κ,ξ},
		{α,κ,ξ}=GPprm[POT,threshold];
		1-(((1-κ(((x-ξ))/α)))^(1/κ))] ;

GPXq::usage="GPXq[POT,G,threshold]は閾値超過資料のリストPOTに対してthresholdを踏まえて母数推定し,非超過確率Gに対応する確率水文量を返す.";
GPXq[POT_List,G_,threshold_:-9999]:=Module[{α,κ,ξ},
{ α,κ,ξ}=GPprm[POT,threshold];
		ξ+(α/κ)(1-(((1-G))^κ))];

POT  標本平均超過関数(SMEF)

テキストの使い方の説明部分が誤っている.下記のように訂正!

この指標による判断は，式 (41) の値が一定となる下限の u である．すなわち下限の u の値を閾値に使うことになる．

SampleMEFatu::usage="SampleMEFatu[POTraw_List,u]は閾値が小さい閾値超過資料のリストPOTrawのうちu以上の要素の平均を返す.POTrawの最大値よりuが大きい場合はNullを返す.";
SampleMEFatu[POTraw_List,u_]:=If[Sort[POTraw,Greater][[1]]≥u,Mean[Select[POTraw-u,#≥0&]],Null];(*for Real u*)

具体的事例

Gumbelprm[AMS日吉2day[[All,2]]]
GEVprm[AMS日吉2day[[All,2]]]

GumbelF[AMS日吉2day[[All,2]],500]

GumbelXq[AMS日吉2day[[All,2]],100]

GumbelXq[AMS日吉2day[[All,2]],{100,200,600}]

GEVF[AMS日吉2day[[All,2]],500]

GEVXq[AMS日吉2day[[All,2]],100]

GEVXq[AMS日吉2day[[All,2]],{100,200,600}]

十分低いレベルの閾値 (62 mm/day) で抽出したPOT資料

年毎POTcr62

DateListPlot[Length/@年毎POTcr62,{1952},PlotLabel->"日吉ダム2日雨量 50mm/2day超過日数 : "<>ToString[MannKendall[Length/@年毎POTcr62]]]

DateListPlot[Total/@年毎POTcr62,{1952},PlotLabel->"日吉ダム2日雨量 62mm/2day超過雨量年合計 : "<>ToString[MannKendall[Total/@年毎POTcr62]]]

POTExpPlotMark[Flatten[年毎POTcr62],"日吉ダム2日POT",0.4,"Rainfall [mm]",{m2r,0.03}]

SMEFを描く

Plot[SampleMEFatu[Flatten[年毎POTcr62],u],{u,50,200},
Frame->True]

この図から閾値は109mmと設定する.

threshold=109

pot=Select[Flatten[年毎POTcr62],#≥threshold&]

Length[%]

Expprm[pot,threshold]

GPprm[pot,threshold]

年数=Length[年毎POTcr62]

AMSとPOTの関係

豪雨や洪水の発生がポアソン過程に従う場合，G(x)は年最大値の非超過確率F(x)と次式で結びつけられる．

F (x)=exp {-λ (1-G (x))}

ここに，λは1年あたりの事象の平均的な発生率である．

G (x)として，モデルの確率分布を用いるとAMSの確率分布と比較できる．また，標本から得られる経験分布を用いるとAMSの確率プロットと比較できる．ここでは，両者を1枚の年最大値の確率プロットに描いて比較する．

1)AMSの確率分布をAMS標本と比べる

Gumbelplot[AMS_List,col_:Orange]:=Module[{α,ξ,xstrt,Fgumbel,x},
	Fgumbel[x_]:=Exp[-Exp[-((x-ξ)/α)]];
		{α,ξ}=Gumbelprm[AMS];
			xstrt=Min[Min[AMS]-(Max[AMS]-Min[AMS])/8,0];
		Plot[ -Log[-Log[ Fgumbel[x]]],{x,xstrt,Max[AMS]+(Max[AMS]-Min[AMS])/4},PlotStyle->col]];

GEVplot[AMS_List,col_:Purple,prmsw_:False]:=Module[{κ,α,ξ,prm,xstrt,GEV,x},
		GEV[x_]:=Exp[-(((1-((κ(x-ξ))/α)))^(1/κ))];    (* κ≠0 18.2.18 *)
		{α,κ,ξ}=GEVprm[AMS];
				If[prmsw,Print["GEV parameters  : α: ",α," κ: ",κ," ξ: ",ξ]];
			xstrt=Min[Min[AMS]-(Max[AMS]-Min[AMS])/8,0];
	Plot[ -Log[-Log[ GEV[x]]],{x,xstrt,ξ+(α/κ)(1-(((-Log[.9975]))^κ))},PlotStyle->col]];

2) POT資料をGumbel確率紙にプロットする

GumbelpaperPOTtoAME[x_List,ny_,α_:0.4,unit_:"Rainfall [mm]",mark_:{m2r,0.03},xrange_:{0,0},imgsz_:4 72,title_:"",tpos_:{0.3,0.83}]:=Module[{n=Length[x],ynorm,Tnen,grf,orderx,ytick,lab1,lab2,yp},
ynorm[data_]:=-Log[-Log[data]];
		orderx=Sort[x];ytick=0.01{99.5,99,98,95,90,80,70,50,30,20,10,5,1};
lab1=Transpose[{ynorm[ytick],ytick}];
Tnen={200,100,50,20,10,5};
lab2=Transpose[{Take[ynorm[ytick],6],Tnen}];
yp= PlottingPosition[Range[n],n,α];    
				yp=Exp[-(n/ny)(1-yp)];  grf=ListPlot[Transpose[{orderx,ynorm[yp]}],PlotMarkers->mark,PlotRange->{If[xrange=={0,0},{Max[0,Min[x]-0.1(Max[x]-Min[x])],Max[x]+0.2(Max[x]-Min[x])},xrange],{-1.54,6}},AspectRatio->1.3,GridLines->{Automatic,ynorm[ytick]},Frame->True,Axes->None,FrameTicks->{Automatic,lab1,Automatic,lab2},FrameLabel-> {unit, "F",None,"T"},RotateLabel->False,ImageSize->imgsz,PlotLabel->"Gumbel Probability Paper",Epilog->{Black,Text[title,ImageScaled[tpos],{-1,0},Background->White]}]
	];

(*nyの情報を基に指数分布をGumble確率紙にプロット*)ExpPlotonGumPaper[POT_List,ny_,threshold_:-9999,col_:{Orange,Dashed}]:=Module[{n=Length[POT],α,ξ,prm,xstrt,Gexp,sp,x},
	Gexp[x_]:=1-Exp[-(((x-ξ)/α))];
		sp[p_]:=-Log[-Log[p]];
	{α,ξ}=Expprm[POT,threshold];
	Plot[ -Log[(n/ny)(1-Gexp[x])],{x,ξ,ξ+α sp[0.9999]},PlotStyle->col]];

(*nyの情報を基に一般パレート分布をGumble確率紙にプロット 描画範囲を形状母数で適正に判断するよう修正20131007*)GPPlotonGumPaper[POT_List,ny_,threshold_:-9999,col_:{Purple,Dashed}]:=Module[{n=Length[POT],α,κ,ξ,Ggp,sp,x},
	Ggp[x_]:=1-(((1-((κ(x-ξ))/α)))^(1/κ));
		sp[p_]:=-Log[-Log[p]];
{ α,κ,ξ}=GPprm[POT,threshold];
Plot[ -Log[(n/ny)(1-Ggp[x])],{x,ξ,Which[κ>0,ξ+(α/κ) ,κ≤0,3Max[POT]]},PlotStyle->col]];

potonGumbelPP=GumbelpaperPOTtoAME[pot,年数]

AMSのプロットにGumbel分布とGEVを描いて比較する

Show[Gumbel確率紙,Gumbelplot[AMS日吉2day[[All,2]]],GEVplot[AMS日吉2day[[All,2]]]]

Show[Gumbel確率紙,potonGumbelPP,Gumbelplot[AMS日吉2day[[All,2]]],GEVplot[AMS日吉2day[[All,2]]],ExpPlotonGumPaper[pot,年数,threshold],GPPlotonGumPaper[pot,年数,threshold]]

別に配布したCDFファイル”閾値によるAMSとPOTの関係変化.cdf”をみてください．このファイルはMathematicaを持っていなくてもCDFplayerをダウンロード(Free)すればMathematicaの計算機能が働いて動きます．閾値と年のスライダ - を動かしてみてください

スライダ - の右の[+] をクリックするとコントロ - ラ - が出ます.

閾値の違いによるPOTとAMSの差がクリアになります

まとめ

水文頻度解析の歴史的経過踏まえてトレンド評価およびL積率法による水文頻度解析手法を紹介した.

AMSは簡便であるが，AMSのみではわからない問題があり，これは，AMSだけを見ていると判断ができない．

POTとAMSを比較して同じ時系列資料を見るとその違いが明瞭になる．

POTの閾値は標本平均超過関数により決める方法が有効である．

POTを用いる方法は，年数が限られているGCM資料を用いる際にも有効である．