ECOH/YG 技術論文 No. 13

2010年10月11日版

L-モーメント法を用いた高波の極値統計解析におけるプロッティングポジション公式の選択と各種統計量の信頼区間
Selection of Plotting Position Formulas in the L-moment method applied for Extreme
Wave Analysis and Confidence Interval of Statistical Estimates

合田良実
Yoshimi GODA
土木学会名誉会員

工博 横浜国立大学名誉教授(株) エコー（〒110-0014 東京都台東区北上野 2-6-4）

A Monte Carlo simulation study has been carried out on the reliability of the parameter and quantile
estimation of the GEV, GPA, and Weibull distributions by means of the L-moment method for application to
extreme wave statistics. The three population distributions are assigned with five to six values of the shapeparameter each, and the scale and location parameters are so specified as to yield the 10-year and 50-year wave
heights of 8.33 and 10.0 m, respectively. Samples of extreme wave heights are generated from 16 populations
at the size varying from 10 to 1000, and 50,000 samples for respective conditions are analyzed for statisticalcharacteristics. The unbiased plotting position recommended by Hosking and Wallis (1997) is confirmed toyield minimal biases on the parameter estimates, but it yields small biases on quantile estimates. A new plotting
position estimator is proposed, which yields minimal quantile biases. A set of empirical formulas are derived
for estimation of the coefficient of variations of the L-moment ratios and the quantile of the three distributions.

Key Words: Extreme wave statistics, GEV distribution, GPA distribution, Weibull distribution, Lmoments, plotting position, quantile variability(1) を示し，定数値として  = 0,  = 0 の組み合わせと，

１．まえがき = 0.35,  = 0 の二つを提示している．
i 
Pi 
n

先に著者ら (2009a, 2009b) は極値波浪統計解析に対する L-モーメント法の適用法について検討し，Hosking・Wallis (1997) が取り扱っていないワイブル分布への適用について解説した．そして，日本沿岸の波浪観測地点のうち 8 個所を選び，L-モーメント 法を使って確率
波 高 の 推 定 を 行 っ た ． さ ら に 著 者 ら (2010) は
NOWPHAS のデータベースから日本海沿岸11地点を選
び，地域頻度解析の手法を適用してこの沿岸域の確率
波高の推定を行った．
Hosking・Wallis (1997) が紹介している L-モーメン
ト法による極値解析では，L-モーメントの計算に当たって使用するプロッティングポジションとして式(1)

ここに，i は昇順に並べ替えられた極値データの順位，
n は極値データの総数，すなわち標本の大きさである．
一般化極値分布 (Generalized Extreme Value distribution: GEV) に対して確率重み付き積率法 (Probabilityweighted moments: PWM) を適用する場合については，
Hosking ほか (1985) が = 0.35,  = 0 の組み合わせが良い結果を与えると報告している．また，一般化パレー
ト分布 (Generalized Pareto distribution: GPA) に対しても，

 = 0.35,  = 0 の組み合わせが適切であることを
Hosking・Wallis (1987) が報告している．しかし，その-1-

後の Hosking・Wallis (1995) の検討では，  = 0,  = 0 のモーメント と2 および L-モーメント比3 と 4 の理組み合わせのほうが優れているとして，この組合せを論値ならびに母数の推定式を導いた。なお，1 次 L-

不偏公式 (unbiased plotting position) と呼んでいる．このモーメント  は平均値であり，1 次 L-モーメントと 2検討結果に基づいて，Hosking・Wallis (1997) は L-モー

次 L- モーメントとの比  = 1 /  は L-Coeffcient of

メント法の適用に際してプロッティングポジションとVariation (L-CV) と名付けられている．また，L-モーメ

して不偏公式を用いることを推奨している．

ント比3 は L-skewness，L-モーメント比4 は L-kurtosis

極値統計解析に最小2乗法を用いる場合には，推計と呼ばれるので，本論文でもこれらの名称を使用する．

値の偏りが生じないようにプロッティングポジション

L-モーメントおよび L-モーメント比ならびに分布関数について吟味し，極値データに当てはめる分布関数ごの母数推定式は， Hosking・Wallis (1997) の解説書にもとに最適のプロッティングポジション公式を選択する． 付録としてまとめられている。そのうちから，高波のこれについては合田 (2008, 340-343頁) が詳しく述べて極値統計解析に使用する分布関数に関わるものを以下いる．

に再録する。ただし，ワイブル分布については前論文

L-モーメント法の母体となった PWM 法では，分布（ 合 田 ら 2009a, 2009b ） に お い て Greenwood ほ か関数の母数や再現確率量の推計値の精度について多く(1978) を参照して独自に導いたものである。

の研究が行われており，推計値の真値からの偏りなど(A) 指数型分布が吟味されてきた．それらによると，標本の大きさが100程度以下ではかなり大きな偏りを示すことが珍しa) 分布関数：F ( x )  1  exp{( x  B ) / A},

くない（たとえば，Landwehr ほか (1979a, 1979b)，b) 確率密度関数：

Boes ほか (1989) など）．

Bxf ( x)  1 / A  exp{( x  B) / A}

しかしながら， PWM 法ならびに L-モーメント法に

(2)
(3)

c) 平均値と分散：関する研究では，プロッティングポジション公式の選

E[ x]  B  A, Var[ x]  A 2

定についてあまり吟味されていないようである．そこd) 確率統計量 (quantile)：

で本研究では， L-モーメント法で使用するのに適切な

x  F 1 ( P)  B  A ln(1  P)e) L-モーメント：

プロッティングポジション公式を見いだすことを目的として，モンテカルロ法による大規模な数値実験を行

(4)
(5)

1  B  A,  2  (1 / 2) A,  3  1 / 3,  4  1 / 6 (6)うものである．

f) 母数推定式：

全体の構成は次の通りである．まず， 2.においてこAˆ  2l2 ,こで使用する 3 種類の極値分布の関数形ならびにL-

Bˆ  l1  Aˆ(7)

モーメント比を紹介する．L-モーメント法そのものについては先の著者らの論文 (2009a) で記述したので，

ここで，F ( x ) は非超過確率を表し，A は尺度母数，B

ここでは割愛する．次に 3.では，L-モーメント 法を用は位置母数と呼ばれるものである．指数型分布は後述

いて極値統計解析を行う手順を説明し，4.で数値実験の一般化パレート分布で形状母数が k = 0，ワイブル分の対象とする極値分布や計算条件などの概要を説明す
る．数値計算の結果は，まず 5. において L-モーメント

布で k = 1.0 の場合に相当する．なお， Â , B̂ の ˆ の記
号は推計値であることを示す．以下においても，推計の推計値の信頼性について記述し，極値分布の母数推値には ˆ の記号を付す．また，l1, l2 は標本について計

計については 6.，確率波高の推計値の偏り率について算した1 次と 2 次の L-モーメントを表す．

は 7.，確率波高の信頼区間については 8.で述べる．

(B) ガンベル分布（極値Ⅰ型分布）

また，既往の諸研究との比較を 9.で論じ，さらに 10.

a) 分布関数：F ( x )  exp[  exp{( x  B ) / A}],    x   (8)において地域頻度解析における地域分割の際の不均質性の判断基準を新しく提案し，先に検討した日本海沿b) 確率密度関数：

岸の高波への適用例を示す．最後に全体とまとめと今f ( x) 

後の課題を11. に述べる．

 xB1
 x  B exp exp
AA
 A 
(9)c) 平均値と分散：

２．極値分布関数と L-モーメント比の関係

E[ x]  B  A , Var[ x]   2 A2 /6

Hosking (1990) は，各種の極値分布関数に対する L-

2

(10)

また，式 (19) における形状母数 k は近似推定式であd) 確率統計量 (quantile)：
1(11)x  F ( P)  B  A ln(1  ln P)e) L-モーメント：
1  B  A ,  2  A ln 2,  3  0.1699,  4  0.1504 (12)f) 母数推定式：Aˆ  l / ln 2, Bˆ  l  Aˆ 
2り，Hosking (1990) が導いたものである。

(D) 一 般 化 パ レ ー ト 分 布Distribution : GPA)a) 分布関数：(13)1

1  [1  k ( x  B) / A]1 / k 1  k  0 : B  x  B  A/k
F ( x)  k  0:
Bx
1  exp{( x  B) / A}, k  0 :Bx

ここに  はオイラーの定数, 0.5772....である。ガンベル
分布は次項の一般化極値分布で形状母数が k = 0 の場合に相当する．
(C) 一 般 化 極 値 分 布(20)(Generalized Extreme-Valueb) 確率密度関数：

Distribution: GEV)a) 分布関数：f ( x)  [1  k ( x  B) / A]1 / k 1 / A : B  x, k  0 (21)

exp[{1  k ( x  B) / A}1 / k ], k  0 :    x  B  A/k
F ( x)  k  0 : B  A/k  x  exp[ exp{( x  B) / A}],k  0:    x  c) 平均値と分散：

E[ x]  B  A /(1  k )Var[ x]  A /[(1  k ) (1  2k )]

(14)b) 確率密度関数：f ( x)  [1  k ( x  B) / A]1 / k 1
 exp{[1  k ( x  B) / A]

1/ k

/A : k 0

2

E[ x]  B  A [1  (1  k )] / kVar[ x]  A [(1  2k )   (1  2k )] / k 222 B  A{1  (1  P) k } / k , k  0
x  F 1 ( P)  k 0 B  A ln(1  P),(16)

 B  A{1  ( ln P ) k } / k , k  0
x  F ( P)  k 0 B  A ln( ln P ),(17)f) 母数推定式：kˆ  (1  3 t 3 ) /(1  t 3 ),Aˆ  (1  kˆ) (2  kˆ) l 2 , Bˆ  l1  (2  kˆ) l 2

大・極小値に対する漸近関数ではないけれども，ポア(18)

ソン分布に従う確率変量のうち，ある閾値を超えるものに着目する場合などを表示するのに適切な分布関数f) 母数推定式：kˆ  7.8590c  2.9554c 2 , c (1  2 ) (1  kˆ)(25)

一般化パレート分布は一般化極値分布のような極 4  {5 (1  4 k )  10 (1  3 k )  6 (1  2 k )} /(1  2 k )

 kˆ

(23)

(24)

e) L-モーメント：
1  B  A{1  (1  k )} / k , 2  A(1  2 k ) (1  k ) / k , 3  2 (1  3 k ) /(1  2 k )  3,

Aˆ 

(22)

e) L-モーメント：
1  B  A /(1  k ),  2  A /{(1  k )(2  k )}, 3  (1  k ) /(3  k ),  4  (1  k ) ( 2  k ) /{(3  k ) (4  k )}d) 確率統計量 (quantile)：
12d) 確率統計量 (quantile)：(15)c) 平均値と分散：

l 2 kˆ

(Generalized Pareto

2
log 2

3  t 3 log 3といわれる。北野ほか (2002) は Kodiak 沖の高波デー
タについてこの分布の適用性を検討している。なお，(19)

この分布関数の特性等については Coles (2001) などを, Bˆ  l1  Aˆ {1  (1  kˆ)} / kˆ

参照されたい。

ここに，t3 は標本について計算した L-skewness である． (E) ワイブル分布
一般化極値分布は，式 (14) の表現では k < 0 のときにa) 分布関数：
極値 II 型分布，k > 0 のときに極値 III 型分布となり，F ( x)  1  exp[{( x  B) / A}k ], B  x  

前者では確率統計量が下限値を持ち，後者では上限値を持つ。また，k = 0 のときは極値 I 型分布である。な(26)b) 確率密度関数：

お，高波の極値統計で使われている極値 II 型は式 (14)f ( x)  k / A  [( x  B) / A] k 1の形状母数を k’ = –1/k と変換したものに等しい。

 exp[{( x  B) / A}k ] : B  x  

3

(27)

c) 平均値と分散：
E[ x]  B  A (1  1/ k )標本に当てはめる分布関数として複数の候補を取り22Var[ x]  A (1  2 / k )   (1  1/ k ) 上げるときには，標本にもっとも適合する関数を選ぶ(28)必要がある．著者ら (2009a, 2009b) は先に TUD 指標(Twenty-Up Deviation) を提案した．これは，標本中でd) 確率統計量 (quantile)：x  F 1 ( P)  B  A[{ ln(1  P)}1 / k ]大きなほうから20個のデータを取り上げ，当てはめた(29)

分布関数を用いたこれら 20 位のデータに対する確率統e) L-モーメント：

計量の推計値と観測値の相対誤差の 2 乗平均の平方根
値のもっとも小さいものを最適合と判断するもので

1  B  A(1  1 / k ),  2  A(1  2 1 / k ) (1  1 / k ),あった．

 3  3  2 (1  3 1 / k ) /(1  2 1 / k ),ここで，分布関数によるデータの推計値というのは，

 4  {5 (1  4 1 / k )  10 (1  3 1 / k )  6 (1  2 1 / k )} /(1  2 1 / k )一般に次のように表される．

(30)xˆi  F 1 ( Pi )  Bˆ  Aˆ y if) 母数推定式：k̂ = 285.3t36 – 658.6 t35 + 622.8 t34 – 317.2 t33+ 98.52 t32 – 21.256 t3 + 3.5160

式 (33) 右辺の yi は順位 i の非超過確率 Pi に対応する



Aˆ 

l2
(1  2 1 / k )(1  1 / kˆ)ˆ, Bˆ  l1  Aˆ (1  1 / kˆ)

(33)

(31)

基準化変量と称するもので，分布関数毎に次のように(32)与えられる．
1) 一般化極値分布 (GEV)：

なお，著者ら (2009b) の論文では，その中の表－2 に{1  ( ln Pi ) k } / k ,yi  
 ln( ln Pi ),記載した尺度母数 A の推定式が式 (32) の分子に k を乗じた形となっているけれども，それは校正漏れによる間違いである．ただし，論文中の計算結果はすべて式

k 0
k 0

(34)

2) 一般化パレート分布 (GPA)：(32) に基づいているので，訂正の必要はない．

[1  (1  Pi ) k ] / k , k  0
yi  
k 0
ln(1  Pi ),３．L-モーメント法による極値統計解析の手順(1) データの並べ替えとL-モーメントの計算(35)3) ワイブル分布：

高波の資料（標本）が得られたならば，標本中のyi  { ln(1  Pi )}1 / k

データを昇順に並べ替え，1～4 次の L-モーメント l1～(36)

l 4，とその比 t = l2/l1, t3 = l3/l2, および t4 = l4/l2 を計算する．

前述のように，GEV分布で k = 0 というのはガンベル

この計算には，Hosking が開発した FORTRAN プログ

分布であり，GPA分布で k = 0 というのは指数分布であ1

ラム lmoments をベースにしてプログラムを作成し，そるので，式 (34) と (35) はそれぞれガンベル分布と指数れを使用する．

