ECOH/YG 技術論文 No. 14

2010年12月12日版
2011年3月29日追記

年最大値の極値統計量に及ぼす標本サイズの影響に関する数値実験
Numerical Experiments for Examination of Sample Size Effect on Return Wave Heights
Estimated from Samples of Annual Maximum Series Data

合田良実
Yoshimi GODA
土木学会名誉会員 工博 横浜国立大学名誉教授 (株) エコー（〒110-0014 東京都台東区北上野 2-6-4）

Distributions of annual maximum wave heights theoretically approach to the generalized extreme-value (GEV)
distribution when the sample size is sufficiently large. A series of Monte Carlo simulations have confirmed it for samples of
10,000 years long extracted from the peaks-over-threshold data from the parents of the GEV, generalized Pareto (GPA),
and Weibull distributions. However, samples of annual maximum data with duration shorter than several hundred years
deviate from the GEV considerably and estimated return wave heigthts are smaller than the parent values for the return
period longer than 50 years. Empirical prediction formulas for the coefficient of variation of return wave heights arederived for the three parent distribution functions.

Key Words: Extreme wave statistics, annual maximum series data, GEV distribution, GPA distribution, Weibull
distribution, quantile estimation

１．まえがき著者は極大値の事象がポアソン分布することを前提として年最大値の分布関数を極大値の分布から導き，極水文統計をはじめとする自然現象の極値事象を扱う

大値の分布から推定される確率波高の再現期間が年最分野では，ある閾値を超える極値のすべてを対象とす

大値分布のものよりも 1/2 年短いだけであることを示

るのか，あるいは各年のうちの最大値に限定して解析している（合田 [2008], p.335 参照）．

すべきであるかについて，早くから議論されてきた．

波浪の分野では計器観測の歴史が短いこともあって，前者のデータは極値時系列データあるいは極大値資料

極値解析はもっぱら極大値資料に基づいてきた．しかと呼ばれ，後者は年最大値資料と呼ばれる．

し，わが国の沿岸波浪観測のデータが蓄積されるにつ統計解析の作業量から言えば毎年最大値資料のほうれて 40 年に近い期間にわたる年最大波高のデータが解が楽であるためか，日雨量や高水流量などでは年最大

析可能となってきた．こうしたデータについて試算し値資料を扱うことが多い．しかし，ある年の年最大値てみると，極大値資料と年最大値資料とでは分布関数として取り上げられなかった極大値であっても，他のに有意な違いが見られ，100 年確率波高も後者のほう

年の年最大値よりも大きいことは珍しくない．そのたが低いことが少なくない．このため，一連の数値実験

め，そうした極大値を無視するのは不適切であり，極を行って年最大値資料に基づく極値分布ならびに確率

大値資料に基づいて100年確率統計量などを推定すべ

波高の特性を解明することにした．

きあるとの意見がある．これについては，Langbein

この研究は，先に取りまとめた数値実験の成果（合(1949) が簡単な確率計算に基づいて極大値と年最大値

田 [2010]，以下では前報と称する）を引き継ぐものでに対する T 年確率統計量の再現期間の差が小さいことあり，データ解析は Hosking (1990) と Hosking・Wallis

を例示し，実用的には無視できると説明した．また，-1-

(1997) による L-モーメント法を用いている．したがって，解析法の詳細については前報を参照されたい．

1.0

本報の全体構成は次の通りである．まず，2. におい
Non-exceedance probability, Pて極大値資料から年最大値資料を導くときの前提となる高波の発生回数のポアソン分布について，酒田港およびアラス沿岸 Kodiak 島沖の波浪データを用いて吟味する．次に，極大値分布と年最大値分布の理論的関係について 3. で説明する．今回の数値実験の条件や計算方法については 4. で述べる．波浪条件としては 50 年
確率波高を H50 = 10.0 m に設定し，裾長度パラメータは

γ50 = 1.2 に固定する．極大値分布の母集団としては，一

0.8

0.6

0.4

Obs. (Hc = 3.0 m)
Theory ( λ = 28.0)
Obs. (Hc = 4.0 m)
Theory ( λ = 15.1)
Obs. (Hc = 5.0 m)
Theory ( λ = 7.50)
Obs. (Hc = 6.0 m)
Theory ( λ = 3.09)

0.2

般化極値分布 (GEV)，一般化パレート分布 (GPA) およ0びワイブル分布の 3 種類を用い，高波の発生数が年平

0

10
20
30
Occurrence number, n

均 3 および 10 回を想定して毎年の最大波高を抽出する．
数値実験では期間が十分に長いデータとして期間長40

図-1 酒田港の高波の発生回数のポアソン分布との比較

が K = 10,000 年の年最大値を扱った場合を 5. で記述する．現実ではあり得ないけれども，このような超長期

1.0

あっても一般化極値分布で近似できることを例示する．
しかし，期間長が数十年程度の実際の資料であれば，標本サイズの影響によって一般化極値分布への収束が保証できないことを 6. で説明し，こうした小標本の年最大値資料の統計的性質の幾つかについて紹介する．
確率波高の偏りについては 7. で記述し，年最大値資料から推定される確率波高の変動係数について導いた経験的推定式を8. で提示する．さらに，9. で幾つかの現

Non-exceedance probability, P

間にわたる年最大値資料であれば，母集団がいずれで
0.8

0.6

0.4
Obs. (Hc = 6.0 m)
Theory (λ = 3.9)
Obs. (Hc = 7.0 m)
Theory (λ = 2.25)
Obs. (Hc = 8.0 m)
Theory (λ = 1.3)

0.2

地観測波浪データを用いて，本研究の適用妥当性を検0討する．最後に 10. において全体のまとめを述べる．

0

2

4

6

8

10

Occurrence number, n

２．高波の発生回数のポアソン分布

図-2 Kodiak 島沖の高波の発生回数のポアソン分布との比較

統計の分野では，ポアソン分布という確率分布がし

トワーク (NOWPHAS, Nationwide Ocean Wave informa-

ばしば用いられる．これは，ある事象が一定期間内にtion network for Port and HArbourS) のうちの酒田港の高発生する平均回数が分かっている場合，その期間内に

波資料ならびにアラスカ沿岸 Kodiak 島沖の高波の追算

その事象が発生する回数 m がどのような確率で分布す資料を用いて検討した．酒田港の資料は 1970 〜2008 年

るかを表すものである．高波の場合であれば，年間のの 39 年間の観測データであり，有効統計年数 34.4 年の平均発生回数 (λ で表示) が対象である．

間に有義波高 3.0 m 以上の波高極値を 955 回記録してい一定期間内に事象が m 回発生する確率 P(m) は次のる．一方，Kodiak 島沖の波浪追算は 1956 〜1975 年のように表される．

1
P ( m) = e − λ λ m
m!

20 年間を対象として米国陸軍工兵隊海岸工学研究セン: m = 0, 1, 2, L

(1)

ターが実施したものであり，各年の有義波高極大値のある期間内に 1 回も発生しない確率は P(0) = e-λ であり，

λ = 1, 2, 3 に応じてそれぞれ P(0) ≈ 0.368, 0.135, 0.050

時系列は合田 (2008, p. 329) に採録されている．
酒田港の資料について高波極大値選定の閾値を Hc =3.0, 4.0, 5.0, および 6.0 m と高めていくと，極大値の個数となる．

が 955, 518, 258, および 106 と減少する．それぞれの波

高波の発生回数がポアソン分布で近似できるか否か

高閾値を超える高波の個数を年ごとに数え，発生回数について，わが国の港湾関連部局の全国沿岸波浪ネッの累積度数分布を作成して非超過確率を計算した．た2

だし，年によって欠測期間の長さが異なるため，年間

式 (1) で表されるので，年最大値の非超過確率は次のの発生回数は欠測期間を考慮して補正している．このように求められる．
∞

補正された高波極値の平均発生率は，4 通りの閾値にF ( x) = Pr[ H max < x] = ∑ P(m) G m ( x)対して λ = 28.03, 15.06, 7.50 および 3.09 回 / 年と算定され

m =0

(3)

∞

1
= ∑ e −λ [λG (m)]m = exp{−λ[1 − G ( x)]}m = 0 m!

る．そこで，式 (1) のポアソン分布の確率を求めて非超過確率を算定し，実測値と比較したのが図-1 である．
図中の凡例で Theory と記したのはポアソン分布の計

ここで Pr[X] は事象 X の起こる確率を表す．なお，式

算値を表している．ポアソン分布は離散関数であると(3) の第 2 行の演算には指数関数のべき級数展開を利用

ころからこのような非超過確率の形で比較したものでしている．

ある．実測値と理論値の両者とも累積回数は整数であ

著者が既に紹介しているように（合田 [2008], p.

り，整数値の ±0.5 に広がる水平線としてデータを表示335），R 年確率統計量は非超過確率 1 – 1/R に対応するしている．たとえば，閾値が Hc = 4.0 m のときには年間値として定義される．すなわち，発生数が 8 個の年が 1回，9個の年が 2 回あり，年間発1(4)
R
生数が 10 個の年はなかった．したがって，発生回数 9
この非超過確率を式 (3) と等値すると，極大値分布のと10 の累積度数が 3，非超過確率が 0.094 であり，横軸の 9 〜10 の個所にこの確率の相当する水平線が描かれ非超過確率が次のように表される．
ている．実測値と理論値の非超過確率に若干のずれは1 ⎛ 1⎞(5)
G ( x R* ) = 1 + ln⎜1 − ⎟λ ⎝ R⎠あるものの，全体として両者は合致しているといえる．
図-2 は Kodiak 島沖の高波のデータについて発生回数
ここに， x R* は年最大値としての R 年確率統計量 xR にの非超過確率を比較したものである．波浪追算資料が
波高閾値を Hc = 6.0 m として提示されたものであり，20対応する，極大値としての値である．式 (5) はさらに R
年間で 78 個の高波が報告されていて，そのうち 1 年はが 1 よりも十分に大きい場合について，自然対数のべ
無発生である．平均発生率は λ = 3.90 回 / 年である．閾
き級数展開を使って近似計算することによって，次のように書き直される．1
値を 7.0 および 8.0 m に高めると，平均発生率が 2.25 お1
よび 1.3 回 / 年に低下する．酒田港に比べて高波の総数(6)
G ( x R* ) ≅ 1 −λ ( R − 1 / 2)が少ないために実測回数と理論予測回数との間の差異
F ( xR ) = 1 −がやや大きく見えるけれども，全体としてはほぼ一致

すなわち，極大値としてのR 年確率統計量 x R* は年最しているといえる．
この 2 地点の事例から見ても，高波の発生回数をポ大値よりも 1/2 年だけ短い再現期間に対応する値であアソン分布で近似することには妥当と考えられる．

る．

３．極大値と年最大値の分布関数の相互関係(2) 極大値がワイブル分布に従う場合の年最大値分布(1) 極大値と年最大値の一般的関係

ワイブル分布の分布関数は次のように表される．

⎡ ⎛ x − B ⎞k ⎤G ( x) = 1 − exp⎢− ⎜⎟ ⎥ : B<x<∞⎣⎢ ⎝ A ⎠ ⎦⎥いま1 年間に発生する高波の極大値を発生順にH1, H2,
L, Hm と表示すると，その年の最大波高 Hmax は次のように定義される．

H max = max{H 1 , H 2 , L, H m }(7)