分布を包含している．
最小 2 乗法では，観測されたデータ xi と基準化変量 yi(2) 当てはめ候補の分布関数の母数値の推定の間に式 (33) を仮定し，直線回帰式を適用して母数 A

標本の L-skewness [t3] が分かれば，式 (19), (25), (31)と B の値を推計している．

によって形状母数 k の値が推定できる．指数分布とガ

TUD 指標では上位 20 個のデータに限定したけれど

ンベル分布は形状が固定されているので，この手順はも，本研究では i = 1, 2, …, n の全データについて観測値省略される．

xi とその推計値 x̂i の間の相関係数 r ( xi , xˆi ) を計算し，尺度母数 Â は推定された k̂ 値と 2 次の L-モーメン

その上限値である 1 からの残差 r  1  r ( xi , xˆi ) が最小

ト l2 から求められ，位置母数 B̂ はこれらの推計値に1
次の L-モーメント l1 の情報を追加することで求められる．

となるものを最適合と判断する．いわば，最小残差法(3) 最適合の分布関数の選定(PPC)” とも呼ばれている（Choudhury et al. 1991）
合田 (1988) が以前に提案した MIR 法では，当てはめ

1(Minimum Residue of Correlation coefficient: MRC) を採用する．この方式は，“Probability Plot Correlation Method

る分布関数と標本の大きさに応じた残差の平均値をあhttp://lib.stat.cmu.edu からダウンロード可能である．

4

ばならない．

らかじめ算定しておき，残差平均値に対する比の大小
で最適分布を判定した．しかし，L-モーメント法では

1  1 / 50
y50  y10






1  1 / 50  
B  H 50 1 
 1  y10 / y50  A  H 50

相関係数の残差の平均値などを算定することができな
いため，残差の絶対値を判定の指標とするものである．
具体的には，観測された極値あるいはモンテカルロ
法によって発生させた極値の標本に対して複数の極値(39)

分布関数を当てはめ，それぞれについて相関係数

ここに，y10 と y50 は再現期間10年と50年に対する基準

r ( xi , xˆi ) の 1 からの残差 r  1  r ( xi , xˆi ) を算定する．

化変量であり，R = 10 および 50 として式 (37) で非超過

そして，当てはめ候補の分布関数の中で残差が最小の確率を求め，それを式 (34) ～ (36) に代入して算定されものを最適合と見なす方式である．

る．たとえば，裾長度を50 = 1.20，平均発生率を  = 2に設定すると，指数分布であれば H50 = 10.0 m に対して自動的に A = 1.036 m，B = 5.23 m と定まる．

(4) 確率波高の推計標本に対して当てはめた分布関数の形状・尺度・位(2) 発生させる高波の極値分布関数

置の 3 母数が推定されると，所定の再現期間 R に対す

各種の極値分布は， L-モーメント比 3 (L-skewness)

る確率波高 HR は式 (33) を用いて推定できる．このときの基準化変量 yi は，分布関数に応じて式 (34)～ (36)と4 (L-kurtoness) の相関関係によってその特性を表すによって計算する．その際の非超過確率 Pi = PR は次式ことができる．図-1 は，GEV，GPA，およびワイブル

で与えればよい．

分布の特性曲線を示したもので，図中の ◆，■，● のPR  11R

記号はそれぞれの分布関数で形状母数 k が特定の値を(37)とる場合を示している．

ここに， は極値解析を行った高波データの年間の平0.4

均発生率である．

４．数値計算の条件と計算方法

GumbelGEV (k < 0)
GEV (k > 0)
Weibull
Exponential
GPA (k > 0)
GPA (k < 0)

L-Kurtosis, 40.3(1) 発生させる高波の特性
数値シミュレーションによって高波の標本を多数作

-0.30

GPA:
k = -0.2

-0.1

-0.15

0.1
2.8

Weibull:k = 0.75

1.0

0.001
0.15

1.4
0.2

2.0

性として重要なのは，合田 (2002) が提唱した裾長度パ

0.4
0.6

ラメータである．これは 10 年確率波高 H10 に対する 5000

年確率波高 H50 の比として次のように定義される．

 50  H 50 /H 10

GEV:
k = -0.45

0.2

成する際には，高波として現地の波浪特性を反映させたものであることが必要である．高波の極値分布の特

L-Skewness vs.
L-kurtosis

0.1

0.2

0.3

0.4

0.5

L-Skewness, 

3

図-1 高波の母分布関数の L-kurtosis と L-skewness の関係(38)

わが国沿岸の高波の極値統計解析の結果では，裾長度

本研究では数値シミュレーションで発生させる高波パラメータ 50 が 1.13～1.27 の値をとっている（合田の母分布関数として GEV分布，GPA分布，およびワ

2008, 375頁）．そこで，本研究では裾長度パラメータ

イブル分布の 3 種類を考え，形状母数をそれぞれ 5 ～6として 50 = 1.10, 1.20, および 1.30 の 3 通りの値を用い通りに変えた，全体で 16 種類の分布関数を使用する．

ることとした．

分布関数と形状母数の組合せを表-1 に示す．これらのまた，50 年確率波高としては H50 = 10.0 m に設定す分布関数の 3 と 4 の関係は図-1 に表すとおりである．

る．これは数値の換算が容易なように区切りの良い数

なお，GEV 分布の k = 0.001 はガンベル分布を近似的値としたものである．

に表現したものである．また，ワイブル分布の k = 1.0

裾長度パラメータをある値に設定すると，分布関数

は指数分布であり，GPA 分布の k = 0 と共通である．たの尺度母数と位置母数は任意の値を与えることができだし，確認のために GPA 分布の k = 0.001 のケースも計ず，次式によって一義的に決定される値を使わなけれ算に加えた．

5

10

1

10

0

表-1 高波の発生に使用する分布関数

-0.45

-0.3

-0.15

0.001

0.15

–

一般化ﾊﾟﾚｰﾄ (GPA)

-0.2

-0.1

0.001

0.2

0.4

0.6

ワイブル

0.75

1.0

1.4

2.0

2.8

–

Probability density, f (H)形状母数 k

分布関数
一般化極値 (GEV)GEV Distribution(3) 高波の平均発生率と発生年数ならびに確率波高算
定の再現期間

10

-1

10

-2

10

-3

10

-4

k = -0.3
k = -0.15
k = 0.001
k = 0.15

極値の統計的変動性はもっぱら標本の大きさ，すな
わち標本中のデータ個数 n に支配される．今回の数値
実験では，標本の大きさを n = 10, 20, 50, 100, 200, 500,0および 1000 と 7 通りに変化させた．平均発生率は = 2

5
10
Wave height, H (m)15(a) GEV 分布

個 / 年に設定したので，高波発生の期間が K = 5, 10, 25,
50, 100, 250, および 500 年に相当する．確率波高を算定10する再現期間としては，R = 1, 2, 5, 10, 20, 50, 100, 200,

1

GPA Distribuition

500, および 1000 年の 10 通りについて検討した．非超10
Probability density, f (H)

過確率では PR = 0.5, 0.75, 0.90, 0.95, 0.975, 0.99, 0.995,
0.9975, 0.999, および 0.9995 にそれぞれ相当する．
統計解析にあたっては，標本の大きさ毎に N = 50,000個の標本を作製し，解析結果の信頼度を高めるようにした．
(4) 計算に使用する尺度・位置母数の値

10

-1

10

-2

10

-3

10

-4

既に (1) で述べたように，極値分布関数が与えられ，裾長度パラメータが設定されたときには，尺度母数 A

0

k = -0.2
k = -0.1
k = 0.2
k = 0.40と位置母数 B が式 (39) によって一義的に決定される．

5
10
Wave height, H (m)15(b) GPA 分布

今回の数値計算では H50 = 10.0 m とし，裾長度を 3 通りに変えているので，この条件の下で表-1 の母分布関数
10の尺度・位置母数の値を式 (39) で計算したところ，表1

Weibull Distribution-2～4 の母数値が得られた．
10
Probability density, f (H)

表で明らかなように，どの分布においても形状母数が大きくなるにつれて尺度母数と L-CV が増加し，位置母数と平均値が減少する．また，裾長度が大きくなると尺度母数と L-CV が増加し，位置母数と平均値が減少する．L-skewness と L-kurtosis は裾長度の影響を受けず，形状母数のみの関数である．
このようにして設定した極大値波高の確率密度関数

0

10

-1

10

-2

10

-3

10

-4

k = 0.75
k = 1.0
k = 1.4
k = 2.0をグラフ表示すると，裾長度が 50 = 1.20 の場合について図-2 のようになる．裾長度をそのように設定したこ

0とによって，波高 H = 8 ～11 m くらいの範囲では，分布関数による差異が小さい．分布関数の差異が大きくな

5
10
Wave height, H (m)(c) ワイブル分布

るのは，もっぱら波高の小さい範囲である．また，ワ

図-2 数値実験に用いる極値分布関数の確率密度関数

イブル分布と GPA 分布ではそれぞれの位置母数の値以（裾長度 50 = 1.20 の場合）

下の波高には確率密度関数が存在しない．

615

表-2 GEV 分布の母数値，平均値，および L-モーメント比 (H50 = 10.0 m, 平均発生率  = 2)
50 = 1.1

50 = 1.2

50 = 1.3

L-skw

L-krt

(3)

(4)

A (m)

B (m)

平均値

L-CV

A (m)

B (m)

平均値

L-CV

A (m)

B (m)

平均値

0.36020.26630.0999

8.458

8.595

0.0153

0.182

7.198

7.424

6.086

6.433

0. 0518

0.179

8.223

8.400

0.0213

0.325

6.775

0.454

5.488

5.939

0.0765-0.15

0.27000.1977

0.319

7.886

8.125

0.0319

0.584

6.124

7.067
6.562

0.03250.0464

0.254

-0.3

0.49400.37790.0723

0.809

4.633

5.240

0.125500.16930.1501

0.567

7.395

7.722

0.0509

1.026

5.289

5.823

0.1236

1.439

3.387

4.216

0.23640.15

0.07710.1209

0.997

6.686

7.131

0.0858

1.800

4.020

4.741

0.2367

2.530

1.589

2.718

0.5717

k
-0.45

L-CV

表-3 GPA 分布の母数値，平均値，および L-モーメント比 (H50 = 10.0 m, 平均発生率  = 2)
50 = 1.1

50 = 1.2

50 = 1.3

L-skw

L-krt

(3)

(4)

A (m)

B (m)

平均値

L-CV

A (m)

B (m)

平均値

L-CV

A (m)

B (m)

平均値

L-CV

-0.2

0.4286

0.2481

0.263

8.012

8.341

0.0219

0.482

6.355

6.958

0.0481

0.668

4.953

5.788

0.0801k-0.1

0.37930.2042

0.386

7.743

8.172

0.0276

0.707

5.863

6.649

0.0622

0.980

4.271

5.360

0.10690.001

0.33290.1663

0.567

7.395

7.961

0.0356

1.040

5.224

6.262

0.0829

1.439

3.387

4.825

0.14890.2

0.25000.1071

1.203

6.381

7.383

0.0617

2.205

3.364

5.202

0.1606

3.053

0.812

3.356

0.34460.4

0.17650.0642

2.539

4.659

6.472

0.1168

4.655

0.207

3.532

0.3922

4.885

0.342

3.831

0.37940.6

0.11110.0338

3.507

4.715

6.906

0.1221

6.429

0.310

4.328

0.3570––––

注：50 = 1.2 の k = 0.6 については平均発生率が  = 1 である．また，50 = 1.3 では k = 0.6 のケースを計算できない．

表-4 ワイブル分布の母数値，平均値，および L-モーメント比 (H50 = 10.0 m, 平均発生率  = 2)k

50 = 1.1

L-skw

L-krt

(3)

(4)

A (m)

B (m)

0.272

7.917

平均値

50 = 1.2

50 = 1.3

L-CV

A (m)

B (m)

平均値

6.180

6.774

L-CV

A (m)

B (m)

平均値

L-CV

4.711

0.75

0.45058.240

0.02375.533

0.08960.3333

0.565

7.399

7.964

0.0355

1.036

5.231

6.267

0.05290.08260.690

1.0

0.23650.1667

0.499

1.434

3.397

4.831

0.14841.4

0.21500.1227

1.155

6.562

7.615

0.0538

2.117

3.698

5.627

0.1339

2.932

1.274

3.945

0.26442.0

0.11400.1054

2.190

5.301

7.241

0.0785

4.015

1.385

4.943

0.2108

5.012

0.087

4.529

0.28722.8

0.04010.1040

3.701

3.614

6.910

0.1046

5.937

0.336

5.623

0.2062––––

注：50 = 1.2 の k = 2.8 および 50 = 1.3 の k = 2.0 については平均発生率が  = 1 である．また，50 = 1.3 では k = 2.8 のケースを計算できない．

(5) 位置母数の制約条件
極値分布関数のうちで k > 0 の場合の GPA 分布とワ

それを式 (34) ～ (36) に代入して計算する．すなわち，

イブル分布は下限値が B であり，波高は負の値をと

基準化変量 y50 と y10 は平均発生率  の関数であり，ることがないので，位置母数 B は正の値でなければ平均発生率が大きな場合には式 (40) の条件を満たさならない．GEV 分布で k < 0 の場合には下限値が B –ないことが起きる．GPA 分布とワイブル分布についA/k であるので，位置母数 B は A/k よりも大きくなけて，裾長度が 50 = 1.2 と 1.3 のときに許容される平均

ればならない．

発生率の最大値を試算すると，表-5 の結果が得られ

位置母数 B が正となる条件は，式 (39) の第二式をる．

書き換えることにより，次のように求められる．
1  1 / 50
B
 1
0
H 50
1  y10 / y 50

1  1 / 50
1
,
1  y10 / y50


表-5 位置母数を正値にとどめるための平均発生率の許容最大値 max

1  y10 / y50  1  1 / 50

ワイブル分布

GPA分布

50 =

50 =

50 =

50 =

1.2

1.3

1.2

1.3

72.1

7.74

25.9

5.59

0.1

k

y50
  50
y10

1.2

(40)
基準化変量 y50 と y10 は，再現期間をそれぞれ R = 50
および 10 年として式 (37) で非超過確率 PR を算定し，

7

k
0.1

GEV分布
k

50 =

50 =

1.2

1.3

26.1

5.72

1.4

25.3

3.75

0.2

7.57

2.59

0.15

13.0

3.75

1.6

11.5

2.17

0.3

3.53

1.55

0.2

7.68

2.69

2.0

3.87

1.03

0.4

2.09

1.07

0.3

3.63

1.63

2.8

1.03

< 1.0

0.6

1.04

<1.0

0.4

2.17

1.14

この表に示した平均発生率の許容最大値は，モンテカルロ法で数値シミュレーションを行う場合に考きの L-モーメントその他の偏りが無視できない大き慮しなければならない数値である．実際の高波の極さであったため，その改良として今回の数値実験の大値データを極値統計解析する場合には，極大値を過程で提案するものである．D は C よりも定数  の規定する波高の閾値によって高波の平均発生率が定値を少し大きくしてみたもので，ワイブル分布の k =まる．すなわち，閾値を低くするほど平均発生率が1.0 のケースに適用したところ結果が不具合であった

大きくなる．高波のデータに対する極値分布がある

ため，他のケースでは試みていない．

程度予測されるときには，表-2～4 の値を参考として

E は，ガンベル分布に対して Gringorten (1963) が提波高の閾値をある程度高めに設定する必要がある．

案したもので，最小 2 乗法で当てはめを行うときには偏りが生じないことを合田 (1988) が数値実験で確認している．F は，正規分布および対数正規分布に(6) 検討対象のプロッティングポジション式

対して不偏性を保持するプロッティングポジション

式 (1) のプロッティングポジションにおいては，として Blom (1965) が提案したもので，不偏性につい定数  と  にどのような値を与えるかによって極値ては同じく合田 (1988) が確認している．

データの L-モーメントの値が異なり，その結果としE と F のプロッティングポジション式の定数は次のて確率再現波高も異なるものとなる．

条件を満たしている．

極値統計において標本に分布関数を当てはめる方

2 +  = 1

法の選択に際しては，不偏性 (unbiasedness) と有効性
(efficiency) の二つの条件について検討しなければな

このように設定すると，式 (1) の昇順の序数 i を降順

らない．不偏性というのは，同一の母集団に属するの序数 m に替えても定数  と  として同じ値を使う

多数の標本について推計した分布関数の母数や確率ことができる．

再現統計量の平均値が母集団の値からずれていない

プロッティングポジション式としては，極値 II 型という条件である．また，それぞれの標本から推定分布を対象として提案された合田・小野澤 (1990) のした統計量は平均値の周りに分散するけれども，そのばらつきの度合いが小さいことを有効性と称する．
プロッティングポジションの与え方によって，不式その他がある．ただし，この両者とも形状母数をが不確定のときには適用がむずかしい．そのため，れるため，本研究では表-5 に示す 5 種類のプロッ

本研究では表-5 の 5 種類のプロッティングポジショ

ティングポジション式を取り上げた．

ン式を比較検討の対象としたものである．

A は，Hosking・Wallis (1995) が Wakeby 分布 (5 母数の分布) と GEV 分布について数値実験を行って B(7) モンテカルロ法による標本の作成法よりも優れていると結論付け，地域頻度解析の解説書で推奨しているものである．

モンテカルロ法では，[0, 1] の範囲に一様に分布する疑似乱数をコンピュータで発生させ，その値を非表-5 検討したプロッティングポジション式
式，ワイブル分布に対して提案された合田 (1988) のあらかじめ仮定することが必要なため，母分布関数

偏性と有効性の条件を満たす様態が異なると考えら



(41)

超過発生確率 Pi として式 (17), (23), (29) に代入して，

備考それぞれ GEV, GPA, ワイブル分布の確率統計量 xi を

A

0

0

Unbiased position

計算する．この操作を繰り返すことによって，所定

B

0.35

0

Landwehr et al.

C

0.45

0

New proposal

D

0.50

0

Temporal trial

似乱数の発生には合田 (2008, p. 294) のアルゴリズム

E

0.44

0.12

Gringortenを使用した．

F

0.375

0.25

Blom

記号

の大きさの極値の標本が作成できる．なお，一様疑発生させた標本が所定の分布に従っていることは，それぞれの標本について平均値と分散を求め，標本

B は Landwehr ほ か (1979) が 提 案した もの で ，の大きさごとに 50,000 回繰り返した結果の総平均とHosking・Wallis (1995, 1997) はこれを Plotting Position

総分散が，統計的変動性の範囲内で理論値に一致す式と呼んでいる．しかし，プロッティングポジショ

ることで確認した．

ンとしては他にもいろいろあるので，これだけをプ
ロッティングポジション式と呼ぶのは適切ではない．
C は，B のプロッティングポジションを用いたと8

今回の対象とした 3 種類の極値分布の平均値と分散は式 (16), (22), (28) で与えられる．このうち，平均

(1) 母集団の L-モーメント比値については表-2 ～4 に記載済みである．なお，各標本の平均値は，これらの表に記載した平均値の周りに  / n の標準偏差を持って変動する．たとえば，
ワイブル分布で裾長度が 50 = 1.2 であって k = 1.0, n =本研究で母集団として用いている 3 種類の分布関数について表-1 の形状母数を使って計算した L-モー
メント比の理論値は，表-2 ～4 に記載したとおりで

10 の標本であれば，式 (28) による標準偏差が  =ある．

1.036 m であるので，標本の平均値 6.27 m に対してそL- モ ー メ ン ト 比 の う ち で ， L-skewness (3) と L-の 90% 信頼区間が 5.73～6.80 m と推定される．

kurtosis (4) は式 (18), (24), (30) で明らかなように形状(8) 推計値の偏り率 (Bias) の定義

しかし，L-CV () は 2 次と 1 次の L-モーメント比であ

母数 k のみの関数であり，裾長度の影響を受けない．
統計解析で偏りというのは，一つの標本から推定るため，尺度母数と位置母数の値によって変化する．

された値が真値から離れている度合いをいう．ただこの辺りの状況は，表-2 ～4 の数値を比較すれば明し，個々の標本は統計的に変動するため，多数の標らかである．

本からそれぞれ得られた推計値の平均を用いて真値(2) 標本の大きさによる偏り率の変化からの偏りを判断する．
本研究では L-モーメント諸量，分布関数の母数値，
確率波高などを推計するので，こうした標本ごとの推計値を i，母集団の真値を 0 で表す．推計値の平均 [] ならびに標準偏差 [] は次のように定義される．

[ ] 1 N
i
N i 11 N [ ]     i2   2 [ ]
N
i



(42)
1/ 2

(43)

ここに，N は標本の総個数であり，本研究ではすべて N = 50,000 である．

標本から得られる各種の推計値は，標本が小さいほど偏りが大きく，変動も激しくなる．図-3 は，ワイブル分布で形状母数が k = 1.0 の場合，すなわち指数分布について L-CV, L-skewness, および L-kurotosisの偏り率（百分率表示）を示したものである．裾長度は 50 = 1.20 の場合である．いずれも表-5 に記した
6 通りのプロッティングポジション式を用いた結果を表示している．
3 種類の L-モーメント比のいずれについても，Unbiased Position ( = 0,  = 0) は偏り率が極めて小さく，その名称にふさわしい結果である．
Landwehr ほか (1979) によるプロッティングポジ
ション ( = 0.35,  = 0) は，L-CV と L-kurotosis に関す偏りの度合いは絶対値ならびに相対比の二つの方る偏り率が大きく，これを用いると不適切な推定結法で定義することができる．すなわち，Absolute Bias  [ ]   0(44)Relative Bias  ( [ ]   0 ) /  0(45)

果をもたらす危険性がある．これに対して新しく提案した定数の式 ( = 0.45,  = 0) は，偏り率を大きく縮減しており，Unbiased Position による結果にかなり

式 (44) の偏りはここでは偏り量と呼び，これは，分近づいている．試みにこの定数の値を  = 0.45 から 

布関数の形状母数などその値が 0 となる可能性のあ= 0.50 に変えてみたところ，L-CV と L-kurotosis に関るものに対して用いる．式 (45) の偏りの相対比はこする偏り率が負の値に大きく振れてしまい，プロッ

こでは偏り率と呼ぶことにする．なお，表示に際しティングポジションとして不具合な結果であった．

ては実数値ではなく，百分率で表すこともある．
れる標準偏差で表されるが，推計値の偏りも含めた = 0.50 のプロッティングポジションはこの分布関数以外には試してみなかったけれども， = 0.35, 0.45,
および 0.50 の比較結果から見て  = 0.45 のプロッ

2 乗平均誤差 RMSE で表されることも多い．これはティングポジションが実用的であると判断した．

推計値のばらつきの度合いは，式 (43) で算定さ

L-skewness の偏り率に関しては，定数の値を  =次のように計算される．





RMSE  [ ]   0 2   2 [ ]1/ 2

0.35 から  = 0.50 に変えてもほとんど差異が見られず，(46)

 = 0 の Unbiased Position よりはやや劣る結果であった．
一方，最小 2 乗法で使われる Gringorten 式と Blom

５．L-モーメント比の推計値の偏り率

式は，L-CV, L-skewness, および L-kurotosis のいずれ

9

以外のプロッティングポジション式は L-モーメント30

= 0= 0
= 0.35, = 0
= 0.45,= 0
 0.50, = 0
= 0.44, = 0.12
= 0.375= 0.25

Bias rate of L-CV, t (%)

Weibull Population
with k = 1.02010

比をはじめとする各種の推計量の偏りが大きいことが確認された．このため，L-モーメント法による極値解析ではこれらのプロッティングポジション式を推奨することができない．したがって以下においては，プロッティングポジション式を  = 0 と  = 0.450の二つに限定して記述する．
-10(3) 正規化された L-モーメント比の偏り率
-20
10

100

1000

図-3 は L-モーメント比の偏り率を標本の大きさにSample size, n

対して表示したものであり，これをみると偏り率が(a) L-CV

ほぼ標本の大きさ n に反比例している．この関係は

Bias rate of L-skewness, t  (%)

0

理論的にも予測されるもののようで， Hosking ほか
Weibull Population
with k = 1.0-50(1985) は一般化極値分布 (GEV) の偏り量に標本の大きさ n を乗じた量についての漸近理論値を示してい= 0, = 0
= 0.35, = 0
= 0.45,  = 0
= 0.50,= 0
= 0.44, = 0.12
= 0.375,= 0.25

-100

る．
そこで，これ以降においては偏り量あるいは偏り率に n を乗じた値を標本の大きさ毎に求め，それを複数の標本の大きさについて平均したものを分布関-150

数の代表値として用いる．ここでは，これを正規化された偏り量・偏り率と称する．平均値を求める際

-200
10

100
Sample size, n1000には，一般に標本の大きさが n = 50 ～500 のものを対象としたが，偏り率と標本の大きさの反比例関係が(b) L-skewness

確認できる場合には n = 20 ～1000 の範囲にまで拡げたケースもある．

Weibull Population
with k = 1.0100

図-4 は，一般化極値分布 (GEV)，一般化パレート
分布 (GPA)，およびワイブル分布のそれぞれについBias rate of L-kurtosis, t  (%)

150

50

て， L-CV, L-skewness, および L-kurtosis の正規化された偏り率と母集団の L-skewness [3] との関係を示し0たものである．プロッティングポジションとしては，

= 0,  = 0
= 0.35,  = 0
= 0.45, = 0
= 0.50,  = 0
= 0.44, = 0.12
= 0.375, = 0.25

-50
-100

 = 0 と 0.45 の 2 種類を用いた結果を区別して示している．前者の計算結果は黒塗りの記号，後者は白抜きの記号で表示している．

-150
10

100
Sample size, n

1000

裾長度を 3 通りに変えて数値計算した結果のうち，

L-skewness は裾長度の影響がほとんど見られないの(c) L-kurtosis

に対し，L-CV と L-kurtosis は裾長度によって異なる

図-3 ワイブル分布 (k = 1.0) に対する L-モーメント比の偏り率値を示すこともあった．ここでは，全体の傾向を示すために，3 通りの裾長度の計算結果の算術平均のについても偏り率が負の大きな値となり，L-モーメ

値で表示している．

ントの計算には不向きであることが例示される．

図-4(a) は L-CV の偏り率を示している．  = 0 のプ

数値シミュレーションでは，裾長度が 50 = 1.20 のロッティングポジションでは偏り率がゼロといって場合について A ～F のプロッティングポジションのよい． = 0.45 では偏り率がやや大きく，特に一般化

うち，D を除くすべてを使って標本を作成して統計極値分布では L-skewness が 0.3 以上のときに偏り率量の解析を行った．しかし，図-3 の例に見られるよが 1.0/n を超えている．

うに， = 0 の Unbiased Position と，新提案の  = 0.45

10Normalized CV of L-CV, n 0.5 x (t) (t)

4.0

2.0

Normalized Bias of L-CV (t )GEV ( = 0)
GEV ( = 0.45)
GPA ( = 0)
GPA ( = 0.45)
WBL (  = 0)
WBL (  = 0.45)

n x Bias of L-CV (t )

1.5
1.0
0.5
0
-0.5

CV of L-CV (t )3.5

GEV (empirical)
GEV ( = 0)
GEV ( = 0.45)
GPA (empirical)
GPA ( = 0)
GPA ( = 0.45)
WBL (empirical)
WBL (  = 0)
WBL (  = 0.45)

3.0
2.5
2.0
1.5
1.0
0.5
0

-1.0
0

0.1

0.2

0.3

L-skewness,

0.4

0

0.5

0.1

3

n x Bias of L-skewness (t3)Normalized CV of L-skewness, n 0.5 x (t3) (t3)4Normalized Bias of L-skewness (t3 )

0
-2
-4
-6
GEV ( = 0)
GEV ( = 0.45)
GPA ( = 0)
GPA ( = 0.45)
WBL ( = 0)
WBL ( = 0.45)

-8
-10
-12

0.5

12

CV of L-skewness (t3 )
10
GEV (empirical)
GEV ( = 0)
GEV ( = 0.45)
GPA (empirical)
GPA ( = 0)
GPA ( = 0.45)
WBL (empirical)
WBL (  = 0)
WBL (  = 0.45)

8

6

4

2

0
0

0.1

0.2
0.3
L-skewness,  3

0.4

0

0.5

0.1



0.2

0.3

L-skewness,

3

0.4

0.5



(b) n1/2 × (t3)(t3)(b) n × L-skewness (t3)10Normalized CV of L-kurtosis n 0.5 x (t4) (t4)8Normalized Bias of
L-kurtosis (t4 )4n x Bias of L-kurtosis (t4)0.4(a) n1/2 ×(t)(t)(a) n × L-CV (t)

2

0.2
0.3
L-skewness,  3




0

-4

GEV ( = 0)
GEV ( = 0.45)
GPA ( = 0)
GPA ( = 0.45)
WBL ( = 0)
WBL ( = 0.45)-8

CV of L-kurtosis (t4 )GEV (empirical)
GEV ( = 0)
GEV ( = 0.45)
GPA (empirical)
GPA ( = 0)
GPA ( = 0.45)
WBL (empirical)
WBL (  = 0)
WBL (  = 0.45)

8

6

4

2

0
0

0.1

0.2
0.3
L-skewness,  3

0.4

0

0.5



0.1

0.2
0.3
L-skewness,  3

0.4