ここに，k は形状母数，A は尺度母数，B は位置母数で(2)ある．式 (7) を式 (3) に代入すると，極大値がワイブル
分布に従う母集団から抽出された年最大値の分布関数

ここで高波の極大値の分布関数すなわち非超過確率
を G(x) で表す．1 年間に m 個発生した高波の極大値の1

うちの最大波高 Hmax がある値 x を超えない事象は，H1,
H2, L, Hm のいずれもが x を超えない事象であり，そ

1
1
1⎛
1 ⎞
−
= − ⎜1 +⎟R 2R 2R ⎝ 2R ⎠
1
1
1
1
1
=− ⋅
≅− ⋅
=−
R 1 /(1 + 1 / 2 R)
R 1 − 1 / 2R
R − 1/ 2ln(1 − 1 / R) ≅ −の確率は [G(x)]m で表される．極大値の発生個数 m は 0〜∞ のうちの任意の値を取り得る．高波の発生個数がポアソン分布に従うとすると，それぞれの生起確率は3

上式中の尺度母数 A* と位置母数 B* は一般化パレート分が次のように導かれる．

⎧⎪⎡ ⎛ x − B ⎞ k ⎤⎫⎪F ( x) = exp⎨− λ exp⎢− ⎜⎟ ⎥⎬⎪⎩⎣⎢ ⎝ A ⎠ ⎦⎥⎪⎭

布の母数から次のように算定される．

⎧⎪ A * = Aλ − k⎨ *⎪⎩ B = B + A(1 − λ − k ) / A

(8)

(14)

⎧⎪
⎡ ⎛ x − B ⎞k⎤⎫⎪
= exp⎨− exp⎢− ⎜⎟ + ln λ ⎥⎬⎪⎩⎢⎣ ⎝ A ⎠⎥⎦⎪⎭Madsen et al. は数値シミュレーション手法を用い，式

ワイブル分布の中で形状母数が k = 1 のときは指数分(11) の分布に従う極大値資料と式 (13) による年最大値

布を表す．極大値が指数分布の場合には，年最大値は資料のそれぞれから推定される確率統計量の分散の大次のガンベル分布に従うことになる．

きさを比較検討している．ただし，あくまでも式 (11),⎧⎡ x − B* ⎤⎫F ( x) = exp⎨− exp⎢−
⎥⎬
A ⎦⎭
⎣
⎩

(13) の分布関数に基づく検討であって，一般化パレー
ト分布から年最大値の標本を作成してそれを統計解析(9)する作業は行っていない．

ここに，
B* = B + A lnλ４．数値計算の条件と計算方法(10)(1) 発生させる高波の特性

すなわち，尺度母数は極大値分布と同じ値を保ち，位置母数が母集団の値よりも A lnλ だけ大きな値を取るこ

数値シミュレーションによって高波の標本を多数作とになる．

成する際には，高波として現地の波浪特性を反映させ

式 (9), (10) は先に Cunnane (1973) が提示したものであ

たものであることが必要である．高波の極値分布の特り，Cunnane は指数分布と式 (9) のガンベル分布の確率性として重要なのは，合田 (2002) が提唱した裾長度パ

波高推定値の分散の大きさを比較している．ただし，ラメータである．これは 10 年確率波高 H10 に対する 50指数関数で表される母分布関数から抽出される個別の年確率波高 H50 の比として次のように定義される．

標本について吟味しているわけではない．

γ 50 = H 50 /H 10(3) 極大値が一般化パレート分布に従う場合の年最大

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

値分布パラメータ γ50 が 1.13〜1.27 の値をとっている（合田2008, 375頁）．前報では裾長度パラメータとして γ50 =一般化パレート分布 (GPA) の関数形は次のようなも

1.10, 1.20, および 1.30 の 3 通りの値を用いたけれども，

のである．
1/ k
⎧ ⎡
x − B⎤⎪1 − ⎢1 − k
A ⎥⎦⎪G ( x) = ⎨ ⎣⎪1 − exp ⎡− x − B ⎤
⎢⎣
⎪⎩
A ⎥⎦

(15)

本研究では γ50 = 1.20 に固定した．その代わりに，後述: k ≠0のように平均発生率として λ = 3 および 10 回 / 年の 2 通(11)りの値を用いている．
また，50 年確率波高としては H50 = 10.0 m に設定する．

: k =0

これは数値の換算が容易なように区切りの良い数値としたものである．

このうち k = 0 は指数分布を表し，このときの年最大値

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

分布は式 (9) となる．

の尺度母数と位置母数は任意の値を与えることができk ≠ 0 の場合の年最大値分布は次のように導かれる．
1/ k⎧⎪ ⎡ k ⎛x − B ⎞ ⎤ ⎫⎪F ( x ) = exp ⎨− ⎢λ ⎜1 − k⎟⎬A ⎠ ⎥⎦ ⎪⎭⎪⎩ ⎣ ⎝

ず，次式によって一義的に決定される値を使わなければならない．

(12)

Madsen et al. (1997) はこの式をさらに書き換え，年最大値分布が次のような一般化極値分布となることを示し
ている．
1/ k
⎧⎪ ⎡
x − B * ⎤ ⎫⎪F ( x ) = exp ⎨− ⎢1 − k⎥ ⎬ :k ≠ 0
A* ⎦ ⎪⎪⎩ ⎣⎭1 − 1 /γ 50
y50 − y10

⎫
⎪
⎪
⎬
⎡
1 − 1 /γ 50 ⎤ ⎪B = H 50 ⎢1 −⎥
⎣ 1 − y10 / y50 ⎦ ⎪⎭A = H 50(16)

ここに，y10 と y50 は再現期間 R = 10年と50年に対する基(13)

準化変量である．

4

本研究では，極大値分布として一般化極値分布

なお，GEV 分布の k = 0.001 はガンベル分布を近似的(generalized extreme-value: GEV) ， 一 般 化 パ レ ー トに表現したものである．また，ワイブル分布の k = 1.0(generalized Pareto: GPA) 分布，およびワイブル分布のは指数分布であり，GPA 分布の k = 0 と共通である．た

3 種類の極値分布関数を用いる．これらに対する基準

だし，確認のために GPA 分布の k = 0.001 のケースも計

化変量は次式で計算される．

算に加えた．

1) 一般化極値分布 (GEV)：

表-1 高波の発生に使用する分布関数⎧{1 − [− ln(1 − 1 / λR )]k } / k , k ≠ 0
yR = ⎨
k =0⎩− ln(− ln(1 − 1 / λR )),

分布関数

(17)

形状母数 k

一般化極値 (GEV)

-0.45

-0.3

-0.15

0.001

0.15

一般化ﾊﾟﾚｰﾄ (GPA)

-0.2

-0.1

0.001

0.10

0.15

ワイブル

0.75

1.0

1.4

2.0

–

2) 一般化パレート分布 (GPA)：⎧[1 − 1 /(λR) k ] / k , k ≠ 0
yR = ⎨k =0⎩ln(λR),(3) 母分布関数の尺度・位置母数および統計量(18)

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

3) ワイブル分布：

数 A と位置母数 B が式 (16) によって一義的に決定され

y R = {ln(λR)}

1/ k

(19)

る．今回の数値計算では H50 = 10.0 m とし，平均発生率
を 2 通りに変えているので，この条件の下で表-1 の母(2) 年最大値の母集団である極大値の分布関数

分布関数の尺度・位置母数の値を計算したところ，表2〜4 の母数値が得られた．また，分布関数の平均値

各種の極値分布は， L-モーメント比 τ3 (L-skewness)µ[x]，標準偏差 σ [x]，歪み度 β1 ならびに尖鋭度 β2を求とτ4 (L-kurtoness) の相関関係によってその特性を表すことができる．図-3 は，GEV，GPA，およびワイブル分

めた結果も併せて示している．理論値が計算できるも布の特性曲線を示したもので，図中の ◆，■，● の記のはそれを記入したが，歪み度と尖鋭度については極号はそれぞれの分布関数において形状母数 k が特定の大値の全データの集計結果から計算した値である．
なお，これらの統計量の定義は次の通りである．ま

値をとる場合を示している．これらの形状母数は前報

た，尖鋭度を除く統計量の理論値は表-5 に提示した．
0.4
GumbelGEV (k < 0)
GEV (k > 0)
Weibull
Exponential
GPA (k > 0)
GPA (k < 0)

L-Kurtosis, τ4

0.3

L-Skewness vs.
L-kurtosis

GEV:k = -0.45

-0.30

GPA:
k = -0.2

-0.1

-0.15

0.2
0.15
2.8

Weibull:k = 0.75

0.2
0.4

0
0.1

1/ 2

(21)

1 ⎡1 N3⎤∑ (ξi − µ[ξ ]) ⎥⎦σ 3 ⎢⎣ N i

(22)

β2 =

1 ⎡1 N(ξ i − µ[ξ ])4 ⎤⎥∑4 ⎢σ ⎣N i⎦(23)

表-2 は極大値の母集団がGPA分布の場合である．この

0.6
0

(20)

β1 =

1.4
2.0

1 N
∑ξi
N i =1⎡1 N⎤
σ [ξ ] = ⎢ ∑ (ξ i − µ[ξ ])2 ⎥
⎣N i
⎦

1.0

0.001

0.1

µ[ξ ] =

0.2

0.3

0.4

うち，形状母数が -0.45 のケースは分布の裾が非常に長0.5

L-Skewness, τ

く広がるために，データ個数が 5 千万個と多くなって

3

図-3 高波の母分布関数の L-kurtosis と L-skewness の関係も尖鋭度を安定して計算できない．表中で *** の記号で用いたものである．

を付したのは，信頼できる数値が得られなかったケー

本研究では，年最大波高の標本を抽出する母分布関スである．形状母数が -0.30 のケースでも尖鋭度につい数として GEV分布，GPA分布，およびワイブル分布のては信頼度が低い．

3 種類を取り上げている．それぞれの分布の形状母数

さらに，GEVを母分布関数として抽出された年最大は表-1 に示すような 4 〜5 通りのものを標準として用値の標本データについてその平均値等も記入してある．

いたが，ときに応じてこれ以外の形状母数の分布につ

各項で POT の列が母分布関数の統計値，AM の列が年いても数値シミュレーションを行っている．

最大値の統計値である．表から分かるように，年最大5

表-2 GEV 分布を母分布関数とする極大値ならびに年最大値の母数値と統計量(H50 = 10.0 m, 裾長度 γ50 = 1.2)λ10µ (x) (m)B (m)σ (x) (m)

β2

β1

(τ3)

POT

AM

POT

AM

POT

AM

POT

AM

POT

AM

POT

0.2723

7.122

7.300

7.330

7.653

0.8875

1.324

207.1

140.9

***

***

-0.30

0.49400.37790.15200.28880.4517

6.639

6.951

6.926

7.366

0.7029

1.008

10.41

9.393

843.6

587.0
16.38

-0.45
3

A (m)

L-skw

k

AM

-0.15

0.27000.54590.7476

5.926

6.473

6.336

6.978

0.8947

1.118

2.530

2.217

19.27

0.001

0.16931.0312

1.2373

4.849

5.810

5.444

6.428

1.321

1.453

1.133

0.824

5.366

4.580

0.15

0.0771

1.922

2.025

3.232

4.904

4.090

5.645

2.102

2.060

0.436

0.074

3.140

3.047

-0.45

0.49400.37790.08870.2586