(c) n1/2 ×(t4)(t4)(c) n × L-kurtosis (t4)図-4 正規化された L-モーメント比 の偏り量図-5 正規化された L-モーメント比の変動係数

11

0.5

L-skewness に関しては図-4(b) に示しており，偏りにも依存するので，その値は形状母数のみでなく裾率は L-CV よりも大きく，負の値を取るのが大半で長度によって変化する．しかし，標準偏差と平均値ある．二つのプロッティングポジションを比較するとの比，すなわち変動係数として整理すると，裾長と， = 0 のケースが全体に小さく， = 0.45 では偏度に依存しないことが数値計算の結果として確認さり率が大きくなる．ワイブル分布の  = 0 のケース以れた．そこで，標準偏差に n 1/2 を乗じたものを平均

外は，L-skewness が大きくなる（確率密度分布の右値で除した値を求め，これを正規化された変動係数

裾が長くなる）ほど偏り率が大きくなっている．まと呼ぶ．

た， = 0.45 のプロッティングポジションは  = 0 に正規化された変動係数に関しては， = 0 と = 0.45

よる偏り率の 2 倍以上の場合がある．

の二つのプロッティングポジションによる結果にほ

図-4(c) は L-kurtosis の偏り率を示したものである．

とんど差がなかったので，分布関数毎に両者の平均

これに関しては，二つのプロッティングポジションとして整理した．一般化極値分布 (GEV)，一般化パによる結果には有意な差がなく，ほぼ同じ偏り率を

レート分布，およびワイブル分布のそれぞれについ示 す ． 分 布 関 数 の 中 で は ， 一 般 化 極 値 分 布 で L-て整理した結果を図-5 に示す．

skewness が t3 = 0.494 (k = -0.45) のときには約 8 /n の大3 通りの裾長度による計算結果による変動係数の値は，L-CV と L-skewness に関してはほとんど差異がきな偏り率を示している．

なく，最大値と最小値の差は多くは 1% 以下，最大でも 5% 未満であった．ただし，L-kurtosis について

(4) 正規化された L-モーメント比の変動係数

は裾長度によって数 % の差異が見られることが多く，

L-モーメント比の推計値に関しては，平均値の偏最大では 20% 近い差が生じた．しかしデータ整理とりのみならず，推計値の変動の幅も検討する必要がしては，3 通りの裾長度による計算結果の平均値をある．先に Hosking ほか (1985) は，一般化極値分布用いることとした．

(GEV) の母数推計値の分散についても，それが標本

このようにして求めた  = 0 のプロッティングポジ


の大きさに反比例して減少することを前提にした漸ションによる結果は黒塗りの記号， = 0.45 のを用い

近理論値の図表を提示している．標準偏差は標本のた結果は白抜きの記号で示してある．図-5(a), (b), (c)

大きさの平方根，すなわち n 1/2 に反比例すると考えは，それぞれ L-CV, L-skewness, および L-kurtosis のられる．そこで，L-モーメント比の推計値の標準偏差（絶対値）については，それと n 1/2 の積を計算し，それを標本の大きさのある範囲（n = 50 ～1000）につ変動係数を母集団の L-skewness (3) に対してプロットしたものである．二つのプロッティングポジショ
ンによる計算結果にほとんど差が見られないところ

いて平均して代表値とすることとした．

から，両者の平均値に対して 2 次曲線を当てはめた．

L- モーメント比のうちで L-skewness と L-kurtosis

それによって得られた算定式は表-6 に示すとおりでは形状母数 k のみの関数であるので，裾長度には依ある．図-5 に描かれている曲線は，この近似推定式

存しない．しかし，L-CV は尺度母数と位置母数の値によるものである．

表-6 L-モーメント比の変動係数の推定式
L-モーメント比
L-CV (t)L-skewness (t3)L-kurtosis (t4)

変動係数 CV の推定式

分布関数
一般化極値分布 (GEV)

CV (t) = (21.6745 t32 – 6.7226 t3 + 1.256) /n0.5

一般化パレート分布 (GPA)

CV (t) = (7.2791 t32 – 0.4957 t3 + 0.4973) /n0.5

ワイブル分布 (WBL)

CV (t) = (6.4538 t32 – 1.4805 t3 + 0.7362) /n0.5

一般化極値分布 (GEV)

CV (t3) = (40.1193 t32 – 29.8590 t3 + 7.3288) /n0.5

一般化パレート分布 (GPA)

CV (t3) = (35.5486 t32 – 24.7527 t3 + 5.6864) /n0.5

ワイブル分布 (WBL)

CV (t3) = (28.5170 t32 – 23.0420 t3 + 5.8494) /n0.5

一般化極値分布 (GEV)

CV (t4) = (5.8445 t32 – 3.5959 t3 + 3.4001) /n0.5

一般化パレート分布 (GPA)

CV (t4) = (93.1882 t32 – 68.4313 t3 + 15.0430) /n0.5

ワイブル分布 (WBL)

CV (t4) = (-6.0637t32 + 2.2134 t3 + 2.6893) /n0.5

12表-7 形状母数の推計値の標準偏差 [ (k)] の変動係数の推定式


分布関数

ワイブル分布 (WBL)

ここで， L-skewness の変動係数のうちのワイブ
ル分布に関しては， k = 2.8 のケースが k = 2.0 以下のケースからかけ離れて大きいな変動係数を示したため，これを除外して2 次曲線の当てはめを行っている．
なお，図-4, 5 においては横軸の L-skewness が3 =0.3333 の 値 を 取 る 縦 軸 上 で に お い て ， 一 般 化 パ
レート分布 (GPA) とワイブル分布 (WBL) の線な
らびに計算点が交差している．ここではGPA が k =k

一般化パレート分布 (GPA)

0.25

Estimated value of shape parameter,

一般化極値分布 (GEV)

標準偏差の推定式

 (k) = (2.4457 t32 – 0.0543 t3 + 0.6998)/n0.5
 (k) = (7.1472 t32 – 3.9073 t3 + 1.6025)/n0.5
 (k) = (26.8575 t32 – 22.8583 t3 + 5.8874)/n0.5

0.20

GEV Population with k = 0.001
plotting by  = 0.45

0.15
0.10
0.05
0
-0.05
-0.10
-0.15
10

0.001 の値を取り， WBL が k = 1.0 の値であって，ともに指数分布を代表している．図-5(c) で両者が100
Sample size, N

1000(a) 一般化極値分布 (GEV)：k = 0.001
Estimated value of shape parameter,k

完全に一致しないのは，標本の統計的変動性が強く現れたためと思われる．また，推定式がこの点だけでなく，他の形状母数の点も含めた当てはめによって設定されていることも関係している．
表-6 に記載した推定式を使うと，図-5 に示した

L-skewness (3) の値の範囲で，標本から推計された
L- モーメント比の変動係数，ひいては L- モーメン
ト比の信頼区間を推計できる．たとえば，大きさ

n = 500 の標本から L-モーメント比として t = 0.1096,

0.6

GPA Population with k = 0.2
plotting by  = 0.45

0.5
0.4
0.3
0.2
0.1
0
10

t3 = 0.2589, t4 = 0.1299 の値が得られたとする．この

100

1000

Sample size, N

標本にワイブル分布を当てはめると，その形状母(b) 一般化パレート分布 (GPA)：k = 0.2れる（式 (31) による）．一方，表-6 中の近似式を用いて L- モーメント比の変動係数を計算し，それにそれぞれのて L- モーメント比の値を乗じると，標準偏差が次のように推定される．

 (t) = 0.0039,  (t3) = 0.0208,  (t4) = 0.0166.
推計された L- モーメント比の信頼区間を標準偏差の ±1.64 倍に取ると， L-CV, L-skewness, および L-

kurtosis についてそれぞれ次のように推定される．

Estimated value of shape parameter,k

数は L-skewness である t3 の値から k = 1.23 と推計さ

1.50

Weibull Population with k = 0.75
plotting by  = 0.45
1.25

1.00

0.75

0.50
10

(t)90% = 0.1033～0.1159,

100

1000

Sample size, N(t3)90% = 0.2248～0.2930,(c) ワイブル分布：k = 0.75(t4)90% = 0.1027～0.1571

図-6

標本の大きさによる形状母数の推計値の変化

６．分布関数の母数推計値の偏り量と偏り率計する．標本が小さいときには形状母数の推計値(1) 標本の大きさによる形状母数の推計値の変化

k̂ が真値 k0 から偏るだけでなく，推計値が標本毎に統計的に大きく変動する．

一つの標本に対して極値分布関数を当てはめる

図-6 は一般化極値分布 (GEV) ，一般化パレート

ときには 3(2) に述べたように，まず形状母数を推分布 (GPA) ，およびワイブル分布のそれぞれにつ13

いて形状母数の平均値とその 68% 信頼区間 ( 標準ことに関係している．すなわち，表-2 ～4 で明ら

偏差の ±1 倍の範囲 ) を示したものである．標本のかなように，形状母数と L–skewness は逆相関の関大きさが n = 10 のデータについては信頼区間の幅係にあるため，後者が小さくなると前者が大きくが大きすぎるため，表示から除外した．

なるのである．

3 種類の分布関数のうち，一般化極値分布と一般化パレート分布では L-skewness [3] の値の増加，(2) 正規化された形状母数の推計値の偏り量と標準偏差

すなわち形状母数の値の減少（正値から負値への変化）につれて形状母数の偏り量が増大する．ワ

図-6 から分かるように，標本の大きさが 50 未満

イブル分布では逆に，母集団の形状母数が k = 0.75

では信頼区間の幅が大きすぎ，形状母数の値を適から 2.8 へと増大するにつれて，形状母数推計値

切に推定することが困難である．また，形状母数の推計値と真値との差 k  kˆ – k0 は，標本の大きの偏り量が増加する．
形状母数の推計値の標準偏差については L- モー

さに反比例して減少しているように見える．そこ

メント比と同じように，標準偏差の値に標本の大で，これまでと同じように標本の大きさ n を乗じ

きさ n の 1/2 乗を乗じて正規化し，これを標本の大て正規化された偏り量をとりあげ，標本の大きさ

きさが n = 50 ～500 の範囲で平均化した．計算の結が n = 50 ～500 の範囲での平均値を求めた．その結果は図-8 のようになった．二つのプロッティング

果は図-7 に示すとおりである．この図には，形状

ポジションによる結果にはほとんど差がないので，

母数 k の偏り量についてプロッティングポジショ

その平均値を対象として 2 次曲線を当てはめ，形ンとして  = 0 と  = 0.45 の 2 種類を用いた結果を状母数の標準偏差の推定式を求めた．この推定式

併示しており， 3 種類の分布関数の形状母数推計

は表-7 に記載したようなものである．また，図-8

値の偏り量を一括して示した．

中の分布関数ごとの曲線は，この推定式で計算し形状母数の推計値に関しては，裾長度の影響をたものである．なお，ワイブル分布に対する当て

ほとんど受けないことが数値計算の結果で確認さはめに際しては， k = 2.8 のケースが他の形状母数れている．これは，図-3(b), 4(b) に示した L–skew-のケースと大きく外れる結果を示したので，これness が裾長度の影響を受けないことに対応していを除外して推定式を導いた．

る．

7

14

Normalized Bias of ShapeParameter, k
GEV ( = 0)
GEV ( = 0.45)
GPA ( = 0)
GPA ( = 0.45)
WBL (  = 0)
WBL ( = 0.45)

10

8

6

4

shape parameter, n 0.5 x (k)Normalized standard deviation ofn x Bias of shape parameter, k

12Standard Deviation of
Shape Parameter,  (k)6

GEV (empirical)
GEV ( = 0)
GEV ( = 0.45)
GPA (empirical)
GPA ( = 0)
GPA ( = 0.45)
WBL (empirical)
WBL ( = 0)
WBL ( = 0.45)

5

4

3

2

1

2
0

0

0

0

0.1

0.2
0.3
L-skewness,  3

0.4

0.5

0.1

0.2
0.3
L-skewness,  3

0.4

0.5





図-7

図-8

正規化された形状母数推計値の偏り量正規化された形状母数推計値の標準偏差

図-8 においては，一般化極値分布 (GEV) と一般

形状母数の偏り量 k はすべて正の値である．す化パレート分布 (GPA) は形状母数推計値の標準偏

なわち，形状母数の推計値は母集団の値よりも常差が (0.7 ～ 1.6) /n0.5 程度であって， L-skewness の値に大きく見積もられる．この傾向は図-6 でも例示にかかわらず，ほぼ一定である．これに対してワされている．偏り量が正値であるのは，図-4(b) にイブル分布では形状母数推計値の標準偏差が (1 ～おいて L–skewness の偏り量の大半が負の値である

14

6) /n0.5 の範囲で大きく変化している．ただし， L-

しては負値であり， (-0.1 ～ -13) /n の範囲にある．

skewness の値の減少に伴って形状母数の値が k =標準偏差は尺度・位置母数ともに L–skewness の値

0.75 から 2.8 に増大するので，変動係数としてみれの増加につれて減少する傾向を示しており，変化ば変化の幅は 2 倍以下である．

の範囲は尺度母数では (0.5 ～7.6) /n0.5 で，位置母数
では (0.1 ～7.4) /n0.5 である．裾長度の影響は一様でなお，形状母数推計値の標準偏差を指数分布にほぼ相当する k = 0.001 の一般化パレート分布と k =はなく，規定しがたい．

1.0 のワイブル分布とで比較すると，後者のほうがこのように，尺度母数と位置母数の推計値につ

約 24% 大きい．本来ならば同一の値になるべきもいては数値実験によってあらかじめその統計的挙のが食い違う理由については，今のところ不明で動を見定めることができない．与えられた標本にある．

対して推計された母数値を使い，モンテカルロ法
による数値実験で変動範囲を調べるのが唯一の方法と思われる．

(3) 尺度・位置母数の推計値の偏り量と標準偏差
既に表-2 ～4 に示したように，尺度母数と位置

７．確率波高の推計値の偏り率

母数の値は裾長度によって大きく変化する．このため，尺度・位置母数の推計値に関しては真値か(1) 当てはめ分布関数を母集団分布に固定した場らの偏りについては位置母数のような一般化が困

合の確率波高の偏り率難であり，分布関数ごとに定性的な記述にとどま標本に対して分布関数を当てはめ，形状・尺らざるを得ない．

度・位置母数の値を推計すると，各種の再現期間

一般化極値分布では，尺度母数の推計値 Â の偏り量は (-2 ～ 2) /n の範囲で L–skewness の増加につに対する確率波高を推定することができる．一般には高波極値の母集団が分からないため，複数のれて負から正の値に変化する傾向があり，一方，分布関数をあてはめの候補とし，最適合の分布を

裾長度が大きくなると偏り量が小さくなる傾向がもって母集団分布とみなしている．本研究では，ある．位置母数の推計値 B̂ の偏り量は (-0.4 ～ 1.0)
/n の範囲にあって L–skewness の値にあまり影響さ

あらかじめ設定した母集団から抽出した標本についてその統計的性質を調べているので，ここではれないが，3 通りの裾長度のうちでは 50 = 1.2 のと母集団の分布関数に当てはめて確率波高を推計しきの偏り量が小さい．標準偏差は，尺度・位置母た場合の偏り率について記述する．

数とも L–skewness の値が大きくなると減少する傾確率波高の偏り率は，分布関数の形状母数，裾向があるけれども，それらの値は裾長度によって

長度，標本の大きさ，プロッティングポジション

異なる．変化の範囲は (0.1 ～2.8) /n0.5 である．

式，および再現期間によって異なる．裾長度が大一般化パレート分布では，尺度母数の推計値のきな場合には，確率密度分布の右裾が長く伸び，偏り量が (0.5 ～ 8.7) /n の範囲にあるが，その変化を裾長度や L–skewness に関連づけることが難しい．
位置母数の推計値はおおむね負の値であり，最小再現期間の長い確率波高として大きな値が出やす
くなる．このため，裾長度の増大に伴って確率波
高の偏り率が大きくなる．

値は約 –2.3 /n である．尺度母数の推計値の標準偏標本の大きさの影響は， L- モーメント比の場合

差は， L–skewness の値の増加につれて減少する傾と同様に偏り率が標本の大きさ n にほぼ反比例す向があり，範囲としては (0.5 ～7.4) /n0.5 である．位ると見なすことができる．そこで，プロッティン

置母数の推計値の標準偏差も L–skewness の値の増

グポジション式毎に偏り率と標本の大きさの積の加 に つ れ て 減 少 す る 傾 向 が あ り ， 範 囲 と し ては平均値を求めて整理した結果が図-9 である．ただ(0.1 ～ 1.6) /n0.5 である．尺度・位置母数の推計値のし， n の小さい範囲や大きい範囲では偏り率と n -1

標準偏差は，裾長度が大きくなるにつれて増大すの比例関係が崩れる場合があり，特に再現期間のる傾向が認められる．

長いときに比例関係が崩れやすい．そのため，こワイブル分布では，尺度母数の推計値の偏り量

こでは N = 50～500 の範囲での平均値を示している．

が L–skewness の値の増加につれて減少する傾向が

なお，ここでは裾長度が 50 = 1.2 の結果を示しておあり，変化の範囲は (0.2 ～13) /n である．位置母数

り，裾長度による偏り率の変化については後述すの推計値の偏り量も同じ傾向を示すものの，値とる．

15プロッティングポジション式としては， = 0 と

4
n x Bias rate of return wave height0.45 の 2 種類を用いたものを併示している．ここに示しているのは式 (45) で定義される相対比としての偏り率である．たとえば，再現期間 R = 100 年に対する図の縦軸の値が 0.5 であったとし，確率波高が大きさ n = 100 の標本から推計されたものであれば，確率波高の推計値は 0.5% の偏り率を持つことになる．なお，本研究では高波の平均発生率を  = 2 に設定して標本を作製しているので，再現

k = -0.30 ( = 0)k = -0.15 ( = 0)
k = 0.001 ( = 0)
k = 0.15 ( = 0)k = -0.30 ( = 0.45)
k = -0.15 ( = 0.45)
k = 0.001 ( = 0.45)
k = 0.15 ( = 0.45)

3

2

GEV Population

1

0

-1

期間 R = 1000 年は非超過確率 PR = 0.9995 に相当す

1

る．

10

100

1000

Return period, R (year)

確率波高推計値の偏り率は 3 種類の分布関数の(a) 一般化極値分布 (GEV)

いずれにおいても，再現期間が 50 年 ( 非超過確率

PR = 0.99) 以下であれば ±0.5/n 程度以下であり，ほ5n x Bias rate of return wave height

ぼ不偏性を満足しているといえる．しかし，再現
期間が 50 年を超えるようになると偏り率が急速に増大し，再現期間 R = 500 年 (非超過確率 PR = 0.999)
では偏り率が 2/n 以上となるケースも見られる．
分布関数の中では，ワイブル分布の偏り率が最小であり，一般化パレート分布の偏り率が最も大きい．また， = 0 のプロッティングポジション式
は偏り率が正の値のことが多いのに比べて，  =0.45 のプロッティングポジション式は偏り率が負のことが多い．

k = -0.2 ( = 0)k = -0.1 ( = 0)
k = 0.2 ( = 0)
k = 0.4 ( = 0)k = -0.2 ( = 0.45)
k = -0.1 ( = 0.45)
k = 0.2 ( = 0.45)
k = 0.4 ( = 0.45)

4

3
2

GPA Population

1
0

-1

-2
1

偏り率の絶対値で比較すると，再現期間の長い範囲で = 0 のプロッティングポジション式のほう

10100Return period, R (year)

1000(b) 一般化パレート分布 (GPA)の偏り率が大きくなっている．これまでの L- モー
メント比や形状母数の推計値の場合には = 0 のほ

k = 0.75 ( = 0)
k = 1.0 ( = 0)
k = 1.4 ( = 0)k = 2.0 ( = 0)
k = 0.75 ( = 0.45)
k = 1.0 ( = 0.45)
k = 1.4 ( = 0.45)
k = 2.0 ( = 0.45)

2.0
n x Bias rate of return wave height

うが小さな偏りを示していたことを考えると，確率波高の推計値に関して二つのプロッティングポ
ジション式の優劣が逆転するのはやや不思議な気もする．
もっとも， Hosking ・ Wallis (1995) が  = 0 と  =0.35 のプロッティングポジションを比較検討した
ときには，地域頻度解析などで分布関数を適正に当てはめるためには  = 0 を使うべきであると結論する一方で，確率波高の推計が主目的であれば  =

1.5
1.0

WBL Population

0.5
0.0
-0.5
-1.0
-1.5
1

0.35 の使用も許容できるとしている．今回新しく

10

100

1000

Return period, R (year)

提案した = 0.45 のプロッティングポジションは (c) ワイブル分布= 0.35 の方式よりも推計値の偏りがはるかに小さ

図-9いので，これを確率波高の推定に使用するのは当正規化された確率波高推計値の偏り率（50 = 1.2 の場合）を得たものと思われる．
裾長度の影響を 3 種類の分布関数のそれぞれか

数分布である k = 1.0 である．図で明らかなように，

ら 1 ケースについて調べた結果が図-10 である．一

裾長度が大きくなるにつれて確率波高の偏り率の般化極値分布はガンベル分布に相当する k = 0.001 ，

絶対値が増大している．

一般化パレート分布は k = 0.2，ワイブル分布は指16その際には偏り率ができるだけ小さいことが望ま

n x Bias rate of return wave height

2.5

50 = 1.3 ( = 0)

2.0

50 = 1.3 ( = 0.45)

GEV Population(k = 0.001)

れる．その意味で，波高の極値統計解析ではプ
ロッティングポジションとして新しく提案する =50 = 1.2 (= 0)0.45 の定数を使用するのが適切と考えられる．

50 = 1.2 ( = 0.45)

1.5

50 = 1.1 ( = 0)

50 = 1.1 ( = 0.45)

1.0(2) 最適合分布関数が母関数に的中する割合

0.5

前項では当てはめ分布関数を母集団分布に固定した場合の確率波高の偏り率を記述した．実際の0極値統計解析では，標本がどのような母集団から

-0.5
1

10
100
Return period, R (year)1000

得られたものかが未知であり，複数の分布関数のうちで標本にもっとも適合すると考えられるもの(a) 一般化極値分布 (k = 0.001)

を母集団とみなしている．
今回の数値シミュレーションでは 一般化極値分

n x Bias rate of return wave height

4.0

50 = 1.3 ( = 0)

3.0

GPA Population(k = 0.2)

布，一般化パレート分布，およびワイブル分布が高波極値の母集団分布の候補であると考え，その50 = 1.3 ( = 0.45)
50 = 1.2 ( = 0)

2.0

中から最小残差法に基づく最適合の分布を当ては

50 = 1.2 ( = 0.45)

め候補として取り上げている．すなわち，3(3) で

50 = 1.1 ( = 0)

50 = 1.1 ( = 0.45)

1.0

述べたように，標本のデータを昇順に並べ替えた
値 xi とその推計値 x̂i の間の相関係数 r ( xi , xˆi ) を計0算し，その上限値である 1 からの残差r  1  r ( xi , xˆi ) が最小となるものを最適合と判断し

-1.0
1

10

100

1000

ている．
ここで，一つの母集団から得られた極大値の標Return period, R (year)(b) 一般化パレート分布 (k = 0.2)

n x Bias rate of return wave height2.0

50 = 1.3 ( = 0)
50 = 1.3 ( = 0.45)

1.5

50 = 1.2 ( = 0)本に対して， L- モーメント法によって最適分布と認定された分布関数が元の母集団分布に合致する

WBL Population(k = 1.0)

割合をここでは当てはめ的中率を呼ぶことにする．
この的中率は，標本が大きくなるにつれて向上す50 = 1.2 ( = 0.45)

ることが予想され，このことは今回の数値計算の50 = 1.1 ( = 0)

1.0

50 = 1.1 ( = 0.45)

結果でも確認されている．

0.5

当てはめ的中率が標本の大きさによって変化する状況を図-11 に示す．裾長度を 3 通りに変えて計0算した結果には有意な差が見られなかったので，-0.5

当てはめ的中率としては 3 通りの裾長度の結果の
-1.0
1

10

100

1000

Return period, R (year)と小さい場合には的中率が低くなるため，図には(c) ワイブル分布 (k = 1.0)
図-10

平均値を用いている．また，標本の大きさが n = 10
表示していない．

正規化された確率波高推計値の偏り率（裾長度の影響の比較）

当てはめ的中率は標本が大きくなるにつれて向上する．例外は一般化パレート分布で形状母数の値が負の場合であり，標本の大きさが n = 100 の辺

また，プロッティングポジションの影響をみる

りで的中率が 30% 前後に落ち込んでいる．この付と，  = 0.45 よりも  = 0 のほうが偏り率の絶対値

近では，一般化パレート分布よりもワイブル分布が大きいことが確認される．特に，裾長度が 50 =が最適合と判断される割合が多くなっている．図-1.3 の場合に差異が明瞭である．波高の極値統計の

1 の L-kurtosis と L-skewness の相関図に見られるよ

最終目的は，与えられた標本に基づいて所定の再うに，形状母数の値が負のときの一般化パレート

現期間における確率波高を推計することであり，

17

分布はワイブル分布と近接した位置を占めており，

100

標本の大きさが十分でないときには，母集団の識GEV - to - GEV return rate (%)

GEV Population

別がむずかしいのであろう．

k = -0.1580

また，一般化極値分布で形状母数が k = -0.45 の60

k = -0.15

ケースでは，プロッティングポジションとして  =k = -0.3

0 を用いたときには，標本が n = 1000 と非常に大きくなっても当てはめ的中率が 50% にとどまる．こk = -0.4540のときには一般化パレート分布を最適とする割合
が 46% に達している．ただし，プロッティングポ

20

ジションとして  = 0.45 を用いると，標本が大き

Closed symbols:  = 0
Open symbols:  = 0.45

くなるにつれて一般化極値分布を最適合とする割

0
10

100
Sample size, n合が増大し，n = 1000 では的中率が 74% となる．

1000

なお，図-11(b) の一般化パレート分布の k = 0.001のケースと図-11(c) のワイブル分布の k = 1.0 のケー(a) 一般化極値分布 (GEV)スは，ともに指数分布を代表している．両ケース
とも，標本が大きくなっても的中率が 50% にとど

100

GPA Populationk = 0.4

GPA-to-GPA return rate (%)

80

まるのは，この場合には二つの分布関数を識別することが理論的に不可能なためである．

k = 0.2

著者は以前に極値統計解析の信頼度手法について，最小 2 乗法を用いた場合の結果を検討した

60
k = 0.001
40（合田, 1988）．そのときには，ワイブル分布で形状母数を 4 通りに固定したものとガンベル分布を

k = -0.1

組み合わせた 5 種類の 2 母数型分布を当てはめ対k = -0.2

20

象として適合率をしらべた．今回のように 3 母数

Closed symbols:  = 0
Open ymbols:  = 0.45型分布の中からの選択でないため状況が異なるけ
れ ど も ， 標 本 の 大 き さ が 100 程 度 で は 的 中 率 が

0
10

100
Sample size, n

1000

50% ～ 80% であり，ガンベル分布の場合には的中率 53% であって，残りは k = 1.4 あるいは k = 2.0のワイブル分布のほうによく適合した．

(b) 一般化パレート分布 (GPA)

分布関数の当てはめ的中率は最適合の判定基準に依存するので，今回用いた最小残差法以外の基

100

Weibull-to-Weibull return rate, (%)

Weibull Population

準を使うと当てはめ的中率としてやや異なる数値
k = 2.0が得られるであろう．しかしながら，標本の大き80

さが 100 程度未満では，標本の背後に隠れている

k = 1.4

母分布関数を的確に見つけだすのはかなり難しい

60
k = 0.75

k = 2.8といえよう．

k = 1.040(6) 母集団が未知の場合の確率波高推計値の偏り率

20
Closed symbols:  = 0
Open symbols:  = 0.45

実際の極値統計解析では母集団が未知であり，

0
10

100
Sample size, n

複数の当てはめ候補の中から最適合の分布関数を1000選んだとしても，それが真の母集団を表しているという保証はない．

(c) ワイブル分布
図-11

今回の数値実験では，高波極値の母集団である母集団への当てはめ的中率と標本の大きさの可能性の高い 3 種類の分布関数，すなわち一般化

関係18極値分布，一般化パレート分布，およびワイブル

5

標本毎に最適合分布関数を選び，4(3) に述べた 10通りの再現期間について確率波高を推計し， N =50,000 回の繰り返し計算で得られた標本毎の確率
波高推計値の平均値ならびに標準偏差を算定した．
そして，推計値の平均について偏り率を計算している．なお，この際に 50,000 個の標本に当てはめられた分布関数の形状・尺度・位置母数はそれぞ

GEV Population with k = -0.3
plotting by  = 0.45

Bias rate of return wave height (%)

分布を母集団の代表として選んでいる．そして，

0-5
n = 20n = 50n = 100n = 200
n = 500
n =1000-10-15

れに異なる．

1

10
100
Return period, R (year)

これまでの検討結果から明らかなように，式 (1)のプロッティングポジション式において定数を  =1000(a) k = -0.3 のケース

0 と  = 0.45 に設定したときの推計値の信頼度には

4Bias rate of return wave height (%)