6.968

7.312

7.089

7.657

0.45881.2495

160.6

78.1

***

***

-0.30

0.20200.4219

6.330

6.976

6.531

7.381

0.49130.9879

9.708

8.570

650.8

363.1
18.14

-0.15

0.27000.45790.8240

5.330

6.580

5.643

7.010

0.74931.0752

2.503

2.368

18.97

0.001

0.1693

1.039

1.114

3.567

5.891

4.165

6.496

1.330

1.375

1.133

1.004

5.370

5.068

0.15

0.0771

2.319

1.789

0.628

5.041

1.663

5.777

2.536

1.896

0.436

0.270

3.127

3.155

表-3 GPA 分布を母分布関数とする極大値ならびに年最大値の母数値と統計量(H50 = 10.0 m, 裾長度 γ50 = 1.2)
λ

3

10

A (m)

L-skw

k

µ (x) (m)B (m)σ (x) (m)

β2

β1

(τ3)

POT

AM

POT

AM

POT

AM

POT

AM

POT

AM

POT

-0.20

0.42860.44460.5122

6.167

6.751

6.723

7.202

0.71750.9998

4.595

3.635

60.82

40.43

-0.10

0.37930.67930.6924

5.582

6.421

6.336

6.937

0.8438

1.104

2.812

2.068

17.77

11.74

AM

0.001

0.33291.03990.9368

4.802

6.012

5.841

6.609

1.0378

1.281

1.995

1.311

8.956

5.787

0.10

0.2903

1.5753

1.2572

3.792

5.521

5.224

6.210

1.3073

1.5316

1.517

0.845

5.778

3.748

0.15

0.26981.9414

1.4574

3.161

5.233

4.850

5.974

1.4806

1.692

1.337

0.662

4.868

3.207

-0.2

0.4286

0.34950.5536

5.692

6.173

6.128

7.168

0.56391.0123

4.648

3.467

63.1

38.1
10.89

-0.1

0.37930.60220.7578

4.811

6.370

5.480

6.890

0.7481

1.131

2.811

1.909

17.77

0.001

0.33291.0412

1.0383

3.550

5.943

4.590

6.542

1.039

1.330

1.994

1.135

8.953

5.372

0.10

0.2903

1.77691.4107

1.776

5.429

3.391

6.116

1.475

1.615

1.517

0.639

5.776

3.571

0.15

0.26982.32561.6456

0.600

5.126

2.622

5.862

1.774

1.800

1.337

0.438

4.867

3.137

表-4 ワイブル 分布を母分布関数とする極大値ならびに年最大値の母数値と統計量(H50 = 10.0 m, 裾長度 γ50 = 1.2)
λ

3

10

A (m)

L-skw

k

µ (x) (m)B (m)σ (x) (m)

β2

β1

(τ3)

POT

AM

POT

AM

POT

AM

POT

AM

POT

AM

POT

AM

0.75

0.45050.48180.6010

5,869

6.481

6.443

6.982

0.7761

1.0515

3.121

2.129

19.00

10.81

1.0

0.33331.0356

0.9341

4.811

6.016

5.847

6.611

1.0356

1.2784

2.000

1.315

9.001

5.819

1.4

0.21502.1805

1.3981

3.106

5.497

5.093

6.178

1.43841.6002

1.198

0.681

4.838

3.695

2.0

0.11404.22791.9221

0.536

5.015

4.283

5.747

1.9586

1.9779

0.631

0.194

3.244

3.029

0.75

0.45050.44280.7708

4.941

6.298

5.468

6.832

0.7132

1.1450

3.122

1.712

19.00

8.403

1.0

0.33331.0356

1.139

3.564

6.547

4.600

6.547

1.03581.3284

2.000

1.141

9.001

5.406

1.4

0.21502.3453

1.341

1.352

5.595

3.490

6.249

1.5474

1.5468

1.198

0.699

4.839

4.002

表-5 母集団の確率密度関数，平均値，標準偏差，および歪み度
分布関数

平均値µ(x)

確率密度関数f(x)1 / k −1

一般化極値
(GEV)
一般化パ
レート
(GPA)