有意な差が見られない．しかし，  = 0 の使用には一つの問題が残る．すなわち，式 (1) に i = n を代入すると Pn = 1 となり，プロッティングポジショ
ンが定められない．このため，データ中の最大値を図面にプロットして当てはめた分布関数と比較することができない．また，最適合関数を定めるための相関係数 r ( xi , xˆi ) を i = 1 ～n に適用できない

GEV Population with k = -0.15
plotting by  = 0.45

2
0
-2
-4
n = 20n = 50n = 100n = 200
n = 500
n =1000

-6
-8
-10
-12

ため，相関係数の計算は i = 1 ～ (n - 1) に限定して

1

10
100
Return period, R (year)

最適合の判定を行なわざるを得ない．

1000(b) k = -0.15 のケース

このようなところから，高波の極値統計解析に
Bias rate of return wave height (%)

L- モーメント法を適用する際には，プロッティン
グポジションとして  = 0.45,  = 0 を使用するのが適切と考えられる．したがって，以下に述べる母集団未知の場合の確率波高の推計値の偏り率なら
びに推計値の変動係数の検討においては，すべて
このプロッティングポジション式を用いて計算した結果を使用する．
母集団が未知の場合の確率波高の推計値の偏り

4

GEV Population with k = 0.001
plotting by  = 0.45

2
0
-2
-4

n = 20n = 50n = 100n = 200
n = 500
n =1000

-6
-8
-10
1

率については，裾長度が 50 = 1.2 の場合の計算結果

10100
Return period, R (year)

1000

標本に当てはめられた最適合分布から推計された
確率波高の偏り率を百分率で表示したもので，形状母数の値毎に 4 枚のグラフとして，標本の大きさ別にプロットしてある．これまでの図表のように標本の大きさ n 乗じて正規化したものではない．
たとえば，一つの標本に対して形状母数が k = 0 の一般化極値分布が最適合と判定されたとき，標本

Bias rate of return wave height (%)(c) k = 0.001 のケース

を分布関数毎に図-12～14 に示す．
図-12 は一般化極値分布の母集団から抽出されたの大きさが n = 20 であれば 100 年確率波高の推計値

GEV Population with k = 0.15
plotting by  = 0.452
0-2
n = 20n = 50n = 100n = 200
n = 500
n =1000

-4
-6
-8
1

10

100

1000

Return period, R (year)(d) k = 0.15のケース
図-12が -6% の偏り率，標本の大きさが n = 100 であれば100 年確率波高の推計値が -2.3% の偏り率を持つと読み取られる．

19

母集団が未知の場合に一般化極値分布の標本への最適合分布から推計された確率波高の偏り率(50 = 1,2,  = 0.45)

15

GPA Population with k = -0.2
plotting by  = 0.45

Bias rate of return wave height (%)

Bias rate of return wave height (%)

50

N = 20N = 50N = 100
N = 200
N = 500
N =1000

-5

-10
1

10
100
Return period, R (year)

1000

Weibull Population with k = 0.75
plotting by  = 0.45
10

5

0

-5

(a) k = -0.2 のケース

1

10

GPA Population with k = -0.1
plotting by  = 0.45

3
2

1000

(a) k = 0.75 のケース

1
0n = 20n = 50n = 100n = 200
n = 500
n =1000

-1
-2
-3
-4

100

Return period, R (year)

4

Bias rate of return wave height (%)

Bias rate of return wave height (%)

n = 20n = 50n = 100n = 200
n = 500
n =1000

1

10

100

1000

Return period, R (year)(b) k = -0.1 のケース

10

Weibull Population with k = 1.0
plotting by  = 0.45

n = 20n = 50n = 100n = 200
n = 500
n =1000

5

0

-5
1

10
100
Return period, R (year)

1000

GPA Population with k = 0.2
plotting by  = 0.45

8(b) k = 1.0 のケース
Bias rate of return wave height (%)

Bias rate of return wave height (%)

10

n = 20n = 50n = 100n = 200
n = 500
n =1000

6
4
2
0
-2
-4
1

10
100
Return period, R (year)

1000

Weibull Population with k = 1.4
plotting by  = 0.45

42

n = 20n = 50n = 100n = 200
n = 500
n =1000

0
-2
-4
1

10
100
Return period, R (year)

1000(c) k = 0.2 のケース(c) k = 1.4 のケース
GPA Population with k = 0.4
plotting by  = 0.45
10Bias rate of return wave height (%)

Bias rate of return wave height (%)

15

n = 20n = 50n = 100n = 200
n = 500
n =1000

5

0

-5
1

10
100
Return period, R (year)

(d)
図-13

Weibull Population with k = 2.0
plotting by  = 0.452
0-2

n = 20n = 50n = 100n = 200
n = 500
n =1000

-4
-6
-8
1

1000

k = 0.4のケース

10
100
Return period, R (year)

1000(d) k = 2.0 のケース

母集団が未知の場合に一般化パレート分布の標本への最適合分布から推計された確率波高の偏り率 (50 = 1,2,  = 0.45)

図-14

20

母集団が未知の場合にワイブル分布の標本への最適合分布から推計された確率波高の偏り率(50 = 1,2,  = 0.45)一般化極値分布から抽出された標本の場合には，

2Bias rate of return wave height (%)

確率波高の偏り率は負の値であるのが大半である
が，一般化パレート分布から抽出された標本の場
合には図-13 に示すように偏り率が正の値のことが多い．
ワイブル分布から抽出された標本について，母集団が未知であるとして最適合の分布関数から確率波高を推計した場合の偏り率は，図-14 に示すとおりである．偏り率の絶対値は，一般化極値分布

GEV Population(k = 0.001)
0-2
n = 50 (50 = 1.3)
n = 50 (50 = 1.2)
n = 50 (50 = 1.1)

-4

n = 500 (50 = 1.3)
n = 500 (50 = 1.2)-6および一般化パレート分布の場合よりもやや小さ

n = 500 (50 = 1.1)
1めであるものの，標本の大きさならびに再現期間の長さによる確率波高推計値の偏り率変化の様相

10
100Return period, R (year)

1000(a) 一般化極値分布 (k = 0.001)

は類似している．

12
Bias rate of return wave height (%)

こうした母集団が未知の場合の確率波高の偏り
率は裾長度によっても影響される．図-15 は特定のケースについて裾長度の影響を示したもので， (a)は一般化極値分布の k = 0.001 のケース，(b) は一般
化パレート分布の k = 0.2 のケース，(c) はワイブル
分布の k = 1.0 のケースである．図の混乱を避けるため，標本の大きさが n = 50 および 500 の場合を示している．
図-15 に示した例では，裾長度が 50 = 1.1 から 1.3

n = 50 ( 50 = 1.3)

10

n = 50 ( 50 = 1.2)

8

n = 50 ( 50 = 1.1)

6

n = 500 (50 = 1.2)

GPA Population(k = 0.2)

n = 500 (50 = 1.3)
n = 500 (50 = 1.1)

4
2
0
-2

1

へと大きくなるにつれて，確率波高の偏り率が増大している．特に，標本の大きさが n = 50 のとき

10100
Return period, R (year)

1000(b) 一般化パレート分布 (k = 0.2)に裾長度の影響が顕著に認められる．ただし，図-10Bias rate of return wave height (%)15(a) の一般化極値分布の場合には裾長度が 50 =

1.2 のときよりも 50 = 1.3 のときの偏り率がやや小さめになっている．
極値統計解析の実務においては，今回の数値実験のように与えられた標本の母集団を知ることができない．したがって，厳密にいえば図-12～14 に提示した情報が使えるわけではない．また，裾長度の影響も定式化がむずかしい．しかし，波高のn = 50 (50 = 1.3)

8

n = 50 (50 = 1.2)

Weibull Population(k = 1.0)

n = 50 (50 = 1.1)

n = 500 (50 = 1.3)

6

n = 500 (50 = 1.2)
n = 500 (50 = 1.1)

4
2
0
-2

極値統計解析に限っていえば，高波の極値分布関1数の裾長度は平均的に 1.2 程度であるので，まずは図-12～14 の偏り率を参考にすることができよう．
関数に応じてこれらの図表から偏り率を読み取り，
確率波高の推計値を補正するのが良いと思われる．
なお，図-12～15 では横軸を再現期間としている

1000(c) ワイブル分布 (k = 1.0)
図-15

すなわち，標本に対して最適合と判定された分布10100Return period, R (year)

母集団が未知の場合の確率波高の偏り率に及ぼす裾長度の影響

らの図表の再現期間 1000 年に相当する．
図-12～15 での結果でも分かるように，確率波高が，これは平均発生率を  = 2 に設定したときの非の偏り率はかなり大きなものである．偏り率を下超過確率から換算したものである．非超過確率はげるためには，十分に大きな標本を使うことが不式 (37) で算定するので，たとえば平均発生率が 可欠である．少なくてもデータ数が 100 以上の標= 20 であれば，そのときの再現期間 100 年はこれ本を使うことが望ましいといえる．

21

８．確率波高の推計値の変動係数よる違いがほとんど認められなかった．このため，
以下においては式 (1) に  = 0.45,  = 0 を代入した新極値統計においては，推定された確率波高の絶提案のプロッティングポジション式を用いて得ら対値だけでなく，その信頼区間も必要な情報であれた結果について記述する．

る．今回の数値計算では 4(3) に述べたように，

式 (47) の右辺に導入されている標本の大きさの

10 通りの再現期間毎に 50,000 個の標本から求めた影響については，標準偏差あるいは 2 乗平均誤差

確率波高の平均値と標準偏差を算出している．こと同じように n にほぼ逆比例することが認められた．また，裾長度の影響に関しては大略の傾向の結果に基づいて，確率波高の変動係数を次のように定義した．

CV [ H R ] として (50 – 1)0.7 に比例することが認められた．そ [H R ] f (k , n,  50 , y R )[ H R ]こで，変動係数に対して次のような基準化変量の(47)

2 次式を当てはめることを試みた．

CV [ H R ]  (a y R2  b y R  c) ( 501) 0.7 / n

ここに，  [HR] と [HR]は再現期間 R 年の確率波高の標準偏差と平均値， yR は R 年に対する基準化変(48)

そして，母集団の分布関数の形状母数毎に定数 a,量である．式 (47) は，変動係数が形状母数 k ，標b, c を算定し，さらにその結果を形状母数の関数

本の大きさ n ，裾長度 50，および基準化変量 yR のとして表示することを試みた．これによって，分関数であることを表している．この変動係数の大布関数毎に表-8 のような経験式を得ることができ

きさについては，プロッテイングポジション式にた．

表-8 確率波高の変動係数の近似推定式の係数表
分布関数
一般化極値分布 (GEV)

一般化パレート分布 (GPA)

ワイブル分布係数の推定式
a = 2.7037 k3 + 2.4351 k2 + 0.7325 k + 0.0830
b = 10.953 k3 + 10.509 k2 + 3.1359 k + 0.0372
c = 11.311 k3 + 12.526 k2 + 5.6088 k + 0.7719
a = exp[-5.3844k3 – 0.6089k2 +7.7316k – 1.9805]
b = –22.438k3 + 24.678 k2 + 6.1046 k + 0.1062
c = –28.287k3 +18.860 k2 + 6.9797 k + 0.6776
a = 0.9742 (lnk)2 + 0.6949 (lnk)+ 0.1517
b = 1.9049 (lnk)2 + 1.8142 (lnk)+ 0.3677
c = 0.6349 (lnk)2 + 1.8561 (lnk) + 0.9095

これらの経験式による確率波高推計値の変動係数の予

次の図-16(b) の一般化パレート分布の場合には，予測測値を数値実験で得られたデータと比較した結果の例を

値は数値計算データをかなり下回っており，適合度はあ

図-16 に示す．ここでは，一般化極値分布 (k = 0.001)，一

まりよくない．この分布関数では，形状母数が大きく

般化パレート分布 (k = 0.2)，およびワイブル分布 (k = 1.0)

なったときに式 (48) のような 2 次曲線を当てはめると，のケースを示している．横軸は式 (34)～(36) で定義され基準化変量が小さくなるにつれて変動係数が増加する傾る基準化変量であり，連続した曲線が式 (48) に表-8 で計向が過大に表現されてしまう．当てはめる近似曲線を変算される係数値を代入して求めた予測値である．記号でえることによって適合度が向上するかもしれないけれど示しているデータは，標本の大きさが n = 20 ～1000 の 6も，3 種類の分布関数のすべてについて広範囲の形状母通りの計算結果であり，3 種類の裾長度ごとに記号を変数に共通して適用できるような近似式を見いだすことはえて示している．

むずかしい．このケースは，近似式の適用限界を示す資

まず，図-16(a) の一般化極値分布の場合には，数値実料として参照すべきであろう．それでも，数値計算デー

験によるデータは式 (48) の予測値の周辺に散らばってタは裾長度の増大に応じて変動係数も大きくなっており，おり，予測式はほぼ妥当である．なお，同一の基準化変

その傾向は予測値でも同様である．すなわち，式 (48) に量の値でデータが上下に散らばっているのは，変動係数

組み込まれている (50 – 1)0.7 の項は裾長度の効果を適切

が標本の大きさの –1/2 乗に比例するとの仮定からのずにぼ表現しているといえる．

れによるものである．

22

本の大きさによるばらつきが大きい．これは，標本の大2.5

GEV Population with k = 0.001
plotting by = 0.45

n 1/2 x CV of HR

2.0

きさの影響を n-1/2 で近似したことによる誤差である．

50 = 1.3 (est)

こに示した表-8 の係数を用いた式 (48) の推定式は，確50 = 1.3 (cal)

率波高推計値の信頼区間の概略値を求めるためのもので

50 = 1.2 (cal)あるので，図-16 に見られるような推定式とデータの差

50 = 1.2 (est)1.5

50 = 1.1 (est)1.0

異は，実用的には許容される範囲と言えよう．実際の極50 = 1.1 (cal)

値統計解析では，与えられた標本に対して最適合の分布関数を定め，所定の再現期間に対する確率波高を推計す

0.5

る．それに対する信頼区間はその概略値を式 (48) と表-8
0
0

1

2

3
4
5
Reduced variate, yR

67を用いて算定し，さらに詳しくは推計された分布関数の8母数を使ってモンテカルロ法による数値シミュレーショ

(a) 一般化極値分布 (k = 0.001)
2.5

50 = 1.3 (est)

ンを実施して信頼区間を見いだすことが望ましい．
ここで，確率波高の推計値の信頼区間の推定の試算例

GPA Population with k = 0.2
plotting by = 0.45を示す．いま，大きさが N = 120 で平均発生率が = 2 の極大値波高の標本について統計解析を行い，形状母数が

50 = 1.3 (cal)2.0
n 1/2 x CV of HR

50 = 1.2 (est)k = 1.2 のワイブル分布が最適合と判定されたとする．尺50 = 1.2 (cal)1.5度・位置母数はそれぞれ A = 1.69 m, B = 4.47 m であった．

50 = 1.1 (est)50 = 1.1 (cal)

この分布が母集団であるとみなして，再現期間 R = 100

1.0

年に対する確率波高とその信頼区間を試算してみる．
まず，この分布関数における確率波高と裾長度を求め0.50ておく必要がある．再現期間 10, 50, および 100 年に対する非超過確率は，式 (37) によってそれぞれ P10 = 1 –

0.0
0

1

2
Reduced variate, yR

3

4

1/(2×10) = 0.95, P50 = 1 – 1/(2×50) = 0.99, P100 = 1 – 1/(2×100) =0.995である．したがって，基準化変量は式 (36) で y10 =(b) 一般化パレート分布 (k = 0.2){-ln(1 – 0.95)}1/1.2 = 2.495, y50 = {-ln(1 – 0.99)}1/1.2 = 3.570, y100 =
3.5

3.0

れ H10 = 4.47 + 1.69×2.495 = 8.687m, H50 = 4.47 + 1.69×3.570= 10.503m, H100 = 4.47 + 1.69×4.013 = 11.25m である．これ
50 = 1.3 (est)によって，裾長度は 50 = H50/H10 = 10.503/8.687 = 1.209 と算2.0

50 = 1.2 (est)

定される．

1.5

50 = 1.1 (est)2.5
n 1/2 x CV of HR{-ln (1 – 0.995)}1/1.2 = 4.013 となるので，確率波高はそれぞ

Weibull Population with k = 1.0
plotting by = 0.45

50 = 1.3 (cal)50 = 1.2 (cal)

確率波高の推計値の変動係数を求めるための定数は，
表-8 中の推定式によって次のように計算される．

50 = 1.1 (cal)a  0.974 (ln1.2) 2  0.695 ln 1.2  0.152  0.311
2b  1.905 (ln 1.2)  1.814 ln 1.2  0.368  0.762
c  0.635 (ln 1.2) 2  1.856 ln 1.2  0.910  1.269


1.0

0.5
0
0

1

2

3
4
5
Reduced variate, yR

6

7

したがって，変動係数は次のように計算される．

8CV [ H R ]  (0.311 4.013 2  0.762  4.013  1.269)(c) ワイブル分布 (k = 1.0)

 (1.209  1.0) 0.7 / 120  0.0982

図-16 母集団が未知の場合に標本への最適合分布から推計された確率波高の変動係数

すなわち，100 年確率波高の標準偏差は 11.25×0.0982 =

1.105 m である．90% 信頼区間として標準偏差の 1.64 倍
最後の図-16(c) のワイブル分布の場合は，図-16(a)の一を用いると，この条件における 100 年確率波高 H100 =般化極値分布と同じ程度に数値実験データが予測曲線の11.25 m の信頼区間が 9.44～13.06 m と推定される．

周りに散らばっている．裾長度が 50 = 1.1 のデータは標23

９．L-モーメント法のプロッティングポジション等に関する既往研究との比較

で直接に比較できないけれども，同じオーダーの偏り率が得られている．）また，高水流量の頻度解析データの事例として地点数 21，地点毎のデータ数が 10 ～30 の(1) Unbiased Position について

ケースについて PP と unbiased の偏り率を比較している．

最初にまえがきで紹介したように，式 (1) のプロッ

その結果の表から計算すると，非超過確率 0.99 およびティングポジション式で定数として  = 0,  = 0 の組合せ

0.999 において PP の偏り率が 0.6% と 1.2%，unbiased によを推奨したのは Hosking・Wallis (1997) の解説書 (33-34頁)

る偏り率が 1.1% と 1.5% であり，PP の偏り率のほうが低である．その根拠としているのは Hosking・Wallis (1995)

くなっている．

の論文である．

こうした計算結果からHosking・Wallis (1995) は，単一

この論文では，Landwehr ほか (1979b) が用いた  =地点の標本について確率波高などを推定する際には PP

0.35,  = 0 の組合せを“plotting position estimator” (以下ではに若干の優位性があると認めている．今回新しく提案しPP と略記) と呼び， = 0,  = 0 の組合せを “unbiased

た = 0.45,  = 0 のプロッティングポジション式は 7(1) でestimator” と呼称している．後者は厳密には unbiased でな示したように，Unbiased Position よりも確率波高の推計

いけれども，L-モーメントの計算式の構成に基づいて値の偏り率が小さい．また，L-モーメント比や形状母数

unbiased として引用するのが便利であるとしている．

の推計においても，Unbiased Position よりも偏り率がや
や大きいものの，実用的には無視できる量である．

Hosking・Wallis (1995) が PP が適切でないとしている一つの理由は，PP を使うと標本のデータ座標軸を移動（各データにある一定値を加除）したときに t2, t3, t4 など(2) GEV 分布の偏り量の漸近理論値問題の L-モーメントおよび L-モーメント比の値が変化する，

L- モーメント法の母体である Probability Weighted

すなわち不変性が保持されないということである．しかし，データの座標軸を移動することの意義が不明であり，それに対する L-モーメントの不変性を保持すべき理由も

Moments (PWM：確率重み付き積率法) では，母数推計
値の偏り量と分散についての漸近理論値が利用でき，

Hosking et al. (1985) が紹介している．理論値はグラフの理解しにくい．波浪の極値統計では 4(1) で説明したよ

形で掲載されているので，図から読み取った概略値で検うに，代表波高の絶対値によって尺度母数や位置母数が討する．

規定されるので，データの座標軸の移動という操作はま形状母数の偏り量 (Bias) は，それに標本の大きさ n をず考えられない．

乗じた値が母集団の形状母数の値が大きくなるにつれて第 2 の理由として，カッパ分布から抽出したときの 3,減少する．[n×Bias] の値は，k = -0.3, -0.15, 0, および 0.15

4 の L-skewness, L-kurtosis の平均値の偏り（母集団の真(t3 = 0.38, 0.27, 0.17, および 0.08) でそれぞれ約 2.0, 0.9, 0.4,

値との差）について数値実験を行ったところ，n = 20 の0.2 である．これは今回の数値実験で得られた図-7 とほ

小標本では PP のほうが Unbiased よりもかなり大きな偏

とんど同じ値であり，数値シミュレーションの信頼度をり量を得たこと挙げている．この数値実験では，母集団保証するものである．

の諸元を 1 = 1, 2 = 0.2 に固定し，3 は 0～0.5 の間を 0.05

また，形状母数の偏り量の分散については，標本の大間隔，4 は 0 ～0.3 の間を 0.05 間隔で変化させている．

きさを乗じた値がk = -0.3, -0.15, 0, 0.15 でそれぞれ約 1.43,

カッパ分布は，尺度・位置母数の他に分布形状に係わる

0.75, 0.55, 0.48 である．これを標準偏差と標本の大きさ

二つの母数を持つ 4 母数分布であり，3, 4 の変化に応じの平方根の積に換算すると，それぞれ 1.2, 0.87, 0.74, 0.69て分布形状も大幅に変化していたと想定される．

となり，図-8 の結果とほぼ一致する．

今回の数値実験においても，L-モーメント比の推定に

なお，Hosking ほか (1985) の論文では一般化極値分布

関して = 0,  = 0 の Unbiased Position が極めて優れた不偏において尺度・位置母数をそれぞれ A = 1, B = 0 に設定し，性を示しており，Hosking・Wallis (1995) の結論を裏付け

形状母数を幾つかに選んだときの確率波高の推計値の偏ている．

り量や分散の漸近理論値が示されている．さらに，尺Hosking・Wallis (1995) は確率波高の偏り率についても

度・位置母数の推計値の偏り量や分散についての数値実幾つかの検討結果を示している．形状母数 k = -0.2 の

験結果も示している．

GEV 母集団 (A = 0.231, B = 0.81) から抽出された n = 50 の

ただし，尺度・位置母数は次元量であり，たとえば k

標本に関し，非超過確率 0.999 に対する確率波高= 0 のときには Hmean = 0.5772，平均発生率を  = 1 とする(quantile) の偏り率が PP では 3.5%，unbiased では 4.5 % のと H10 = 1.195, H50 = 1.592 であり，裾長度は 50 = 1.33 と高例を示している．（今回の数値実験は母数値が異なるの

24

波の極値としては例を見ない大きな値である．PWM 法の信頼度を高めるとともに，全体としてのデータ数を使や L-モーメント法を扱った研究では A = 1, B = 0 にに設定

うことによって推計値の信頼区間の幅を狭めることを目しても一般性を失わないとする論文 ( たとえば，的としている．Hosking・Wallis (1997) による地域頻度解

Landwehr ほか [1979], Hosking ほか [1985] ) をみかけるけ

析の手順は次の通りである．

れども，裾長度の観点からはそうした取り扱いには一般

1)性を認めることはできない．したがって，Hosking ほか

地理的な位置その他を勘案し，対象とする地域を選定してその地域内の極値データを収集する．

(1985) の論文に記載されている確率波高および分布関数の尺度・位置母数の推計値の偏り量について比較検討す
ることは困難である．

2)

各地点のデータをその平均値で除し，無次元化する．

3)

無次元化された標本について，単一地点と同じ手順
で L-モーメント法による解析を行う．

4)

１０．高波の極値統計解析に L-モーメント法を適用する際の手順について

各地点の L-CV, L-skewness, および L-kurtosis の値について各地点のデータ数を重みとする加重平均を行い，地域としての 3 種類の L-モーメント比を算出する．

(1) 単一地点の資料の解析

5)

各地点の 3 種類の L-モーメント比の地域全体の値（重心）からの距離を計算し，その距離が著しく大

L-モーメント法による極値統計解析の手順については

きければその地点は不整合であるとして当該地域か3.で述べたけれども，高波資料の標本に対する解析が行ら除外する．これは Discordance Test (不整合性テス

われ，所定の再現期間に対する確率波高が推計された後ト)と名付けられており，その作業指針は Hosking・の処理について補足する．

Wallis (1997) が示している．

これまでに述べてきたように，極値統計解析で得られ6)

た推計値は常に真値からのなにがしかの偏りを伴っておさらに，各地点の 3 種類の L-モーメント比を用いて地域内の標本の均質性をチェックし，地域の均質性

り，またある幅の分散を持つ．高波の母集団が既知でが確保できないときには幾つかの地点を除外あるいあったとしても真値そのものは把握不能である．ましては入れ替えて組み直す．これは Heterogeneity Test

母集団自体が不確定な現時点においては，本研究で取り(不均質性テスト) と名付けられている．

扱ったように，一般化極値分布，一般化パレート分布，

7)およびワイブル分布を母集団の候補として考え，そのう

均質性が確認された地域について，地域全体の 3 種
類の L-モーメント比を使って地域共通の分布関数をちで標本に最も適合した分布関数を母集団とみなすこと

選定する．

になる．

8)

この場合の確率波高の偏り率は図-12～14 の図表を参地点毎の極値分布関数は，形状母数は地域共通の値を用い，尺度母数と位置母数は地域共通の母数値

照し，分布関数毎に形状母数と標本の大きさに応じた適に各地点の平均値を乗じて換算する．

切な内挿によって推定する．これらの図表は裾長度が9)50 = 1.2 に対するものなので，裾長度がこれから外れる

地点毎に換算された分布関数の母数を使い，各種の再現期間に対する確率波高を算定する．

場合には一つの参考値として扱わざるを得ない．

10) 確率波高の信頼区間等については，地域共通の分布

確率波高の推計値の信頼区間については，まず裾長度関数を使って数値シミュレーションを行い，変動範を式 (38) によって算定する．そして，式 (48) と表-8 を

囲を調べて推定する．

用いて確率波高の変動係数を計算し，それから標準偏差を推計する．確率波高の推計値が正規分布に従うと仮定(3) 新たな不均質性指標の提案することにより，90% 信頼区間であれば推計値の上下に標準偏差の 1.64 倍を取った幅として定められる．標準偏Hosking・Wallis (1997) は不均質性指標として 3 種類の差の1 倍であれば，68% 信頼区間となる．当然ながら，指標を提示している．

この信頼区間は推計値の偏り量だけ補正したものである．
なお，信頼区間の試算例は 8.で紹介している．

H

(V  V )

V

(49)

上式中のV は次式で定義される量であり，V とV はカッ
パ分布を当てはめて行う数値シミュレーションによって(2) 複数地点の資料に対する地域頻度解析

求められるV の平均値と標準偏差である．

地域頻度解析は，均質とみなされる地域内の複数地点の資料を共通に取り扱い，その地域の極値の母集団推定

25





 N n t (i )  t R 2  i 1i i

V 

N
n



i 1i i

不均質とみなされるのかについては疑念が抱かれる．

1/ 2

Hosking・Wallis (1997) の不均質性指標は 4 母数型の極(49)値分布であるカッパ分布の数値シミュレーションに基づいているけれども，カッパ分布を導入する意義は明瞭でここに，N は地点数，ni はi 番目の地点のデータ数， t

(i)

ない．

R

はその地点の L-モーメント 比に係わる統計量，t はそのむしろ，今回の数値実験によって L-モーメント比の変統計量に関してカッパ分布当てはめのシミュレーション

動係数の推定式が表-6 のようにとりまとめられたので，で得られた地域平均の値である．

これを用いて各地点の L-モーメント比の標準偏差を算定

不均質性指標は，3種類の統計量について検討する．第する．そして，各地点の L-モーメント比の地域平均値か

1は，2 次の L-モーメントと1次の L-モーメントの比であらの差をこの標準偏差で除したものを偏心率と名付け，るL-CV [t] である．すなわち，ある地点の t の値が他より地点毎の偏心率の分散の度合いが正規分布に近ければそも著しくかけ離れていると，その地域は不均質と見なさの地域を均質とみなし，偏心率の分散が正規分布から外れる．第 2 の統計量としては，L-CV と L-skewness の相関れていれば不均質とみなすことを提案する．
偏心率は，L-CV, L-skewness, および L-kurtosis について

プロット図において，ある地点が全体の重心からどれだ
け離れているかを検討する．第 3 の統計量としては，L-は次式で計算する．
 k ,i  [ k ,i   ( k )] /  ( k )kurtosis と L-skewness の相関プロット図における重心から: k  1, 2, 3; i  1, 2,  , Nの距離を用いる．なお，第2 と第3 の不均質指標について

ここに，1, 2, 3 はそれぞれ L-CV[t], L-skewness[t3], およ

は，V として2 次元座標上の距離を用いて算定する．

びL-kurtosis[t4] を代表する．さらに，2 次元座標上の距離地域が均質であるか否かを判断する不均質性指標の閾値としてHosking・Wallis (1997) はH が1以下であれば均質に関する偏心率として次の二つを定義する．