[1 − k ( x − B) / A]× exp{−[1 − k ( x − B) / A]

1/ k

/A

: B ≤ x, k ≠ 0⎡Γ(1 + 2k ) ⎤A⎢ 2⎥⎣− Γ (1 + k )⎦/k

B + A /(1 + k )A(1 + k )(1 + 2k )1 / 2

B + A Γ(1 + 1 / k )⎡Γ(1 + 2 / k ) ⎤A⎢ 2⎥⎣− Γ (1 + 1 / k )⎦k / A × [( x − B) / A] k −1
ワイブル

× exp[−{( x − B) / A} k ]: B<x<∞

E[ X ] = [1 − Γ(1 + k )] / k :歪み度 β1

1/ 2⎡1⎤
B + A⎢⎥/k⎣− Γ(1 + k )⎦: k≠0

[1 − k ( x − B) / A]1 / k −1 / A

標準偏差σ(x)32
A 3 ⎛⎜ E[ X ] − 3E[ X ] ⎞⎟σ 3 ⎜⎝ × E[ X ] + 2 E 3 [ X ] ⎟⎠2 (1 − k )(1 + 2k ) 1 / 2(1 + 3k )1/ 2

A 3 ⎛ Γ(1 + 3 /k ) − 3Γ(1 + 2 /k ) ⎞⎜⎟σ 3 ⎜⎝ × Γ(1 + 1 / k ) + 2Γ 3 (1 + 1 /k ) ⎟⎠E[ X 3 ] = [1 − 3Γ(1 + k ) + 3Γ(1 + 2k ) − Γ(1 + 3k )] / k 3

E[ X 2 ] = [1 − 2Γ(1 + k ) + Γ(1 + 2k )] / k 2 :6

値の標本の統計量は母分布の統計量からかなり変化す波の総個数に相当する極大値をその母分布関数から抽る．これについては，まとめて 5. で論述する．

出する．これは前報で用いたモンテカルロ法である．

表-3, 4 は極大値の母集団がGEV およびワイブル分布

そして，極大値の時系列を各年の発生個数毎に区切り，の場合である．記述項目は表-2 と同じである．ワイブ

その中での最大値を判定し，これを年最大値のデータ

ル分布で平均発生率がλ = 10 のときには，形状母数を kとして採択する．

,= 2.0 に設定すると位置母数 B が負の値となる，すなわ
ち波高の極大値が負となってしまい現実的ではない．

(4) 年最大値の統計期間と標本の総数

このため，年最大値を抽出するためのワイブル分布の年最大値の統計期間としては，K = 10 〜10,000 年にわ母集団としては，k = 2.0 のケースを計算しなかった．

たる 12 種類を対象とした．発生させる極大値の個数は

なお，表-2 〜4 中の尺度母数と位置母数の AM の列

平均発生率を λ として，n = λ K* がその平均値である．

の値は，統計期間が K = 10,000 年の標本 500 組に対して

ただし，前項に述べたように試行毎に若干の変動があ当てはめた一般化極値 (GEV) 分布の母数の平均値である．K = 10,000 年，λ = 10 個 / 年の場合には平均で n =る．5(2) で述べるように，平均発生率が λ = 3 のときに99,995 個である．

は GEV への適合率が 0% であるけれども，参考のため

年最大値の極値統計解析では，N = 50,000 組の標本のにここに記載したものである．

作製を標準とした．ただし，統計年数が 100 年を超える場合には計算時間を短縮するため，N = 5,000,000 / K(3) 高波の平均発生率と年最大値の抽出方法に設定した．すなわち，K = 10,000 年に対する年最大値

今回の数値実験では，平均発生率はλ = 3 および 10 個の標本は N = 500 個である．

/ 年の 2 通りに設定した．前者は台風等による高波，後者は冬季季節風による高波などを想定したものである． (5) 年最大値の標本の解析方法
年最大値の抽出は次のようにして行った．

所定の母集団から抽出された年最大値の標本は，ままず，平均発生率を入力として式 (1) のポアソン分布に基づいて 1 年間の発生個数毎の出現確率を計算し，その累積分布確率をあらかじめ求めておく．確率を計ず個々のデータをすべて集積して波高区分 0.25 m を級とする度数分布を作成して確率密度分布のデータとした．また，データの 1 〜4 乗値の値を集計し，年最大値算する発生個数の上限値は経験式で設定した．この累の平均，標準偏差，歪み度，尖鋭度を計算した．極大

積分布確率は 0 から 1 まで段階的に増加する関数であり，その計算事例は既に図-1, 2 に示したところである．
次に，ある年に対応してコンピュータで疑似一様乱値の統計量についても同様である．表-2 〜4 に記載し,た統計量はこのようにして求めたものである．
年最大値の個別の標本については L-モーメント法を数を発生させ，その乱数値がポアソン分布の累積分布

適用し，GEV, GPA およびワイブル分布のそれぞれに確率を超える一つ前の発生個数を見つけだす．これが当てはめて形状・尺度・位置の 3 母数を推計した．こ

その年に割り当てられる高波の極大値の発生個数であの際には標本のデータを昇順に並べ替え，式 (24) のプる．このようにして年最大値の統計期間の各年についロッティングポジションの定数としては，前報で推奨て極大値の発生個数を割り当てる．

された α = 0.45, β = 0 の値を用いた．

また，乱数の値によっては発生個数 0 のポアソン確率 P(0) よりも小さいことが起きる．その年は閾値を超Pi =える高波が発生しなかったのであり，その年は極値統

i −α
n+β

(24)

計解析の対象から除外する．たとえば， λ = 3 に対して

ここに，i は昇順に並べ替えられた極値データの順位，

は P(0) = 0.04976 であるので，K = 100 年のシミュレーn は年最大値のデータの個数，すなわち平均的に n = λションを行うと，平均して 4.98 年は高波が発生しない

K* である．

*

3 種類の分布関数の母数値が推定されると，昇順のことになる．すなわち，有効統計期間は K = 95.024 年
として統計解析を行う．

順位毎の確率統計量が計算される．その計算値と標本全体の極大値の数は，平均的には有効統計期間 K* 年のデータとの間の相関係数を算定し，その相関係数がに平均発生率を乗じた数であるが，乱数の出現の度合

もっとも 1 に近い分布を最適合関数と判定した．すな

いによって K 年の試行ごとに若干の変動がある．

わち，相関係数の最小残差法を採用した．この方法は，“Probability Plot Correlation Method (PPC)” とも呼ばれて各年の発生個数が設定されたところで，K 年間の高7

いる（Chowdhury et al. 1991）．そして，N 組の標本に1

10Parent: GEV

λ = 10, γ50 = 1.2

ついて最適合と判定された回数を分布関数毎に数え，010

各分布関数の適合率を計算した．
Probability density, f (H)また，年最大値の標本に当てはめられた分布関数の母数に基づいて，再現期間 R = 1.01 〜1000 年に対する確率波高を推計し，母分布関数の R 年確率波高との差を調べるとともに，確率波高の変動係数を算定した．
(6) 推計値の偏り率 (Bias) の定義
本研究では，年最大値の標本に基づく各種の推計値

10

-1

10

-2

10

-3

10

-4

10

-5

が標本の抽出母体である極大値の母集団からどれだけ偏っているかを調べるのが主目的である．このため，

k =-0.30 (POT)k =-0.15 (POT)k =0.001 (POT)k = 0.15 (POT)k =-0.30 (AM)k =-0.15 (AM)k =0.001 (AM)k = 0.15 (AM)0偏りおよび偏り率の算定では，極大値分布の母集団の5

10
Wave height, H (m)偏り量の検討に際しては，まず標本ごとの各種統計

1

10

Parent: GPA

量の推計値 ξi について，その平均 µ[ξ] ならびに標準偏λ = 10, γ50 = 1.2

差 σ[ξ] を式 (20), (21) で計算する．その際の標本個数 N

010は，本研究では N = max{50,000 , 5,000,000/K} としていk =-0.20 (POT)k =-0.10 (POT)k = 0.15 (POT)k =-0.20 (AM)k =-0.10 (AM)k = 0.15 (AM)Probability density, f (H),る．
偏りの度合いは絶対値ならびに相対比の二つの方法で定義することができる．すなわち，(25)

Relative Bias = ( µ[ξ ] − ξ 0 ) / ξ 0

(26)

20

(a) 母集団が一般化極値分布の場合

値を偏り判定の基準の真値として扱う．

Absolute Bias = µ[ξ ] − ξ 0

15

-1

10

-2

10

-3

10

-4

10

ここに，ξ0 は極大値の母集団の値である．式 (25) の偏りはここでは偏り量と呼び，これは，分布関数の形

-5

10

0

状母数などその値が 0 となる可能性のあるものに対して用いる．式 (26) の偏りの相対比はここでは偏り率と5

10Wave height, H (m)1520(b) 母集団が一般化パレート分布の場合呼ぶことにする．なお，表示に際しては実数値ではな
く，百分率で表すこともある．

10

５．年最大値の統計値の特性

101Parent: Weibull
λ = 3, γ50 = 1.2
0Probability density, f (H)(1) 極大値および年最大値の出現確率分布
最初に，母集団である極大値とそれから抽出された
年最大値のデータの出現確率分布を図-4 に示す．極大値については，表 - 5 の確率密度関数で計算される値を示し，年最大値については K = 10,000 年の標本 500 組のデータについて級幅 ∆H = 0.25 m で集計した度数分布か

k = 0.75(POT)k = 1.0 (POT)k = 1.4 (POT)k = 2.0 (POT)k = 0.75(AM)k = 1.0 (AM)k = 1.4 (AM)k = 2.0 (AM)

-1

10

-2

10

-3

10

-4

10

ら換算したものである．データは形状母数毎に分けて表示しており，極大値のデータは POT，年最大値の

-5

10

0

データは AM の凡例で区別している．
図-4 で明らかなように，年最大値のデータは全体と5

10
Wave height, H (m)1520(c) 母集団がワイブル分布の場合

して波高の大きなほうへシフトし，ピークの尖り具合

図-4 母集団ならびに年最大値のデータの出現確率分布も緩和されている．このため，年最大値の平均値は極8

大値の平均よりも常に大きくなり，歪み度は逆に小さ

やや低い段階で年最大値の標本が一般化極値分布に適くなる．表-2 〜4 には，こうした統計値が記載されて合する．母集団がガンベル分布のときには，一般化極いるので，POT と AM の列の数値を比較することで，値分布への移行が λ = 3.1 で完了する．

そうした変化を吟味できる．標準偏差は，形状母数が小さいときは年最大値のほうが極大値よりもかなり大(3) 年最大値に対する一般化極値分布の形状母数

きな値であるけれども，形状母数が大きくなるにつれ母集団が一般化パレート分布である極大値から抽出てその差が縮まる．母集団が一般化極値分布の場合にされた年最大値のデータは，標本が十分に大きければは，k > 0 で年最大値の標準偏差が母集団よりも小さく

式 (13) で表されるような一般化極値分布に従い，形状なっている．

母数は母集団と同一であり，尺度母数と位置母数は式(14) に従って変化する．

(2) 年最大値データの一般化極値分布への適合度

今回の数値実験の結果から，平均発生率 λ = 10，統年最大値の標本が十分に大きな場合，すなわち非常

計期間 K = 10,000 年のデータについて母集団と年最大に長い期間にわたって年最大値のデータが収集された値に当てはめた一般化極値分布の形状母数を比較した

場合には，年最大値のデータが一般化極値分布に従うのが図-6 である．形状母数は表-1 に記載したものに加ことが今回の数値実験でも確かめられた．すなわち，

えてその他のケースについても計算してある．なお，

K = 10,000 年の標本 500 組について 3 種類の分布関数の図の横軸には L-モーメント法で計算される L-skewness

適合度を調べたところ，母分布関数がどれであってもの値を用いており，分布関数毎のこの値は表-2 〜4 に

一般化極値分布への適合率が 100% であった．ただし， 記載されている．
平均発生率が低いときには一般化パレート分布への適母集団が一般化パレート分布の場合には，予測通り合率が高く，平均発生率が高くなるとワイブル分布へに一般化極値分布の形状母数は母集団の値と同じであよく適合し，一般化極値分布が最適合分布となるのはる．母集団が一般化極値分布であっても年最大値の平均発生率が 5 程度以上の場合に限られる．

データは一般化極値分布で近似され，そのときの形状
母数は母集団よりもやや大きくなる．図-6(a) でのデー

100

タに直線回帰式を当てはめると，次のような関係が得られる．

Return rate (%)

80

Parent: GPA(k = 0.10)60kAM = 1.064kPOT + 0.036

Best Fit to:40

一方，母集団がワイブル分布の場合には年最大値分

GEV
GPA
WBL

20

布が式 (8) のように導かれており，一般化極値分布に収束することが証明されているわけではない．しかし，数値実験の結果では統計期間が十分に長いときには一

0
2.0

(27a)

2.5

3.0

3.5
4.0
Mean rate, λ

4.5

5.0

般化極値分布が最適合分布となっている．その際に推計された形状母数 kAM は母集団の形状母数 kPOT と図-

図-5 極値分布関数の適合率に及ぼす平均発生率の影響

6(b) に示すような関係になり，これを 2次多項式で近似すると次のように表される．

図-5 は，母集団が一般化パレート分布で形状母数がkAM = –0.230kPOT2 + 0.806kPOT – 0.578

k = 0.10 の場合について，平均発生率を λ = 2 〜5 に漸変(27b)

させてみた結果である．λ ≤ 2 では年最大値が一般化パ(3) 年最大値に対する一般化極値分布の尺度母数とレート分布によく適合し，2.6 ≤ λ ≤ 3.3 ではワイブル分
布の適合率が高く，λ = 4 になってはじめて一般化極値

位置母数

分布の適合率が 100% となる．

母集団が一般化パレート分布である極大値から抽出

母集団が指数分布（ワイブルで k = 1.0）のときには， された年最大値に当てはめられた一般化極値分布の尺
年最大値への最適合分布が一般化パレート分布からワ度母数は，式 (14) によって母集団の尺度母数に関係づ

イブル分布へ移行するのが λ = 2.75 であり，さらに一般けられる．図-7 はこの関係を確認するもので，一般化

化極値分布へ移行するのがλ ≥ 3.25 である．すなわち，

パレート分布のデータは式 (14) の理論値に一致してい母集団が一般化レート分布のときよりも平均発生率がる．また，ワイブル分布については理論式がないもの

93.0
Ratio of scale parameter, AGEV / APOT

GEV-fitted shape parameter, kGEV

0.4

Record of K = 10,000 years
Mean Rate, λ = 10
0.2

0

-0.2
Parent Distribution
GEV
GPA
Linear

-0.4

-0.6
-0.5

Record of K = 10,000 years
Mean Rate, λ = 10

2.5

GPA theory
Parent: GEVParent: GPA
Parent: WBL

2.0
1.5
1.0
0.5
0
0

-0.4

-0.3

-0.2

-0.1

0

0.1

0.2

0.1

0.2

0.3

Parent L-skewness,

0.3

Parent shape parameter, kPOT

0.4

0.5

τ3

図-7 年最大値に当てはめた一般化極値分布の尺度母数と(a) 母集団が一般化極値分布および母集団の尺度母数の比

一般化パレート分布の場合
BGEV / B*

1.4

0.1

Record of K = 10,000 years
Mean Rate, λ = 10

Ratio of location parameter,

GEV-fit shape parameter, kGEV

0.15

0.05
0
-0.05

Parent Distribution
WBL
Empirical-0.1

Record of K = 10,000 years
Mean Rate, λ = 10
1.3
Parent: GEV
1.2

1.1

1.0

0.9
0

-0.15
0.7

0.1

0.2

0.3

Parent L-skewness,

0.9

1.1

1.3

1.5

0.5

図-8 母集団が一般化極値分布の場合の年最大値分布とParent shape parameter, kPOT

極大値分布の位置母数の比(b) 母集団がワイブル分布の場合
図-6 母集団と年最大値に対する一般化極値分布(10) に示したような関係にある．数値実験の結果を調の形状母数の関係

べると，ワイブル分布を母集団とする年最大値に当てはめた一般化極値分布の位置母数についても，理由はの，結果として式 (14) の値でほぼ近似される．

不明であるものの，式 (10) の関係が成立することが確

一方，母集団が一般化極値分布の場合には，年最大認された．

値に当てはめられた一般化極値分布の尺度母数がかな

母集団が一般化極値分布の場合には，式 (10), (14) のり大きな値に変化する．この尺度母数の変化率を経験ような変換式が見いだせない．数値実験で当てはめら式で近似すると次のように表される．

AAM/APOT = 3.668kPOT2 – 2.402kPOT + 1.058

0.4

τ3

れた位置母数 BAM を式 (10) による変換値 B* に対する比(28)の形で表示したのが図-8 である．

また位置母数に関しては，母集団が一般化パレート(4) 大標本の年最大値データによる確率波高推計値の分布のときは式 (14) で変換式が与えられており，年最偏り率

大値の数値実験の結果でもこの式で変換したとおりになることが確認されている．ただし，平均発生率が λ =前項では，年最大値の標本が十分に大きくて平均発

3 のときには図-5 に示したように年最大値が一般化極

生率が高ければ，母集団の分布特性に関わりなく年最値分布に適合しないため，式 (14) が適用できない．ま

大値のデータが一般化極値分布に適合することを紹介た，ワイブル分布に含まれる指数分布を母集団とするした．また，平均発生率が低いときには図-5 のように，

年最大値は，式 (9) のようにガンベル分布に収束し，大標本の年最大値データが一般化パレート分布あるい*そのときの位置母数 B は指数分布の位置母数 B と式はワイブル分布に適合することを示した．しかしなが10

ら，これは当てはめ候補として一般化極値分布，一般

100 年と 1000 年確率波高の偏り率はそれぞれ –5% と –化パレート分布，およびワイブル分布の 3 種類に限定

9% とやや小さくなる．

したためであり，適合率が 100% であっても当てはめ

実際の極値統計解析では統計期間が 1 万年のようなられた分布関数から推計される確率波高が母集団の値

大標本の年最大値データが取得できないので，確率波に一致するわけではない．

高の偏り率については 7. を参照して頂きたい．

例外は母集団が一般化パレート分布の場合であって，６．年最大値の統計分布に及ぼす標本サイズの影響

このときには式 (13) のように年最大値データが一般化
極値分布となる．今回の数値シミュレーション結果を調べると，当てはめられた一般化極値分布から推計される再現期間 R = 20 〜1000 年に対する確率波高の値は，母集団である一般化パレート分布の値からの偏り率が(1) 年最大値の統計期間による極値分布適合率の変化
前節においては年最大値の標本が十分い大きい場合–0.3% 〜+0.04% 異なるだけである．

として，年最大値が K = 10,000 年にわたって収集でき母集団がワイブル分布の場合には，形状母数が k =たという仮想的なケースを対象として議論をしてきた．

1.0，すなわち指数分布の場合だけ年最大値データが一

しかし，自然現象を対象とする実際の極値統計解析で

般化極値分布となる．これは式 (9), (10) で示したとこ

はそのような長期間のデータを取得することは不可能

ろであり，確率波高の偏り率もほぼ 0 である．しかし， であり，数十年程度のデータを用いて極値統計解析を
形状母数が k < 1 のときには当てはめられた一般化極値

行い，100 年確率事象などの予測を行うのが現実であ分布から推計した確率波高がやや過大評価となる．母る．

集団の形状母数が k = 0.75 であれば，1000 年確率波高の一方，数値実験では統計期間の長さを変える，すな

推計値が +4.6% の偏り率を持つ．逆に形状母数が k > 1

わち標本の大きさを変えることには何も差し支えがなであれば確率波高が過小評価となる．母集団の形状母い．そこで，統計期間を K = 10 〜10,000 年の範囲で 12数が k = 1.5 であれば，1000 年確率波高の推計値が –段階に変え，3 種類の極値分布関数への適合率を調べ

4.4% の偏り率を持つ．

てみた．図-9 は形状母数が k = 1.0 のワイブル分布（指数分布）で平均発生率が λ = 10 の母集団から年最大値

一般化極値分布を母集団として持つ極大値デーから
年最大値を抽出した場合には，そのデータに当てはめを抽出した場合の計算結果である．一般化極値分布

た一般化極値分布から推計した確率波高が負の偏り率(GEV)，一般化パレート分布 (GPA)，およびワイブルを示す．母集団の形状母数が k = –0.45 および 0.20 のと分布 (WBL) への最適合判定は，4 (5) に述べた相関係

きの 1000 年確率波高の偏り率がそれぞれ –1.2% と –数の最小残差法に基づいている．

3.5% である．100 年確率波高であればそれぞれ –0.2%

この図-9 の結果で明らかなように，統計期間が 100と –1.6% である．

年程度以下では 3 種類の極値分布に対する当てはめ適

図-5 に示したのは母集団が一般化パレート分布 (k =合率が拮抗しており，最適合となる分布は個別の標本

1.0) の場合である．統計期間が 1 万年という大標本でによって変動する．統計期間が 100 年程度以下では一あっても平均発生率が λ = 2.5 のときには一般化パレー

般化極値分布が適合する割合が高いけれども，統計期

ト分布が最適合となる．このときに当てはめられて一間が 15 年以下では一般化パレート分布も同程度の割合

般化パレート分布から推計された 100 年および 1000 年で適合する．ワイブル分布への適合率は標本が大きく

確率波高はそれぞれ –5.8% と –14.8% の偏り率をもつ．

なるにつれて増大し， 30%〜35% の割合にまで達する．

平均発生率が λ = 3.0 のときにはワイブル分布が最適合

他の分布に紛れることなしに一般化極値分布が当てはとなり，この場合の 100 年と 1000 年確率波高の偏り率

められるのは，統計期間が 5000 年を超える場合である．

れ –1.2% と小さい．

母集団が指数分布の場合に加えて，k = 0.001 の一般

図には示していないが，母集団が一般化極値分布や化極値分布と k = 0.10 の一般化パレート分布の母集団か

ワイブル分布の場合も状況は同様である．いずれも平ら抽出された年最大値のデータが一般化極値分布に最均発生率が低い場合には一般化パレート分布またはワ適合となる割合を調べたものが図-10 である．母集団

イブル分布が最適合となる．一般化パレート分布が当が一般化極値分布の場合には，年最大値のデータが一てはめられたときには，推計された確率波高の偏り率

般化極値分布に適合する割合が比較的に高い．とはが最大で –22% と大きなものとなる．ワイブル分布が言っても，確実に一般化極値分布に適合するようにな当てはめられたときも負の偏り率を生じるけれども，るのは，統計期間が 1000 年を超えるような場合である．

11

Mean (k = -0.30)
Mean (k = 0.001)

100

Skewness ratio, µ (β1) / (β1)0 & σ (β1) / (β1)0Parent: Weibull(k = 1.0, λ = 10)Best-fit rate (%)

80
Best Fit to
60

GEV
GPA
WBL

40

20

0
10

100

1000

10000

St. Dev. (k = -0.30)
St. Dev. (k = 0.001)

1.0

Parent: GEV Distribution(λ = 10)

0.8

0.6

0.4

0.2

0

Duration of record, K (years)

10

図-9 指数分布から抽出された年最大値のデータが

100

1000

10000

Duration record, K (years)

3 種類の極値分布関数に適合する割合(a) 母集団が一般化極値分布の場合

100

Mean Rate, λ = 10

Mean (k = -0.20)
Mean (k = 0.001)

St. Dev. (k = -0.20)
St. Dev. (k = 0.001)

1.0
Skewness ratio, µ (β1) / (β1)0 & σ (β1) / (β1)0Best-fit rate to GEV (%)

90
80
70
60
50
Parent: GEV (k = 0.001)
Parent: GPA (k = 0.10)
Parent: WBL (k = 1.0)

40
30
10

100
1000
Duration of record, K (year)

10000

0.8

Parent: GPA Distribution(λ = 10)

0.6

0.4

0.2

0
10

図-10 3 種類の母分布関数から抽出された年最大値の

100

1000

10000

Duration record, K (years)

データが一般化極値分布関数に適合する割合(b) 母集団が一般化パレート分布の場合

Skewness ratio, µ (β1) / (β1)0 & σ (β1) / (β1)0k = 0.10 の一般化パレート分布の場合には指数分布と同じような挙動である．
(2) 年最大値の統計期間による標本の歪み度の変化
このように年最大値の標本が小さいときに一般化極値分布に収束しない基本的な原因は，母集団から年最大値の分布関数を導く式 (3) が確率の期待値としての計算であるためである．計算の前提としているポアソン分布の式 (1) はデータの数が十分に大きいときのもの

1.0

0.8

Parent: WBL Distribution(λ = 10)0.6

Mean (k = 0.75)
Mean (k = 1.40)
St. Dev. (k = 0.75)
St. Dev. (k = 1.40)

0.4

0.2

0
10

であり，データの数が少なければ既に図-1, 2 に示した
ように，標本毎にポアソン確率とは異なる値を取る．

100
1000
Duration record, K (years)

10000(c) 母集団がワイブル分布の場合

標本毎の統計的変動は，母集団の分布関数の計算に直接に影響する．

図-11 年最大値の標本の歪み度の平均と標準偏差

標本毎の統計的変動の影響を直接に示すのは，標本が小さいときにデータの歪み度が小さめに偏って現れ間によって変化する状況を示したものである．3 種類

る現象である（付録参照）．図-11 は，表-2〜4に記載の母集団について 2 通りの形状母数の値の場合を例示

した年最大値の標本全体の歪み度 (β1)0 を基準として，している．一般に，形状母数が小さいほど歪み度の低標本の歪み度の平均値と標準偏差が年最大値の統計期下が著しい．また，統計期間が短い場合には標準偏差12Difference of shape parameter, ∆k = kAM - kPOT

Sample L-skewness, t 3

0.20

Mean Rate, λ = 10

0.15

0.10

Parent Distribution
0.05

GEV (k = 0.001)
GPA (k = 0.10)
WBL (k = 1.0)

0
10

100
Duration of record,

1000

10000

0.5

Parent: GEV (γ50 = 1.2)

0.4

k =-0.30 (λ = 3)k =0.001 (λ = 3)
k =-0.30 (λ =10)k =0.001 (λ =10)

0.3

0.2

0.1

0
10

20

K (years)

100

50

Duration of record, K (years)(a) 母集団が一般化極値分布の場合

図-12 統計期間による年最大値の標本 L-skewness の変化

Parent Distribution
GEV (k = 0.001)
GPA (k = 0.10)
WBL (k = 1.0)

0.15

Difference of shape parameter,

GEV-fitted shape parameter,

kGEV

Mean Rate, λ = 10

∆k = kAM - kPOT

0.20

0.10

0.05

0
10

100
Duration of record,

1000

10000

K (years)

0.9

Parent: GPA (γ50 = 1.2)

0.8

k =-0.20 ( λ = 3)
k = 0.15 ( λ = 3)
k =-0.20 ( λ =10)k = 0.15 ( λ =10)

0.7
0.6
0.5
0.4
0.3
0.2
10

20

50

Duration of record,

図-13 統計期間による年最大値の標本に当てはめた

100

K (years)(b) 母集団が一般化パレート分布の場合

が平均に比べてかなり大きく，ここに例示したもので

4.0

は統計期間が 10 年の場合など，歪み度の標準偏差が平kAM / kPOT

3.5

Shape parameter ratio,

一般化極値分布の形状母数の変化

2.5

均値よりも大きくなっている．すなわち，標本毎の変動が極めて大きい．

L- モーメント法では，標本に対して算定された Lskewness の値が当てはめられる分布関数の特性を規定し，このパラメータは標本の歪み度と密接に関連している．したがって，標本の歪み度が年最大値のデータ
全体の統 計量から低 落すると， それに応じ て L-

Parent: Weibull (γ50 = 1.2)

3.0

2.0
1.5
k = 0.75 ( λ = 3)
k = 1.40 ( λ = 3)
k = 0.75 ( λ =10)k = 1.40 ( λ =10)

1.0
0.5
0

20

10

skewness の値も低減し，本来の極値分布に当てはめらDuration of record,

れる割合が低下するのである．

50
K (years)

100(c) 母集団がワイブル分布の場合

図-12 は標本について計算した L-skewness の値が標本の統計期間に変化する状況を例示したものである．こ図-14 年最大値の標本を母集団の極値分布に当てはめたとき

うした L-skewness の値の変化に応じ，年最大値の標本の形状母数の値に当てはめた一般化極値分布の形状母数の値は，統計年数が短いほど大きめの値となる．図-13 は，図-12 に当てはめて推計した形状母数の値が増加する状況が明例示したケースについて形状母数 kGEV の値の変化を示らかである．

したものである．統計期間が十分に長い K = 10,000 年における形状母数の値は，既に図-6 に示したところで(3) 統計期間が短い標本の年最大値分布の形状母数あり，統計期間が短くなるにつれて一般化極値分布に

先に図-9 に例示したように，統計期間が数十年程度13である極値統計解析の実務においては，年最大値の標16本に対して一般化極値分布を一律に当てはめるのは標Return wave height, HR (m)本毎の特性を考慮しない方法であり，推奨すべきでないであろう．実際問題では，標本に最も適合する極値分布を選び，推計された形状・尺度・位置母数の値を用いて R 年確率統計量を推計することになる．本項では，母集団と同じ極値分布に当てはめられたときに形状母数の値がどの程度変化するかを検討する．

Parent: GEV
k = -0.15, λ = 10

POT (λ = 10)
AM (K = 10 years)
AM (K = 25 years)
AM (K =100 years)

14
12
10
8
6

図-14 は統計期間が K = 10 〜100 年の年最大値標本に4ついて，母集団と同じ分布関数を当てはめたときの形

0.1

1状母数を母集団の形状母数と比べたものである．図-

10

100

1000

Retrun period, R (years)

14(a), (b) は，一般化極値分布と一般化パレート分布の(a) 母集団が一般化極値分布の例

場合であり，これらの分布は形状母数が 0 を挟んで正16負の範囲に広がるので，形状母数の値の差を表示してReturn wave height, HR (m)いる．ワイブル分布の場合には形状母数が常に正の値なので，図-14(c) のように母集団の形状母数に対する比として表示している．
図で明らかなように，統計期間が短くなるにつれて形状母数の推計値が大きくなっている．また，母集団の形状母数の影響は極値分布によって異なり，一般化
極値分布とワイブル分布では母集団の形状母数が小さ
いほうが推計値の変化が大きくなっている．ただし，Parent: GPA
k = -0.2, λ = 10

14
POT (λ = 10)
AM (K = 10 years)
AM (K = 25 years)
AM (K =100 years)

12
10
8
6
4
0.1

一般化パレート分布では統計期間によってその関係が逆転している．平均発生率の影響は比較的に小さい．

1

10
100
Retrun period, R (years)

1000(b) 母集団が一般化パレート分布の例

図-14 では統計期間が 100 年までしか図示していない
が，形状母数の値の差異は統計期間が長くなっても解14

消するわけではない．統計期間が K = 10,000 年の標本

Parent: Weibull
k = 2.0, λ = 3

12Return wave height, HR (m)の計算結果を調べると，一般化極値分布では ∆k = 0.02
〜0.05 程度の差が残る．一般化パレート分布では，∆k= 0.27 〜0.56 とかなり大きな差である．ワイブル分布では kAM/k POT = 1.4 〜1.7 である．したがって，たとえば年最大値の標本がワイブル分布に適合したとしても，そのときに推計される形状母数は母集団である極大値の分布関数に対するものよりも 1.5 〜2.0 倍程度大きく

108
6POT (λ = 3)
AM (K = 10 years)
AM (K = 25 years)
AM (K =100 years)

420なっている可能性が高い．これについては，実際の波0.1

浪データとの比較を 9. で述べる．

1

10
100
Retrun period, R (years)

1000(c) 母集団がワイブル分布の例

７．年最大値資料に基づく確率波高の偏り率

図-15 母集団と年最大値の標本から推計される(1) 確率波高の推計手順

確率波高の比較

前節で述べたように年最大値の分布が数十年程度のの標本について最適合と判定された分布関数を用いて場合には，その分布関数が一般化極値分布に収束する

R 年確率波高を推計することになる．

という保証がない．図-9, 10 に示したように，一般化パ

本研究では，個々の標本ごとに最適合分布を判定し，レート分布もワイブル分布もかなりの割合で最適合分その分布に当てはめられた形状・尺度・位置母数の数

布として判定される．極値統計解析の実務では，個別値を用いて，R = 1,01, 2, …, 500, 1000 年に対する確率14

5

Parent: GEV
k = 0.001, λ = 10

Bias rate of return wave height (%)

Bias rate of return wave height (%)

5

0

-5

-10

K = 10 years
K = 15 years
K = 25 years
K = 40 years
K = 60 years
K =100 years

-15

-20

Parent: GEV (λ = 10)
Annual Max of 40 years

0

-5

-10
GEV k = -0.3
GEV k = -0.15
GEV k = 0.15
-15

1

10

100

1

1000

(a) 母集団が一般化極値分布の例
5

Parent: GPA
k = 0.15, λ = 10

0

-5

-10

K = 10 years
K = 15 years
K = 25 years
K = 40 years
K = 60 years
K =100 years

-15

-20

1000

Parent: GPA (λ = 10)
Annual Max of 40 years

0

-5

-10

GPA k = -0.2
GPA k = -0.1
GPA k = 0.15

-15

1

10
100
Return period, R (years)

1

1000

(b) 母集団が一般化パレート分布の例

10100Return period, R (years)

Bias rate of return wave height (%)

Parent: Weibull
k = 1.0, λ = 3

0

-5

-10

-15

K = 10 years
K = 15 years
K = 25 years
K = 40 years
K = 60 years
K =100 years

-20

10
100
Return period, R (years)

Parent: Weibull ( λ = 3)
Annual Max of 40 years

0

-5

-10

WBL k = 0.75
WBL k = 1.00
WBL k = 1.40
WBL k = 2.00

-15
1

-25
1

1000

(b) 母集団が一般化パレート分布の例
55Bias rate of return wave height (%)100(a) 母集団が一般化極値分布の場合

Bias rate of return wave height (%)

Bias rate of return wave height (%)

5

10Return period, R (years)

Return period, R (years)

10
100
Return period, R (years)

1000(c) 母集団がワイブル分布の例

1000

図-17 統計期間 40 年の年最大値資料から推計される(c) 母集団がワイブル分布の例

確率波高の偏り率

図-16 年最大値の標本から推計される確率波高の偏り率の比は変動係数であり，これについては 8. で述べる．
波高を推計した．そして，N 組の標本から得られた再3 種類の母数団から抽出された年最大値の標本につ

現期間毎の確率波高の平均値をもって，N 組全体に対いて確率波高を推計し，母集団における確率波高と比する確率波高と定義した．また，N 組の標本に対する較した事例を，母集団毎に図-15 に示す．年最大値とい確率波高の標準偏差を計算した．標準偏差と平均値とうのは，1 年間に出現する極値の中の最大値であるか

15

ら，再現期間としては R < 1 年はあり得ない．したがっえた事例を示している．母集団の分布関数によって若て，年最大値の確率波高は R = 1 年から立ち上がる形と干の差はあるけれども，50 年確率波高で 2% 程度，100

なる．図-15 には母集団毎に統計年数がK = 10, 25, およ

年確率波高で 3% 程度の過小評価が見込まれる．

び 100 年の場合に推計される確率波高を示しており，高波の極値統計では再現期間として 100 年を超える

母集団における確率波高と比較できるようにしてある． 事象を推計する例はあまりないけれども，高潮などの図から明らかなように，再現期間 10 年付近では年最大場合にはオランダで行われているように 1 万年確率偏値と母集団との確率波高に差異はないけれども，再現差を求めることもあり得る．高潮などについては裾長期間が短い場合ならびに長い場合には年最大値が確率度の設定など高波と異なる点もあるので，図-16 で示さ

波高を過小評価している．再現期間が 5 年未満についれているトレンドをそのまま引き延ばすわけにはいかては年最大値という定義に起因するものであり，是正いであろう．それでも，年最大値資料が確率統計量をは不可能である．

過小評価する可能性が大きいことを考えるならば，高しかし，再現期間 50 年以上に対する過小評価は，年潮についても極大値資料を使用することが推奨される．

最大値の標本がそのサイズの影響を受けた結果である．

８．年最大値資料に基づく確率波高の変動係数

すなわち，図-14 に示したように，年最大値の標本の統計期間が数十年程度の場合には形状母数の推計値が大極大値資料から推計される確率波高に関しては，前きめに見積もられる．極値分布関数は形状母数が大き報においてその変動係数の経験的推定式を導いた．本くなるにつれて確率密度関数の裾が短くなり，確率波研究においてもその推定式を参考として，年最大値資高が低めに推計されるようになる．母集団の確率波高

料から推計される確率波高の変動係数について作業をとの差は統計年数が長くなるにつれて減少するけれど行った．ただし，前報では裾長度として 3 通りの値もも，K = 100 年程度ではまだ明らかな差が認められる．

を用いて裾長度の影響を定式化したけれども，今回はγ50 = 1.20 に固定している．そこで，裾長度に関しては(2) 年最大値から推計される確率波高の偏り率

前報の定式化をそのまま使用した．

確率波高の推計値の差異を偏り率として表示した例

さらに，年最大値資料の場合には変動係数が標本の大が図-16 である．図-14 とは異なる形状母数の場合につきさ n の –1/2 乗にほぼ比例したけれども，年最大値資

いて，統計年数が 10 〜100 年の 6 通りの場合を示して料の場合には母集団によって K-0.4 〜K-0.45 に比例する形いる．母集団が一般化パレート分布の場合には再現期となった．そこで，年最大値資料の確率波高推計値の間が 200 年を超えるあたりから偏り率の変化が小さく変動係数に関する推定式を次のように定式化して，そなるけれども，母集団が一般化極値分布とワイブル分

こに用いられる定数 a, b, c を数値計算の結果から定める

布の場合には再現期間が伸びるに従って偏り率が単調ことにした．

に増加する．

CV [ H R ] = (a y R2 + b y R + c) (γ 50−1) 0.7 / K κ

高波の極値統計の場合には実際の波浪データの統計(29)

期間の上限が 40 年程度であるので，K = 40 年のケース

ここに，κ = 0.45（一般化極値分布）および κ = 0.4（一について確率波高推計値の偏り率を提示したのが図-17

般化パレート分布およびワイブル分布）と設定する．

である．3 種類の母数団毎に形状母数を 3 〜4 通りに変定数 a, b, c については表-6 のように設定した．

表-6 年最大値資料から推計された確率波高の変動係数の近似推定式の係数表
分布関数
一般化極値分布 (GEV)

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

ワイブル分布係数の推定式
a = exp(-0.00088k2 + 0.2458k + 0.01295)
b = –4.0490k2 – 1.8910k + 0.0062
c = 4.2782k2 + 2.9640k + 0.5782
a = 24.594k4 + 6.587k3 + 1.0423k2 + 0.3495k +0.0493
b = –68.204k4 – 18.949k3 – 4.1795k2 – 1.5050k + 0.0077
c = 20.605k3 + 6.878k2 + 2.3229k + 0.4161
a = 0.3385k2 – 0.4510k + 0.1280
b = –0.9265k2 + 1.3062k+ 0.3658
c = 0.5342k2 – 0.3864k + 0.2756

16

果のデータは＋および×で示してある．横軸の同じ値1.2の位置に複数のデータがあるのは，年最大値の異なる

Parent: GEV (k = 0.001)

統計期間の計算結果である．図-18(b), (c) で横軸の値数

K 0.45 x CV of HR

1.0
Prediction
Calc. (λ = 3)
Calc. (λ = 10)

0.8

が 0 の付近で数値実験の値が異常に大きいのは再現期
間 1.01 年に対するものであり，無視してよい．
なお，式 (29) の基準化変数 yR および表-6 の定数 a, b, c0.6の算定式は，抽出された極大値の母集団の形状母数を0.4

引数としている．年最大値の標本を解析するときにはその極大値の母集団は確定しておらず，したがってそ

0.2の形状母数も分からない．実際問題への適用にあたっ
0
0

2

4
6
Reduced variate, yRては，図-14 に示す年最大値の小標本における形状母8数のデータを参考とし，年最大値が抽出されたであろう極大値の母集団の形状母数を推定する作業が必要と(a) 母集団が一般化極値分布の例なる．これについては 9 (2) で計算例を示すので参照さ

1.0

れたい．

Parent: GPA (k = 0.15)

K 0.4 x CV of HR

0.8

Prediction
Calc. (λ = 3)
Calc. (λ = 10)

９．現地の高波データへの適用性について

0.6(1) 現地波高データの極大値と年最大値解析の比較
先に合田・久高 (2010) は日本海沿岸 11 地点の高波

0.4

データについて L-モーメント法による地域頻度解析を0.2行った．そのデータは極大値資料であったが，各地点について年最大値のデータを抽出し，極値解析を行っ

0
0

1

2
3
Reduced variate, yR

4てみた．欠測を除いた有効統計年数は K* = 18.0 〜37.2

5

年であった．また，ポアソン分布について検討したア
ラスカ沿岸 Kodiak 島沖の波浪データについても年最大(b) 母集団が一般化パレート分布の例

値を抽出して極値解析を行った．

1.0

これら 12 地点の極大値の標本はいずれもワイブル分

Parent: Weibull ( k = 1.4)

布によく適合し，また年最大値の標本もワイブル分布

K 0.4 x CV of HR

0.8

Prediction
Calc. (λ = 3)
Calc. (λ = 10)に適合した．そこで，ワイブル分布として推計された
形状母数の値の比率を求めた結果を図-19 に示す．

0.6

日本海沿岸の高波の極大値資料では，ワイブル分布の形状母数が k POT = 1.060 〜1.344，平均で 1.225 であり，0.4

Kodiak 沖では k POT = 1.518 であった．数値シミュレー
0.2

ションでは概略値として k POT = 1.25 を用いた．シミュレーションによる形状母数の比は，平均値を実線で表

0
0

1

2
3
Reduced variate, yR

し，その上下に標準偏差の 1 倍をとった値を一点鎖線

4および二点鎖線で示している．形状母数の比が正規分布すると仮定すると，この二つの鎖線は形状母数比の(c) 母集団がワイブル分布の例

68% 信頼区間を表すことになる．現地の高波のデータ
はばらつきが大きいものの，十分に信頼区間内には

図-18 年最大値資料に基づく確率波高の変動係数の入っており，今回の数値シミュレーションの結果が現推定式と数値計算結果の比較

地にも適用できることを示している．
これらの経験式による確率波高推計値の変動係数のまた 100 年確率波高について，極大値資料から推計

予測値を数値実験で得られたデータと比較した結果のされた値を真値 ξ0 とみなし，式 (26) を用いて年最大値

例を図-18 に示す．予測値は曲線で示し，数値実験結

資料による確率波高推計値の偏り率を計算した結果を17

図-20 に示す．各地点の偏り率は –9.9% 〜0.8% とばら

K* = 34.4 年，λ = 15.06 回/ 年 (Hc = 4.0 m)

ついており，算術平均が –4.6% である．12 地点の平均

kAM = 2.313，AAM = 3.230 m，BAM = 4.331 m

統計年数 26.8 年に対応する数値シミュレーションによH10 = 8.96 m，H50 = 10.16 m，H100 = 10.58 m，γ50 = 1.13

る 100 年確率波高の平均偏り率は –3.4% であり，現地8. で紹介した確率波高の変動係数の推定式を使うた

データよりはわずかに小さい．しかし，偏り率の 68%

めには，まず極大値の母集団に対する形状母数を推定

信頼区間の幅が非常に広いため，この程度の平均値のする必要がある．そこで図-14 を参照し，kAM / k POT = 1.8

差は統計的に有意なものではない．したがって，確率と想定する．それによって，kPOT = 2.313/1.8 = 1.29 と推

Ratio of AM to POT shape parameter, kAM/kPOT

波高の推計値についても数値シミュレーションの結果定される（実際の酒田港の極大値の形状母数は波高閾
値によって k = 1.17 〜1.27 の間で変化した）．

5

この形状母数の推定値 kPOT = 1.29 を表-6 の定数算定式

Simulation upper
Simulation mean
Simulation lower
Japan Sea
Kodiak4に代入すると，ワイブル分布の確率波高の変動係数の定数が次のように計算される．

3

a = 0.3385×1.2932 – 0.4510×1.29 + 0.1280 = 0.1095
b = 0.9265×1.292 + 1.3062×1.29 + 0.3658 = 0.5090

2

c = 0.5342×1.292 – 0.3864×1.29 + 0.2756 = 0.66611

この定数を K* = 34.4 年，γ50 = 1.13 とともに式 (29) に代入する．年最大値の場合の再現期間 R = 100 年に対す

0
10

20
50
Record duration, K (years)

る基準化変数は y100 = (ln100)^(1/2.313) = 1.935 である．

100

したがって，100 年確率波高に対する変動係数が次のように計算される．

図-19 日本海沿岸および Kodiak 沖の波高に当てはめられた
ワイブル分布の形状母数の比率

Bias of 100-year wave height (%)

15

CV = (0.1095×1.9352 + 0.5090×1.935 + 0.661)
× (1.13 – 1.00) 0.7/34.40.4 = 0.0761

Simulation upper
Simulation mean
Simulation lower
Japan Sea
Kodiak10

これによって H100 = 10.58 m に対する標準偏差が σ =0.0761×10.58 = 0.805 m と算定される．さらに，図-14(c)5を参照して 100 年確率波高に対する偏り率を –3% と見0積もると，H100 = 10.58 m に対する 68% 信頼区間がH100 =-5

10.09〜11.70 m と推定される．

-10

なお，極大値資料から推計される確率波高の変動係-15数については前報においてその推定式を次のように提示している．

-20
10

20
50
Record duration, K (years)

100

CV [ H R ] = (a y R2 − b y R + c) (γ 50−1) 0.7 / n

図-20 極大値資料から推計された 100 年確率波高を基準(30)酒田港の極大値資料では n = 34.4 × 15.06 = 518 であった．

とする年最大値資料の 100 年確率波高の偏り率

極大値資料に対する上式の定数は，前報の表-8 の算定式を用いて次のように計算される．

が現地データで検証されたといえる．

a = 0.9742×(ln1.29)2 + 0.6949×(ln1.29) + 0.1517 = 0.3918
b = 1.9049×(ln1.29)2 + 1.8142×(ln1.29) + 0.3677 = 0.9532(2) 年最大値から推計された確率波高の信頼区間の試

c = 0.6349×(ln1.29)2 + 1.8561×(ln1.29) + 0.9095 = 1.423

算例
ここで，酒田港の年最大値資料から推計された確率

一方，λ = 15.06 回 / 年のときの再現期間 100 年に対する波高についてその信頼区間の推定を試みる．酒田港の基準化変量は yR = (ln(15.06 × 100))1/1.29 = 4.678 であるので，

年最大値資料に対して L-モーメント法で推計した母数

これらの数値を式 (30) に代入することによって，変動

値と確率波高は次の通りである．

係数が次のように計算される．

18

CV = (0.3918×4.6782 – 0.9532×4.678 + 1.423)

高に比べて過小評価（負の偏り）となり，再現期

× (1.13 – 1.00)0.7/5180.5 = 0.0583

間が長くなるほど過小評価の度合いが強くなる．

8) 上記のことから，高潮偏差など再現期間の長い確すなわち，この事例の場合には，極大値資料から推計率事象を推計する場合には，小標本の年最大値資される 100 年確率波高の変動係数は 5.8% であり，年最料ではなく，極大値資料を使用することが推奨さ

大値資料の場合の変動係数 7.6% よりも 23% 小さい．

れる．

確率波高の変動係数の値を極大値資料と年最大値資料

9) 統計期間40 年の年最大値の標本から推計した 100について試算してみると，平均発生率が低くて再現期

年確率波高は，母集団の値に対して –3% 程度の過間が長い場合には年最大資料のほうが小さいことがあ

小評価となることを例示した．

るけれども，平均発生率が高ければ極大値資料のほう

10) 小標本の年最大値資料の形状母数値が増大し，確が変動係数が小さくなる．この大小関係は条件によっ

率波高が過小評価される傾向は，日本海沿岸 11 地て異なるので，本論文の表-6 および前報の表-8 に記載

点および Kodiak 島沖の高波データで確認された．

した推定式を使って試算してみる必要がある．

11) 年最大値から推計される確率波高の変動係数の経験的推定式が導かれた．

１０．まとめ12) 酒田港の波高データを用いて，年最大値資料によ今回のモンテカルロ法による一連の数値シミュレー

る 100 年確率波高の変動係数の試算ならびに極大ションによって明らかになった諸点をまとめると次の値資料による変動係数との比較の事例を示した．

ようになる．

最初にも述べたように，これまでの高波の極値統計

1) 高波の年間発生回数がポアソン分布に従うことが解析では年最大値資料がほとんど使われていなかった