性に問題はなく，H が 2 を超えると不均質である公算が大きいとしている．
表-9 解析対象地点のデータ個数とL-モーメント比
地点名
留 萌
瀬 棚
深 浦
秋 田
酒 田
新潟沖
輪 島
金 沢
福 井
鳥 取
浜 田
全地域

データ
個数, n
554
330
408
314
518
252
298
387
248
320
400
4029

L-CV

L-skew

L-kurt

0.09640.09270.10070.10670.10960.10610.0826
0.09610.11180.08980.0917
0.09856

0.24770.27490.27800.26630.25890.29730.25140.22410.25740.29520.27720.264280.11830.13930.14340.13230.12990.13330.13060.10020.12370.15340.13260.12978(51)

[   ( 1 )] 2  [ 2,i   ( 2 )] 2  4,i   1,i 2 ( 1 )   2 ( 2 )



1/ 2

[   ( 2 )] 2  [ 3,i   ( 3 )] 2  5,i   2 ,i
 2 ( 2 )   2 ( 3 )



1/ 2

(52)

(53)

式 (52) の4 はL-CV とL-skewness の相関プロット図上の偏心率，5 はL-kurtosis とL-skewness の相関プロット図上の偏心率を表す．
こうした各地点の偏心率の分散の度合いを正規分布と比べるには，偏心率の最大値を用いるのが簡便である．
すなわち，k  max{| k,i | : i  1, 2,, n} : k  1, 2, , 5(54)

正規分布であれば，偏心率の最大値が1.6 を超える確率は11%，2.0 を超える確率は4.6% であるので，5 種類の偏心

著者らが先に行った日本海沿岸の高波の地域頻度解析

率について正規分布と比較しながら均質性を調べること(合田ら, 2010) では，3 通りの不均質性指標のうちで L-ができる．

CV に係わる H (1) が最もシビアであり，表-9 に示す留萌

表-9 の日本海沿岸 11 地点の高波資料について，母集団

から浜田までの 11 地点全体ではH (1) = 5.27 となって，こがワイブル分布であるとして L-モーメント比に関する偏れらを一つの地域に含めることは不可との判断が下され心率を計算すると，表-10 のような結果が得られる．日本

た．また，輪島から浜田までの 5 地点をまとめて取り

海全域を一つの地域とみなした場合には，表-10 に示すよ

扱っても不均質性指標が H (1) = 4.49 と高い値であり，一

うに L-CV に関する偏心率 1 が輪島で最大値 3.52 となり，つの地域に含めることができないと判断された．しかし

やはり全域を一つの地域とみなすことが無理であること

ながら，11 地点の L-CV には輪島の 0.0826 から福井のが示される．

0.1118 までの狭い範囲に入っており，これでこの地域が26

3 種類の L-モーメント比の中では，L-skewness が母集

表-10 日本海沿岸の高波の均質性指標（全域共通）

団の形状母数を規定するだけでなく，表-6～8 にまとめ

不均質性指標としての偏心率
地点名

1

2

3

4

5

た各種の変動係数・標準偏差の独立変数として機能して(L-CV)(L-skew)(L-krt)(CV-skw)(krt-skw)

いる．その意味で，これに関する最大偏心率 2 を最も重留 萌

-0.677-0.834-0.801

0.823

0.823

瀬 棚

視すべきではないかと思われる．今後はさらに多くの極-1.399

0.420

0.438

0.428

0.428

深 浦

0.522

0.604

0.675

0.636

0.636

値データに適用し，最大偏心率の許容限界について検討

秋 田

1.682

0.079

0.117

0.096

0.096

したいと考えている．

酒 田

2.927

-0.266

0.008

0.208

0.208

新潟沖

1.307

1.159

0.150

0.897

0.897

輪 島

-4.314-0.475

0.038

0.372

0.372

金 沢

-0.696-1.666-2.014

1.767

1.767

福 井

2.387

-0.234-0.273

0.249

0.249

鳥 取

-2.018

1.224

0.980

1.114

1.114

浜 田

-1.824
4.314

0.561
1.666

0.147
2.014

0.444
1.642

0.444
1.767

最大値

１０．まとめと今後の課題
今回のモンテカルロ法による一連の数値シミュレー
ションによって明らかになった諸点をまとめると次のよ
うになる．

1) L-モーメント法においてもプロッティングポジショ
ンの選定は重要な課題であり，不適切なものを用い

表-11 日本海沿岸の高波の均質性指標（南北2 分割）

ると分布関数の母数や確率波高に大きな偏り (Bias)

不均質性指標としての偏心率をもたらす．

1

2

3

4

5

(L-CV)(L-skew)(L-krt)(CV-skw)(krt-skw)

2) Hosking・Wallis (1977) が “unbiased position” と称して

留 萌

-1.745-0.961-0.919

0.988

0.947

瀬 棚

推奨するプロッティングポジション式は，L-モーメ-2.203

0.319

0.360

0.477

0.337

深 浦

-0.294

0.492

0.591

0.487

0.538

ントや分布関数の母数の推計値の偏りがきわめて小

秋 田

0.986

-0.019

0.037

0.180

0.027

酒 田

2.042

-0.390-0.096

0.534

0.311

3) ガンベル分布と正規分布に対して最小2 乗法を用い

新潟沖
北最大値

0.729
2.203

1.070
1.070

0.078
0.919

1.058
1.058

0.825
0.825

るときに有効であった Gringorten 式と Blom 式は L-

輪 島-3.010-0.340

0.151

0.528

0.281

金 沢

0.657

-1.514-1.848

1.501

1.611

したがって，この 2 式は L-モーメント法には推奨で

福 井

3.253

-0.110-0.164

0.611

0.133

鳥 取

きない．

-0.905

1.368

1.081

1.357

1.240

浜 田

-0.546
3.253

0.721
1.514

0.277
1.848

0.717
1.501

0.585
1.611

地点名

南最大値

さいことが確認された．

モーメント法に用いると不偏性を満足しなかった．

4) 今回の数値実験で導入した新しいプロッティングポ
ジション式 ( = 0.45,  = 0) は，偏りが僅少であり，“unbiased position” と同様に不偏性をほぼ満足する．
このため，日本海沿岸の11 地点を留萌から新潟東まで

また，確率波高推計値の偏りが “unbiased position” よの北地域6 地点，輪島から浜田までの南地域5 地点に分割

りもやや小さいので，L-モーメント法のプロッティ

して偏心率を計算すると，表-11 のような結果が得られる．

ングポジション式としてこれを推奨する．

これによれば，北地域に関しては1～5 のどの最大偏心

5) L-モーメント法であっても標本が小さいと母集団以率も1.8 以下となり，均質性に問題はないと判断される．

外の分布関数に当てはめられる割合が大きくなり，

ただし，南地域に関しては1 が金沢で2.65 となり，均質

また母数値や確率波高の信頼区間の幅が非常に大き性についての疑問が残る．ただし，2～5 に関してはいくなる．標本のデータ個数 20 では不十分であり，ずれもも1.8 以下であるので，南地域を一体として解析す

できれば 100 以上が望ましい．

ることも許容されると思われる．

6) 極値分布関数が一般化極値分布 (GEV)，一般化パ式 (51) ～(54) による不均質性指標の一つの利点は，L-

レート分布 (GPA)，およびワイブル分布である母集モーメント比の値が地域内の他地点からかけ離れた地点

団から抽出されたそれぞれの標本に対し，L-モーメ

を直ちに特定できることである．Hosking・Wallis (1997)

ント比と分布関数の形状母数の推計値に対する偏りの不均質性指標の場合には，全体として均質性に問題が率の算定図表が作成され，さらにそれらの標準偏差あっても，その是正は地点の入れ替えを試行錯誤によっを推定するための経験式が求められた．

て行わなければならず，問題地点の抽出に手間取ること

7) 3 種類の極値分布関数に対し，それぞれを母集団とが多い．その点，新しく提案する不均質指標としての偏して設定したときの確率波高の推計値の偏り率の図

心率であれば，地域内の値手の組み替えが容易に行える．

表が作成された．

27

布関数に関する実証的研究，海岸工学論文集，第45巻， pp.
211-215.

8) 母集団として 3 種類の極値分布関数の中から最適等分布を選定する場合の当てはめ的中率の図表が作成

合田良実・竹下直樹・永井紀彦 (1999)：太平洋南岸の極値波高の母分布関数について，海洋開発論文集，Vol. 15，土木学会
海洋開発委員会，pp. 327-331.

され，標本の大きさが 100 程度未満では当てはめ的中率が高くないことが確認された

9) 母集団が未知の場合の確率波高の変動係数に関する

合田良実・久高将信 (2009a)：高波の極値統計解析に対する Lmoments 法の適用と分布関数の選択について，ECOH/YG 技術論文 No. 10, 2009年5月7日，12p.

経験式が導かれた．この経験式は，裾長度，標本の大きさ，基準化変量の関数であり，数値実験によるデータの変動係数をかなりの確度で推定可能である．

10) 地域頻度解析を行うための地域分割に際して検討すべき不均質性指標に関して，L-モーメント比の変動

合田良実・久高将信・河合弘泰 (2009b)：L-moments 法を用いた波浪の極値統計解析について，海岸工学論文集，第56巻，pp.161-165.
合田良実・久高将信 (2010)：L-moment法による日本海沿岸の高波の地域頻度解析について，ECOH/YG 技術論文 No. 12, 2009
年4月12日，13p.

係数の経験式を用いて推定される標準偏差に基づく
新しい不均質性指標を提案し，日本海沿岸の高波資

外山奈央子・水野 量 (2002)：L-moments を用いた地域頻度解析
による全国アメダス地点における確率降水量の推定，気象庁
研究時報，54巻 5－6号合併号，pp. 55-100.

料についてその適用事例を説明した．
今回の数値実験では，3 種類の分布関数の中から最適

Blom, G. (1965): Statistical Estimates and Transformed Beta–Variables,
John Wiley & Sons, New York, Chapter 12.

合のものを選ぶ指標として極値データとその推計値との間の相関係数を用い，相関係数が最大のものを最適合とBoes, D. C., J.-H. Heo, and J. D. Salas (1989): Regional flood
quantile estimation for a Weibull model, Water Resources Res.,
Vol. 25, No. 5, pp. 979-990.

する方式を用いた．母集団から抽出された標本に最も適合する分布が母集団の分布関数に的中する割合は，標本

Coles, S. (2001): An Introduction to Statistical Modeling of Extreme
Values, Springer, 208p.
Chowdhury, J. U., J. R. Stedinger, and L.-H. Lu (1991): Goodnessof-fit tests for regional generalized extreme value flood
distributions, Water Resources Res., Vol. 27, No. 7, pp. 17651776.
Greenwood, J A. , J. M. Landwehr, N. C. Matalas, and J. R. Wallis(1978): Probability weighted moments: Definition and relation
to parameters of several distributions expressable in inverse
form, Water Resources Res., Vol. 15, No. 5, pp. 1049-1054.

が大きくなるにつれて向上する．しかし，一般化パレー
ト分布で形状母数の値が負の値の場合に，標本の大きさ
が n = 100 とかなり大きいときでも的中率が 30% 前後にとどまる場合があり，最適合の判断基準をさらに工夫する必要があろう．
また，今回新しく提案した地域内の複数地点の不均質性指標については今後さらに多くの事例について検討を進め，その適用性を向上させたい．

自然現象の極値については未解明の問題が数多くあり， Gringorten, I. I. (1963): A plotting rule for extreme probability paper,
J. Geophys. Res., Vol. 68, No. 3, pp. 813-814.
今後も多様な研究を進めていくことが肝要である．
Hosking, J. R. M. (1990): L-moments: Analysis and estimation of
なお，この数値実験を実行するに当たっては，プログ
distributions using linear combinations of order statistics, J.Roy.
ラムのバグ解消のために（株）エコー参事の久高将信氏
Statistical Soc., Series B, 52, pp. 105-24.
Hosking, J. R. M. and J. R. Wallis (1987): Parameter and quantileに幾度となくご助力頂いた．ここに記して謝意を表する
estimation for the generalized pareto distribution, Technometrics,次第である．
Vol. 29, No. 3, pp. 339-349.
Hosking, J. R. M. and J. R. Wallis (1995): A comparison of unbiased
参考文献
and plotting-position estimators of L moments, Water
Resources Res., Vol. 31, No. 8, pp. 1049-1064.
北野利一・間瀬 肇・喜岡 渉・矢野陽一郎 (2002)：一般化パ
Hosking, J. R. M. and J. R. Wallis (1997): “Regional Frequency
レート分布による極値波浪解析－拡張形状母数の推定－，海Analysis,” Cambridge Univ. Press, 224p.
岸工学論文集，第49巻， pp. 161-165.
Hosking,
J. R. M., J. R. Wallis, and E. F. Wood (1985): Estimation of
合田良実 (1988)：極値統計におけるプロッティング公式ならび
thegeneralized
extreme-value distribution by the method ofに推計値の信頼区間に関する数値的検討，港湾技術研究所報probability-weighted
moments, Technometrics, Vol. 27, No. 3,
告，第27巻 第1号，pp. 31-92.
pp. 251-261.
合田良実 (2002)：設計波高に係わる極値統計分布の裾長度パラ
Landwehr, J. M., N. C. Matalas, and J. R. Wallis (1979a): Probability
メータとその意義，海岸工学論文集，第49巻， pp. 171-175.
weighted moments compared with some traditional techniques
合田良実 (2008)：「耐波工学 – 港湾・海岸構造物の耐波設計–
in estimating Gumbel parameters and quantiles, Water」，鹿島出版会，第13章，pp. 347-354.
Resources Res., Vol. 15, No. 5, pp. 1055-1064.
合田良実・小野澤昌己 (1990)：極値 II 型分布の特性と推計値のLandwehr, J. M., N. C. Matalas, and J. R. Wallis (1979b): Estimation
信頼区間，土木学会論文集，No.417/II-13(ノート), pp. 289-292.
of parameters and quantiles of Wakeby distributions, Water
合田良実・小長谷 修・永井紀彦 (1998)：極値波浪統計の母分
Resources Res., Vol. 15, No. 5, pp. 1361-1379. (Correction,

28

distribution, J. Hydrology, Vol. 99, pp. 271-280.
van Gelder, P. H. A. J. M., J. De Ronde, N. W. Neykov, and P.
Neytchev (2000): Regional frequency analysis of extreme wave
heights: trading space for time, Coastal Engineering 2000 (Proc.
26th ICCE, Sydney), ASCE, pp. 1099-1112.

Water Resources Res. Vol. 15, p. 1672.)
Ma, Q.-S., Y.-B. Li, and J. Li (2006): Regional frequency analysis ofsignificant wave heights based on L-moments, China Ocean
Engineering, Vol. 20, No. 1, pp. 85-98.
Sinclair, C. D. and M. I. Ahmad (1988): Location-invariant plotting
positions for PWM estimation of the parameters of the GEV

29