酒田港およびアラスカ沿岸 Kodiak 島沖の波高デー

けれども，今後は高波データの蓄積とともに年最大値タで例証された．

資料の利用が試みられる可能性がある．しかし本論文

2) モンテカルロ法を用い，一般化極値分布，一般化で述べたように，統計期間が数十年程度の年最大値資

パレート分布，ワイブル分布の 3 種類の母集団か料は確率波高の推計に負の偏りをもたらす危険性が高ら年最大値を抽出し，その極値統計特性を検討しい．既に年最大値資料の利用が一般的である高潮そのた．年最大値の標本は，統計期間が最長 1万年に他の極値統計についても，極大値資料との比較検討を相当するものまでを作成した．

進めることが望まれる．

3) 年最大値の極値データは，統計期間がが数千年以本論文の成果が種々の分野の極値統計解析に活用さ上と十分に長く，平均発生率が年 5 回程度以上でれることを願う次第である．

あれば，母集団が一般化極値分布，一般化パレー
ト分布，ワイブル分布のいずれであっても，一般

付録：非対称分布における標本Skewness の特性

化極値分布に従うことが確認された．

4) 上記の 3 種類の母集団から抽出された大標本の年

本文の 6(2) で述べたように，標本が小さいときには最大値データに当てはめられた一般化極値分布の標本の歪み度 (skewness) が母集団よりも平均的に小さ形状・尺度・位置母数に関して，母集団の母数かくなる．この問題について別途の数値実験を行ったのらの変換式が明示された．

で，以下にその概要を記載する．

5) 実際の統計解析で入手可能な統計期間が数十年程A.1 対象とした非対称分布関数

度の年最大値データ（小標本）の場合には，一般化極値分布への収束は期待できず，3 種類の極値

数値実験では式 (A.1) 〜 (A.3) に示すべき乗関数，一分布関数のいずれもが最適合となる可能性がある． 般化パレート (GPA) 分布関数，およびワイブル分布関6) 上記の減少は，標本の大きさが縮小するにつれて数の 3 種類の関数を対象とした．ただし，取り扱いを

標本の歪み度が低減することに関連している．標

簡単にするため，尺度母数を A = 1，位置母数を B = 0 に本が小さくなるにいつれて L-skewness が減少し，設定した．このうちのべき乗関数は，[0, 1] の範囲での標本に当てはめられる形状母数の値が増大する．

み定義され，確率密度関数 f (x) はk > 1 であれば x = 0 で

7) 小標本の年最大値資料から推計される確率波高は， 最大値を取り，その後次第に減少してx = 1 で 0 となる．
再現期間が 50 年程度以上の場合に母集団の確率波

形状母数が k = 2 のときの f (x) の減少は放物線型である．

19

A. べき乗関数
C. ワイブル分布関数: 0 ≤ x ≤ 1,⎫⎪
k −1
f ( x) = k (1 − x) ,
⎬
⎪
−1
1/ k
x = F ( P) = 1 − (1 − P) .
⎭F ( x) = 1 − (1 − x) k

(A.1)

F ( x) = 1 − exp[− x k ] : 0 ≤ x ≤ ∞,⎫⎪f ( x) = k x k −1 exp[− x k ],
⎬
⎪
−1
1/ k
x = F ( P ) = [− ln(1 − P )] .
⎭B. 一般化パレート (GPA) 分布関数

上記の3 分布関数の平均，標準偏差，歪み度 (Skew-: 0 ≤ x ≤ 1 / k for 0 < k ≤ 1,⎫⎪: 0 ≤ x < ∞ for k < 0,
⎪ (A.2)⎬f ( x) = (1 − kx)1 / k −1
: k ≠ 0,⎪⎪x = F −1 ( P ) = [1 − (1 − P ) k ] / k .
⎭F ( x) = 1 − (1 − kx)1 / k

表-A.1
分布関数

平均
E[x] = µ [x]

べき乗関数

1
k +1

一 般 化 パ
レ ー ト
(GPA)

1
1+ k

ness) などはすべて形状母数 k の関数として表-A.1 のように表される（一部は表-5 と重複）．これに基づいて

3 分布関数の特性値を幾つかの形状母数の値について計算した結果を表-A.2 に示す．

非対称分布関数の特性値の計算式

E[x2]

歪み度 (Skewness)

標準偏差 σ[x]E[x3]2(k + 1)(k + 2)

(A.3)

6
(k + 1)(k + 2)(k + 3)

62
(1 + k )(1 + 2k ) (1 + k )(1 + 2k )(1 + 3k )1 ⎡ k ⎤
k + 1 ⎢⎣ (k + 2) ⎥⎦

β1

1/ 2

2( k − 1)( k + 2)1 / 2
k 1 / 2 ( k + 3)

2(1 − k )(1 + 2k )1/ 2(1 + 3k )1(1 + k )(1 + 2 k )1 / 2
1/ 2

ワイブル

⎛ 1⎞
Γ ⎜1 + ⎟⎝ k⎠⎛ 2⎞
Γ⎜ 1 + ⎟
⎝ k⎠

表-A.2

平均
標準偏差
歪み度

k=2
0.33330.23570.5657

表-A.3
標本サ
イズ n

5
10
20
50
100
1000

A.2

k=3
0.25000.19370.8607

一般化パレート (GPA)k = 0.2
0.83330.7043
1.1832k = -0.2
1.25001.61374.6476ワイブル

k = 1.0
1.0
1.0
2.0

k = 2.0
0.8862
0.46330.6311

標本の大きさによる非対称分布関数の歪み度の簡易計算

べき乗関数

k=2
0.378
0.485
0.533
0.556
0.562
0.566

(E[x3] – 3E[x2]×E[x]+ 2E3[x]) / σ3[x]非対称分布関数の特性値の算定例

べき乗関数

特性値

⎧ ⎛ 2⎞1 ⎞⎫2⎛⎨Γ ⎜ 1 + ⎟ − Γ ⎜ 1 + ⎟ ⎬⎝ k ⎠⎭⎩ ⎝ k⎠⎛ 3⎞
Γ⎜ 1 + ⎟
⎝ k⎠

k=3
0.512
0.685
0.776
0.830
0.847
0.860

一般化パレート (GPA)k = 0.2
0.619
0.858
1.001
1.103
1.142
1.179

k = -0.2
0.919
1.412
1.841
2.346
2.687
3.579

ワイブル

k = 1.0
0.775
1.132
1.394
1.636
1.759
1.949

k = 2.0
0.262
0.392
0.482
0.555
0.587
0.625

とする．これらのデータの非超過確率 Pj は非減少

小標本の歪み度の簡易計算の数値列を形成する．すなわち次の関係を満たす．

ここで標本が小さいときの歪み度を求めること

0 < P1 ≤ P2 ≤ L ≤ Pj ≤ L ≤ Pn < 1を考える．大きさが n の標本中のデータ x が昇順，

(A.4)

非超過確率 Pj の値は標本によって異なるけれど

すなわち x1 ≤ x2 ≤ L x j L ≤ xn と並べ替えられているも，多数の標本についての平均，すなわち期待値

20そこで，E[Pj] がその中間値である (j – 1/2) / n の値を-1

Probability density function of skewness,

とると仮定すると，各データの期待値 E[xj] は F [(j– 1/2) / n] として計算できる．そうすると，標本の平均 µ [x]，分散 σ [x]，および歪み度 β1 はデータ x
2の 1 〜3 次の積率を用いて次のように求められる．

1 nµ[ x] = ∑ x j =E[ x]n j =1(A.5)

1 n∑ ( x j − µ[ x]) 2 =E[ x 2 ] − E 2 [ x]n j =1(A.6)(

n = 10n = 100
n = 1,000
n = 10,000

6

4

2

0
-1.0

)

(A.7)

0

0.5

1.0

1.5

2.0

2.5

β1

1 n 2
1 n
E [ x j ], E[ x 3 ] = ∑ E 3 [ x j ]∑n j =1n j =1

0.8

GPA Distribution (k = -0.2)

Probability density function of skewness,

ここに，

(A.8)

標本の大きさを n = 5 〜1000 で変化させて式 (A.1)
〜 (A.3) の 3 種類の分布関数についてついて歪み度を計算した結果を表-A.3 に示す．標本が小さくなる

0.6

n = 10n = 100
n = 1000
n = 10,000

0.4

0.2

0
-2

につれて歪み度が減少することが明らかに見てとれ
る．n = 1000 の場合，べき乗関数では母集団の歪み

02

4Skewness,

6

8

10

β1

図-A.2 GPA分布 (k = –0.2) の歪み度の頻度分布

い一般化パレート分布の k = –0.2 の場合には， n =1000 であっても母集団よりも小さな値にとどまる．
なお平均 µ [x] および σ [x] についても E[xj] = F [(jProbability density function of skewness,-1

f (β1)

度とほぼ同じ値が得られるが，非対称の度合いが強– 1/2) / n] と仮定して計算を行ったが，標本サイズの影響は僅かであった．
標本が小さいときに歪み度が母集団よりも小さくなるのは，データの最大値 xn の非超過確率の期待値
がPn = 1– 1/2nであるために最大値の期待値も低く抑えられ，その結果として式 (A.8) の 3 次の積率

E[x3] が小さな値になるためである．

5

Exponential distribution(Weibull with k = 1.0)
4
n = 10n = 100
n = 1000
n = 10,000

3

2

1

0
-1

A.3

-0.5

図-A.1 べき乗関数 (k = 3) の歪み度の頻度分布

1
E[ x 3 ] − 3 E[ x 2 ] × E[ x ] + 2 E 3 [ x ]3
σ [ x]

E[ x 2 ] =Algebraic Function
with k = 3.0

Skewness,

1 1 n( x j − µ[ x]) 3
β1 = 3∑σ [ x] n j =1
=

10

8

f (β1)σ 2 [ x] =f (β1)

E[Pj] は (j – 1) / n と j / n の間の値を取ると考えられる．

数値実験による非対称分布の歪み度の算定

01
Skewness,

2

3

4

β1

図-A.3 ワイブル分布 (k = 1.0) の歪み度の頻度分布

次に，モンテカルロ法を用いて 3 種類の分布関数
から多数の標本を抽出し，標本毎の歪み度を計算した結果を述べる．標本の大きさは n = 10 〜10,000 で

数値実験で得られた歪み度の頻度分布を確率密度

段階的に変化させ，標本数は n = 1000 以下では N =関数の形で示したのが図-A.1 〜A.3 である．それぞ

100,000 とし，n = 2,000 〜10,000 では計算時間を節約

れべき乗関数 (k = 3)，一般化パレート分布 (k = –0.2)，

するために N = 10,000,000 / n とした．

およびワイブル分布 (k = 1.0：指数分布) の場合であ

21

る．

度の分布幅が広がる．

頻度分布は級幅 ∆β1 = 0.1 で作成している．なお，

図-A.4 〜A.6 は数値実験で得られた歪み度の平均標本が n = 10 と小さい場合には，級毎の頻度がや

値と標本サイズとの関係を，表-A.3 に示した簡易計や不規則に変動したため，連続する三つの級にわ算の結果と比較して示したものである．簡易計算はたる移動平均を行って平滑化している． 3 分布関

n = 1000 までしか行っていないので，n = 10,000 につ

数いずれのケースも，標本が小さくなるにつれていては表-A.2 の母集団の値を記載している．数値実

頻度分布が左にシフト（歪み度が減少）し，歪み

験の結果は Sim，簡易計算は Cal の凡例で示してある．これらの図に見られるように，両者の値には若1.0

干の差があるものの，標本の大きさに応じて歪み度

Algebraic Functionが変化する様相は同じである．

Sample skewness,

β1

0.9

この結果は，A.2 節で紹介した考え方，すなわち0.8

小標本における歪み度の低下は標本中の最大値が低0.7

めに抑えられることに起因するというメカニズムを0.6

確認するものである．本文の図-11, 12 に示した年最0.5

大値資料の歪み度が統計期間の長さによって変化すk = 2 (Cal)k = 2 (Sim)k = 3 (Cal)k = 3 (Sim)

0.4

る現象は，こうした非対称分布の特性によるもので
ある．

0.3
1

10

100

1000

10000

Sample size, n

参考文献

図-A.4 標本サイズと乗関数の歪み度の関係

合田良実 (2008)：「耐波工学 – 港湾・海岸構造物の耐波設計–」，鹿島出版会，430p.

5.0

GPA Distribution

4.0

合田良実 (2010)：L-モーメント法を用いた高波の極値統計解析におけるプロッティングポジション公式の選択と各種統計量の信頼区間，ECOH/YG 技術論文 No. 13,
2010年9月30日，29p.

Sample skewness,

β1

3.0

2.0

合田良実・久高将信 (2010)：L-moment法による日本海沿岸の高波の地域頻度解析について，ECOH/YG 技術論
文 No. 12, 2009年4月12日，13p.
1.0

Chowdhury, J. U., J. R. Stedinger, and L.-H. Lu (1991):
Goodness-of-fit tests for regional generalized extreme
value flood distributions, Water Resources Res., Vol.
27, No. 7, pp. 1765-1776.

k =-0.2 (Cal)k =-0.2 (Sim)k = 0.2 (Cal)k = 0.2 (Sim)

0.8
0.6
1

10

100

1000

Cunnane, C. (1973): A particular comparison of annual
maxima and partial duration series methods of flood
frequency prediction, J. Hydrology, Vol. 18, pp. 257-271.

10000

Sample size, n

図-A.5 標本サイズとGPA分布の歪み度の関係

Hosking, J. R. M. (1990): L-moments: Analysis andestimation of distributions using linear combinations oforder statistics, J.Roy. Statistical Soc., Series B, 52, pp.
105-24.

3.0

Weibull Distribution
2.0

Sample skewness,β1

Hosking, J. R. M. and J. R. Wallis (1997): “Regional
Frequency Analysis,” Cambridge Univ. Press, 224p.
1.0

Madsen, H., C. P. Pearson, and D. Rosjerg (1997):
Comparison of annual maximum series and partial
duration series methods for modeling extreme hydrogic
events: 2. Regional modeling, Water Resources Res.,
Vol. 33, No. 4, pp. 759-769.

0.8
0.6
k = 1.0 (Cal)k = 1.0 (Sim)k = 2.0 (Cal)k = 2.0 (Sim)

0.4

0.2
1

10

100

1000

10000

Sample size, n

図-A.6 標本サイズとワイブル分布の歪み度の関係

22