数理解析研究所講究録 1308 巻 2003 年 61-74

61 61
上位

$r$

個の観測値に基づく確率点の推定
神戸商船大学高橋倫也 (Rinya Takahashi)

Kobe University of Mercantile Marine
高千穂大学 渋谷政昭 (Masa 硅 Sibuya)
Takachiho University$\mathrm{i}$

いくつかの単位期間または単位領域のそれぞれで測定した確率標本データの, 各標本上位 個を用いる極値解析法について述べる. 上側微小確率点推定の精度とその問題点について議論する.
$r$

1

はじめに
確率標本より, 与えられた期間または領域での最大値を予測するために, いくつかの単位期間

または単位領域ごとの最大値 (極値) データを利用するのが古典的極値解析の手法である. その予測精度を上げる方法として 「上位 個までのデータ」 または 「ある閾値以上の全てのデータ」
を用いることが提案されている. ここでは前者の場合について,「推定精度がどの程度改善するか」$r$と「この方法の問題点」 について議論する.
Weissman (1978) により, 極値理論に基づく上位 個のデータ (一組) を用いる上側微小確率点の推定が初めて議論された. 彼は, 母集団分布が Gumbel 分布の吸引領域に属す場合 (Gumbelモデル) を詳しく調べた. Smith (1986) は, Gumbel モデルの下で Weissman (1978) を発展さの決定法を提案し, 情
せ, 上位 個のデータを $N$ 組用いる場合の推定問題を議論した. 彼は,
$k$

$r$

$r$

報量を計算し, ベニスの毎年上位 10 個の潮位データを位置パラメータに時間依存性を導入して解析し, この手法の有効性を示した. Tawn (1988) は, Smith (1986) の結果を一般極値 (GEV)
モデル
モデルヘ拡張し, Lowestoft の潮位データの解析を行った. Scarf et al. (1992) は,$\mathrm{G}\mathrm{E}\mathrm{V}$の下で位置と尺度パラメータに時間依存性を導入し, 上位$r$個の孔食深さデータの解析を行い手法の有効性を示した.
以下, 2 節で極値理論から導かれる上位 個の順序統計量の漸近同時分布を示しその性質につの決定について述べ, 4 節 上位 個を用いる有
いてまとめる. 3 節でパラメータの推定法と効性について議論する. 5 節で実データの解析例を示す. 付録で上位 個の漸近同時分布の情報$r$$k^{\backslash }\backslash$

$r$

$r$

$r$

量を述べる.

2

極値理論
一般極値 (GEV) 分布の標準型を$G_{\xi}(z)=\exp[-(1+\xi z)^{-1/\xi}]$とする. ここで$\exp(-e^{-z})$$G_{\xi}$は,$\xi<0$ の場合はで Gumbel 分布,$\xi>0$,$1+\xi z>0(\xi\in \mathrm{R})$

,

(2.1)

Weibull 分布, $\xi=0$ の場合は Go(z)=lim\mbox{\boldmath $\xi$}\rightarrow o $G\xi(z)=$の場合は Fr\’echet 分布である.

62とし, 分布
の順序統計量を分布 $F$ からの確率標本$F$ が GEV 分布の吸引領域に属す $(F\in D(G_{\xi}))$ と仮定する : すなわち適当な数列 $a_{n}>0$ ,$\mathrm{Y}_{1},$$\mathrm{Y}_{1:n}\geq \mathrm{Y}_{2:n}\geq\cdots\geq \mathrm{Y}_{n:n}$$\mathrm{Y}_{n}$$\mathrm{Y}_{2},$$\ldots,$$G_{\xi}$$b_{n}\in \mathrm{R},$$n=1,2,$$\ldots$が存在し,$\lim_{narrow\infty}P(.\frac{\mathrm{Y}_{1.n}-b_{n}}{a_{n}}\leq z)=G_{\xi}(z)$このとき, 上位$r$,$\forall z\in \mathrm{R}$.

個の順序統計量の同時分布関数$P(. \frac{\mathrm{Y}_{1.n}-b_{n}}{a_{n}}\leq z_{1},\cdot\frac{\mathrm{Y}_{2.n}-b_{n}}{a_{n}}\leq z_{2},$$\cdots,\cdot\frac{\mathrm{Y}_{r\cdot n}-b_{n}}{a_{n}}\leq z_{r})$,$z_{1}\geq z_{2}\geq\cdots\geq z_{r}$は同時密度関数$g_{\xi,12\cdots r}(z_{1}, z_{2}, \ldots, z_{r})=\frac{g_{\xi}(z_{1})\cdots g_{\xi}(z_{r-1})}{G_{\xi}(z_{1})\cdots G_{\xi}(z_{r-1})}g_{\xi}(z_{r})$を持つ分布関数
分布
布関数,$G_{\xi,12\cdots r}$,$g_{\xi}(z)=dG_{\xi}(z)/dz$(2.2)に収束する.

に従う確率ベクトルを $(Z_{1}, Z_{2}, \ldots, Z_{r})$ とする. このとき,
によらず,
, は

$G_{\xi,12\cdots r}$$G_{\xi,j}$,$Z_{j},$$j\geq 1$

の周辺分

$r$

$\sum_{k=0}^{j-1}(1+\xi z)^{-k/\xi}\exp\{-(1+\xi z)^{-1/\xi}\}/k!$,$\xi\neq 0$

,
(2.3)

$G_{\xi,j}(z)=\{$$j-1 \sum\exp(-kz-e^{-z})/k!$,$\xi=0$,$k=0$

で, その周辺密度関数,$g_{\xi,j}$,は$(1+\xi z)^{-j/\xi-1}\exp\{-(1+\xi z)^{-1/\xi}\}/\Gamma(j)$,$\xi\neq 0$,$\xi=0$,$g_{\xi,j}(z)=\{$$\exp(-jz-e^{-z})/\Gamma(j)$

,

(2.4)

となる.

定理 1. $(Z_{1}, Z_{2}, \ldots, Z_{r})$ は次の性質をもつ. ただし従う独立確率変数, $S_{j}= \sum_{k=1}^{j}E_{k}$ とする.

(A) GEV $(\xi\neq 0)$ モデルの場合 :$=(\Gamma(j-\xi)-\Gamma(j))/(\xi\Gamma(j)),$1)$E_{\xi}(Z_{j})$2)3)$1_{\xi}(Z_{j})$$E_{1},$$E_{2},$$\ldots$を標準指数分布$=(\Gamma(j-2\xi)\Gamma(j)-\Gamma^{2}(j-\xi))/(\xi^{2}\Gamma^{2}(j))$.

$\mathrm{C}\mathrm{o}\mathrm{r}_{\xi}(Z_{j}, Z_{j+1})=\frac{1}{j-\xi}\sqrt\frac{\Gamma(j+1-2\xi)\Gamma(j+1)-\Gamma^{2}(j+1-\xi)}{\Gamma(j-2\xi)\Gamma(j)-\Gamma^{2}(j-\xi)}$$\{Zj, j\geq 1\}=d\{(S_{j}^{-\xi}-1)/\xi, j\geq 1\}$(B) Gumbel $(\xi=0)$ モデルの場合 :1) $E_{0}(Z_{j})=-\psi(j),$ $V_{0}(Z_{j})=\psi’(j)$ .
2).
.

$\mathrm{C}\mathrm{o}\mathrm{r}_{0}(Z_{j}, Z_{j+1})=\sqrt{\psi’(j+1)\psi’(j)}$3)$\{Z_{j}, j\geq 1\}=\{-\log S_{j}d, j\geq 1\}$

ただし,$\Gamma$はガンマ関数,

注 1. (連続性) (A) で$\psi$.

はディ. ガンマ関数,$\xiarrow 0$$\psi’$

はトリ・ガンマ関数である.

とすれば (B) が得られる.

に$(\mathrm{E}\mathrm{x}\mathrm{p}(1))$.

63 63
注 2.

(A) 3), (B) 3) は Nagaraja (1982), (B) 1) は Weissman (1978) 参照.

. ガンマ関数の値は,

注 3. 変数が正整数の場合, デイ. ガンマ関数とトリ$\psi(n)=-\gamma+\sum_{i=1}^{n-1}\frac{1}{i}$から求まる. ただし,,$\gamma=0.57721566\ldots$$\psi’(n)=\frac{\pi^{2}}{6}-\sum_{i=1}^{n-1}\frac{1}{i^{2}}$,$n=1,2,$$\ldots$

は Euler の定数である.

注 4.
0.626,

Coro $(Zj, Zj+1)\#\mathrm{h}j$ に関して狭義の増加関数で 1 に収束する. 例えば, Coro

注 5.

(B) 3) より, Gumbel モデルの下では$\mathrm{C}\mathrm{o}\mathrm{r}_{0}(Z_{2}, Z_{3})=0.783$,$\mathrm{C}\mathrm{o}\mathrm{r}_{0}(Z_{3}, Z_{4})=0.848$$j(Z_{j}-Z_{j+1})=j \log\frac{S_{j+1}}{S_{j}}$は互いに独立にに従う (Weissman, 1978)$\mathrm{E}\mathrm{x}\mathrm{p}(1)$$\frac{j}{\xi}\log\frac{1+\xi Z_{j}}{1+\xi Z_{j+1}}$

が互いに独立に
注 6.

,

である.

,

$j=1,2,3,$

GEV $(\xi\neq 0)$ モデルの場合は

従って,$j=1,2,3,$

$\xi$

(2.6)

$\ldots$の標準一般 Pareto 分布,$1+\xi y>0$ ,$P(y)=1-(1+\xi y)^{-1/\xi}$ ,( $\xi\geq 0$ のとき $y>0,$(2.5)$\ldots$に従う (Tawn, 1988).

$\mathrm{E}\mathrm{x}\mathrm{p}(1)$形状パラメータ$\cdot$$(Z_{1}, Z_{2})=$$\xi<0$のとき $0<y<-1/\xi$ ) からの $n$ 個の順序統計量を$\mathrm{Y}_{1:n}\geq \mathrm{Y}_{2:n}\geq\cdots\geq \mathrm{Y}_{n:n}$とすると,$\mathrm{Y}_{j:n}=\frac{U_{n-j+1-n}^{-\xi}-1}{\xi}d$と表される. ただし,計量である. この,$j=1,2,$

$n$

,

は一様分布

$1\geq U_{1:n}\geq U_{2:n}\geq\cdots\geq U_{n:n}\geq 0$$\mathrm{Y}_{j:n}$$\ldots,$$\mathrm{U}(0,1)$

からの $n$ 個の順序統と (A) 3) のろ $=d(S_{j}^{-\xi}-1)/\xi$ から $(Z_{1}, Z_{2}, \ldots, Z_{r})$ は形状パラメー$p$個の順序統計量と見なせる (Tawn, 1988). 一方, Gumbel モデ
ルの下では, $(Z_{1}, Z_{2}, \ldots, Z_{r})$ は指数分布からの上位 個の順序統計量と見なせる (Weissman,$\xi$の一般 Pareto 分布からの上位

$r$

$r$

1978).
の周辺密度関数である.
図 1 は, それぞれ $\xi=-0.4,0,0.4$ の場合の
図 2, 3, 4 はそれぞれ, $\xi=-0.4,0,0.4$ の場合の $(Z_{1}, Z_{2})$ の同時密度関数$Z_{1},$$Z_{2},$$Z_{3}$$(1+\xi z_{1})^{-1/\xi-1}(1+\xi z_{2})^{-1/\xi-1}\exp\{-(1+\xi z_{2})^{-1/\xi}\}$,$\xi\neq 0$,$\xi=0$,$g_{\xi,12}(z_{1}, z_{2})=\{$$\exp(-z_{1}-z_{2}-e^{-z_{2}})$$z_{1}\geq z_{2}$,とその等高線である.

,

64
Weibull Model

-3

-4

-2

0

-’

Frechet Mode

GumbeI ModeI

’

-t

2

0

4

$[$

......

$-\mathfrak{l}\overline{-}$’c2$-|.\overline{-}\mathrm{a}$

-2

-’

0

図 1. $\xi=-0.4,0.0.4$ の場合の $Zj$ の周辺密度関数. $j=1,2,3$

禍 eibull ModeI

禍.ibuIl$\mathrm{M}\mathrm{o}\mathrm{d}\cdot \mathrm{I}$$\tau$

図 2. $\xi=-0.4$ の場合の $(Z_{1}, Z_{2})$ の同時密度関数とその等高線.

.

65 65
Gumbel Modei

Gumbel Model$\approx$

図 3.

$\xi=0$の場合の $(Z_{1}, Z_{2})$ の同時密度関数とその等高線.

Frechet Model

Frechet Model$\mathrm{w}$

図 4.

3

$\xi=0.4$の場合の $(Z_{1}, Z_{2})$ の同時密度関数とその等高線.

パラメータ推定
パラメータ $\theta=(\mu, \sigma)$ または$(\mu, \sigma, \xi)$と,極値理論で重要な上側微少確率点 $T$-return level$q(T)$ ,$G_{\xi}( \frac{q(T)-\mu}{\sigma})=1-\frac{1}{T}$

,

すなわち,
$\mu+\sigma\{(-\log(1-1/T))^{-\xi}-1\}/\xi$$q(T)=\{$$\mu+\sigma\{-\log(-\log(1-1/T))\}$ ,,$\xi\neq 0$,$\xi=0$,(3.1)の推定法について述べる. これらのパラメータが補助変量の関数である場合に拡張することもで

66
3.1

モデル

上位

$r$

個の確率ベクトル $(X_{1}, X_{2}, \ldots, X_{r})$ の従う結合密度関数は Gumbel モデルの場合は,$\theta=(\mu, \sigma)$

として,$g(x_{1}, x_{2}, \ldots, x_{r}; \theta)=\frac{1}{\sigma^{r}}\exp[-\sum_{j=1}^{r}(\frac{x_{j}-\mu}{\sigma})-\exp(-\frac{x_{r}-\mu}{\sigma})]$$x_{1}\geq x_{2}\geq\cdots\geq x_{r}$でGEV モデルの場合は$\theta=(\mu, \sigma, \xi)$

,

(3.2)

,

として,

.
$g(x_{1}, \ldots, x_{r};\theta)=$ $\frac{1}{\sigma^{r}}\exp\{-(\frac{1}{\xi}+1)\sum_{j=1}^{r}\log[1+\xi(\frac{x_{j}-\mu}{\sigma})]-[1+\xi(\frac{x_{r}-\mu}{\sigma})]^{-1/\xi}\}$

,

(3.3)
$x_{1}\geq\cdots\geq x_{r}$,$j=1,$ $\ldots,$$1+\xi(x_{j}-\mu)/\sigma>0,$

$r$

,

となる.

データは $x_{1}\geq x_{2}\geq\cdots\geq x_{r}$ が$n$組, すなわち,,,$x_{11}\geq x_{12}\geq\cdots\geq x_{1r}$$x_{21}\geq x_{22}\geq\cdots\geq x_{2r}$$\{$$x_{i1}\geq x_{i2}\geq\cdots\geq x_{ir}$,$x_{n1}\geq x_{n2}\geq\cdots\geq x_{nr}$,の $r\cross n$ 個の数値とする.

32Gumbel モデル
Gumbel モデルの場合の対数尤度は$l( \mu, \sigma)=-nr\log\sigma-\sum_{i=1}^{n}[\sum_{j=1}^{r}(\frac{x_{ij}-\mu}{\sigma})+\exp(-\frac{x_{ir}-\mu}{\sigma})]$

,

(3.4)

となる. パラメータ $\theta=(\mu, \sigma)$ を最尤法で推定する. 尤度方程式は, 次の様になる :$\frac{\partial}{\partial\mu}l(\mu, \sigma)$

$=$

-

,

$\sum_{i=1}^{n}[-\frac{r}{\sigma}+\frac{1}{\sigma}\exp(-\frac{x_{ir}-\mu}{\sigma})]=0$

(3.5)

$\{$

$\frac{\partial}{\partial\sigma}l(\mu, \sigma)$$=$.

$- \frac{nr}{\sigma}+\sum_{i=1}^{n}[\sum_{j=1}^{r}(\frac{x_{ij}-\mu}{\sigma^{2}})-\frac{x_{ir}-\mu}{\sigma^{2}}\exp(-\frac{x_{ir}-\mu}{\sigma})]=0$

この連立非線形方程式から,$h_{r}( \sigma):=\sum_{i=1}^{n}(\frac{x_{ir}-x_{r}=}{\sigma}+1)\exp(-\frac{x_{ir}-x_{r}=}{\sigma})=0$,$\overline{\overline{x}}_{r}=\frac{1}{nr}\sum_{i=1j}^{n}\sum_{=1}^{r}x_{ij}$(3.6)

67 67が得られる. そこで,
$h_{r}’( \sigma)=\frac{1}{\sigma}\sum_{i=1}^{n}(\frac{x_{ir}-x_{r}=}{\sigma})^{2}\exp(-\frac{x_{ir}-=x_{r}}{\sigma})$を用いて, ニュートン法で

を求める.

$r$

これから,

(3.7)

$\hat{\mu}_{r}=-\hat{\sigma}_{r}\log[\frac{1}{nr}\sum_{i=1}^{n}\exp(-\frac{x_{ir}}{\hat{\sigma}_{r}})]$

が求まる. このが上位
る. $T$-return level $q(T)$ の推定は,$\hat{\theta}_{r}=(\hat{\mu}_{r},\hat{\sigma}_{r})$$r$個のデータを用いた場合の$\theta=(\mu, \sigma)$,$\hat{q}_{r}(T)=\hat{\mu}_{r}+\hat{\sigma}_{r}\{-\log(-\log(1-1/T))\}$とすればよい. また, 推定値$\hat{\mu}_{r},\hat{\sigma}_{r},\hat{q}_{r}(T)$の最尤推定値であ(3.8)の標準誤差は付録 Al の漸近分散行列を用いて推定する.

の決定 : ここでは上位 個のデータ $\{(x_{i1}, x_{i2}, \ldots, x_{ir}), i=1,2, \ldots, n\}$ が分布うと見なせる最大の の決定法について議論する.
$r$

$r$

に従

$G_{0,12\ldots r}$

$r$

上位

$j$

番目の確率変数$X_{j}$

は次の周辺分布関数$G_{0,j}(z)=P \{\frac{X_{j}-\mu}{\sigma}\leq z\}=\sum_{k=0}^{j-1}\exp(-kz-e^{-z})/k!$を持つ. 従って,(3.9)$U_{\dot{\iota}j}=G_{0,j}( \frac{X_{ij}-\mu}{\sigma})$

により,からの確率変数を一様分布の決定法が考えられる.
このことから, 次の$\mathrm{U}(0,1)$$X_{ij}$$U_{ij}$に変換できる.

$r$

PP plot :を固定し, 上位$r$$r$個のデータから推定値$\hat{\mu}_{r},\hat{\sigma}_{r}$を求める. これらを用いて, 上位番目のデータ $\{x_{1j}, x_{2j}, \ldots, x_{nj}\}(j=1,2, \ldots, r)$ からの順序統計量を $u_{(1)j}\geq u_{(2)j}\geq\cdots\geq u(n)j$ とし,を求める.

$u_{ij}=G_{0,j}((x_{ij}-\hat{\mu}_{r})/\hat{\sigma}_{r}),$$i=1,2,$$\ldots,$

$j$

$n$

$u_{ij}$$(1- \frac{i}{n+1},$をプロットし 個のると見なせる最大の$r$$\mathrm{P}\mathrm{P}$$r$

plot を作成する.

を決定する.

$u_{(i)j})$,$i=1,2,$ $\ldots,n$ ,

ここで,

$r$

個すべての

$\mathrm{P}\mathrm{P}$plot が直線性を示して $1$(Smith, 1986 参照. )一方, 次の決定法も考えられる.

の場を固定する. $j=1$ の場合は通常の Gumbel 確率紙を考える. $j=2,3,$合, 上位 番目のデータ $\{x_{1j}, x_{2j}, \ldots, x_{nj}\}$ の順序統計量を $x_{(1)j}\geq x(2)j\geq\cdots\geq x(n)j$ とする.
数値計算で確率点 $q_{(i)j}=G_{0,j}^{-1}(1-i/(n+1))$ を求め,

QQ plot :$\ldots,$

$r$

$r$

$j$

$(x_{(i)j}, q_{(i)j})$をプロットし$r$

なせる最大の$r$,$i=1,2,$plot を作成する. すべてのを決定する.

個の$\mathrm{Q}\mathrm{Q}$$r$$\ldots$

個の, n》$\mathrm{Q}\mathrm{Q}$plot が直線性を示して 4) ると見68
plot による方法では周辺分布の適合しか見てぃない. 分布の同時性をチェックする方法として, 2 節の注 5 より次のものが考えられる.
これら$\mathrm{P}\mathrm{P},$$\mathrm{Q}\mathrm{Q}$

指数確率紙 : $j=1,2,$

に対して,

{

$n$
$xij+1,$ $i=1,2,$} を指数確率紙にプロットする.
すべての $j(<r)$ の指数確率紙で直線性を示してぃると見なせる最大の を決定する.
$\ldots$$x_{ij}$一$\ldots,$$r$

従って, PP または QQ と指数確率紙をプロットし得られた
よい.

$r$

の中で最小のものを採用すれば33GEV モデル
GEV モデルの場合の対数尤度は$l( \mu, \sigma, \xi)=-nr\log\sigma-.\cdot\sum_{=1}^{n}\{(\frac{1}{\xi}+1)\sum_{j=1}^{r}\log[1+\xi(\dot{.}\frac{x_{j}-\mu}{\sigma})]+[1+\xi(.\cdot\frac{x_{r}-\mu}{\sigma})]^{-1/\xi}\}$(3.10)となる.
パラメータ $\theta=(\mu, \sigma, \xi)$ を最尤法で推定する. 尤度方程式は簡単にはならない. ニュートン法

で連立非線形の尤度方程式を解かなければならないが, 初期値としては極値データから PWM 法(Hosking et al., 1985) で求めた推定値を用いればよい. 得られた最尤推定値をとすると -return level $q(T)$ の推定は,$\hat{\theta}_{r}=(\hat{\mu}_{r},\hat{\sigma}_{r},\hat{\xi}_{r})$$T$(3.11)$\hat{q}_{r}(T)=\hat{\mu}_{r}+\hat{\sigma}_{r}\{(-\log(1-1/T))^{-\hat{\xi_{r}}}-1\}/\hat{\xi}_{r}$とすればよい. また, 推定値
いて推定する.
$r$の決定 : ここでも, 上位従うと見なせる最大のGEV モデルの下で,

$r$

$r$

$\hat{\mu}_{r},\hat{\sigma}_{r},\hat{\xi}_{r\prime}\hat{q}_{r}(T)$の標準誤差は付録 A2 の Fisher 情報行列を用

個のデータ $\{(x_{i1}, x_{i2}, \ldots, x_{ir}), i=1,2, \ldots,n\}$ が分布$G_{\xi,12\ldots r}$にの決定法について議論する.

$X_{j}$

は次の周辺分布関数$G_{\xi,j}(z)=P \{\frac{X_{j}-\mu}{\sigma}\leq z\}=\sum_{k=0}^{j-1}(1+\xi z)^{-k/\xi}\exp\{-(1+\xi z)^{-1/\xi}\}/k!$を持つ. 従って,(3.12)$U_{ij}=G_{\xi,j}( \frac{X_{ij}-\mu}{\sigma})$により,$X_{\dot{l}j}$を一様分布$\mathrm{U}(0,1)$

からの確率変数$U_{ij}$に変換できる.

このことから, 次の決定法が考えられる.

plot:を固定し, 上位 個のデータから推定値
番目のデータ $\{x_{1j}, x_{2j}, \ldots, x_{nj}\}(j=1,2, \ldots, r)$ から$\mathrm{P}\mathrm{P}$

を求める.

$r$

$r$

$u_{ij}$$\hat{\mu}_{r},\hat{\sigma}_{r},\hat{\xi_{r}}$を求める. これらを用いて, 上位$u_{ij}=G_{\hat{\xi_{r}},j}((x_{ij}-\hat{\mu}_{r})/\hat{\sigma}_{r}),$$i=1,2,$$\ldots,$

$j$

$n$

の順序統計量を $u_{(1)j}\geq u_{(2)j}\geq\cdots\geq u_{(n)j}$ とし,$(1- \frac{i}{n+1},$$u_{(:)j})$,$i=1,2,$$\ldots,$をプロットし 個のplot を作成する. ここで, 個すべてのと見なせる最大の を決定する. (Tawn, 1988 参照. )
$r$

$n$

,

$\mathrm{P}\mathrm{P}$

$r$

$r$

$\mathrm{P}\mathrm{P}$plot が直線性を示してぃる

69 69また, 次の方法も考えられる.

QQ plot :を求める. 上位 番個のデータから形状パラメータの推定値
目のデータ $\{x_{1j}, x_{2j}, \ldots, x_{nj}\}(j=1,2, \ldots, r)$ の順序統計量を $x_{(1)j}\geq x_{(2)j}\geq\cdots\geq x_{(n)j}$ とす
る. 数値計算で確率点 $q(\text{由}=G_{\hat{\xi}_{r},j}^{-1}(1-i/(n+1))$ を求め,
$r$を固定し, 上位$\hat{\xi}_{r}$$r$$(x_{(i)j}, q_{(i)j})$をプロットし 個のるとみなせる最大の$r$$\mathrm{Q}\mathrm{Q}$$r$plot を作成する.

$i=1,2,$

,

ここで,

$r$

$\ldots$

$j$

|n》

個すべての

$\mathrm{Q}\mathrm{Q}$plot が直線性を示していを決定する.

分布の同時性をチェックする方法として,

2 節の注 5 より次のものが考えられる.

指数確率紙 :r. を固定し, 上位 個のデータから推定値,
データ $.\{x_{1j}, x_{2j}, \ldots, x_{nj}\}(j=1,2, \ldots, r-1)$ [こ対して.
$r$$\frac{j}{\hat{\xi}_{r}}\log\frac{\hat{\sigma}_{r}+\hat{\xi}_{r}(x_{1j}-\hat{\mu}_{r})}{\hat{\sigma}_{r}+\hat{\xi_{r}}(x_{ij+1}-\hat{\mu}_{r})}.$,$\hat{\mu}_{r},\hat{\sigma}_{r},\hat{\xi}_{r}$$i=1,2,$$\ldots,$$n$を求める. 上位

$j$

番目の

,

を指数確率紙にプロットし $r-1$ 個の指数確率紙を作成する. ここで $r-1$ 個すべての指数確率
紙で直線性を示していると見なせる最大の を決定する.
$r$

従って,

PP または QQ と指数確率紙をプロットし得られた$r$の中の最小のものを採用すれば

よい.

4

有効性
ここでは,

般の GE

.

Gumbel モデルの下で上位 $r(\geq 2)$ 個のデータを用いる有効性について議論する.

.
好皀妊襪硫爾任竜掴世脇颪靴い, 以下の議論は $\xi=$$0$–の場合にも近似的に成り立つと考えられる.

個のデータを用いる場合の $T$-return level $q(T)$ の推定量の有効性について調べる.
の漸個を用いる有効性を示す. 推定量$q(T)$ の推定を行う場合の漸近分散の比で上位
近分散は, 付録 Al の (A1.4) より上位$r$$\hat{q}_{r}(T)$$r$$AV( \hat{q}_{r}(T))=\frac{\sigma^{2}}{n(rC_{r}-B_{r}^{2})}[C_{r}+(g(T))^{2}r+2g(T)B_{r}]$となる. 従って,$\hat{q}_{1}(T)$に対する$\hat{q}_{r}(T)$の漸近相対効率は$e( \hat{q}_{r}(T),\hat{q}_{1}(T))=\frac{AV(\hat{q}_{1}(T))}{AV(\hat{q}_{r}(T))}=\frac{(rC_{r}-B_{r}^{2})[C_{1}+(g(T))^{2}+2g(T)B_{1}]}{(C_{1}-B_{1}^{2})[C_{r}+(g(T))^{2}r+2g(T)B_{r}]}$と表される. これは,

$T$

と

$r$

の関数である,$T=1\mathrm{O}\mathrm{O},$に対する $e(\hat{q}_{r}(T),\hat{q}_{1}(T))$ の値は次の表の様になる :$1,000,10,000,100,000$ と $r=2,3,$(4.1)$\ldots,$

$10$

70
この表より, 例えば 10,000-return level を推定するとき,100 個の極値データは 50 組の上位 3個のデータとほぼ同じ精度, ほぼ等しい漸近分散, を持つと言うことが出来る.
ここでの議論は, 「上位$r$個のデータが正確に想定した漸近同時分布からのものである」 と仮定している事に注意する必要がある.

5

実データ解析
ここでは 「上位 3 個までの孔食深さ測定」 データの解析を行う. データは面積が等しい 30 個の領域内において上位 3 個の孔食深さを測定したものである.
図 5 は, データの (1 位, 2 位), (2 位, 3 位) の散布図である. 相関係数はそれぞれ 0813 と0851 である.
以下,Gumbel と GEV モデルの下で解析した結果を紹介する.
(First. Second)(Second, Third)$\mathrm{s}$$\ovalbox{\tt\small REJECT}$

Dl

D2

図 5. (a) データの (上位 1 位, 上位 2 位). (b) データの (上位 2 位, 上位 3 位).

5.1

Gumbel モデル

先ず, 極値データ $(r=1)$ が Gumbel 分布に従うと見なしてよいかを調べた. すなわち,分布で形状パラメータ$\xi$について, 検定

(Hosking et al., 1985) で行った.
$p$

$\xi$

の

$\mathrm{H}_{0}$: $\xi=0,$$\mathrm{P}\mathrm{W}\mathrm{M}$$\mathrm{H}_{1}$

推定値は

:

$\hat{\xi}=$$\xi\neq 0$

GEV

を PWM 推定値に基づく方法-0.191, 検定統計量の値は 1395 で値は 0.163 であった.

上位 $r=1,2,3$ 個のデータを用いた場合の$r=1$ ,$\hat{\mu}_{r}$$\hat{\sigma}_{r}$$r$の決定を考える.

$\mathrm{Q}\mathrm{Q}$

$\mu$

と

$\sigma$2,0.420 0.423
0.037 0.036の最尤推定値は次の様になった :

3
0.431
0.045

plot を書かせたのが図 6 で, 指数確率紙が図 7 である. これらの図よ

このデータでは上位 2 個まで使えると考えられる.
従って, $T$-return level の推定は次の様になる :り,$\ovalbox{\tt\small REJECT}(T)=0.423+0.036\{-\log(-\log(1-1/T))\}$.

71 71
2nd Gumbel

Gumbel 00 plot$\mathrm{Q}\mathrm{Q}$

3rd GumbeI

plot$0\mathrm{O}$plot.
$.\dot{8\S\dot{\ddagger}}$$.\mathrm{s}^{\mathrm{s}}\dot{\xi i}$$.\dot{8}\in;$

-’

0

’

$t$

’

$-\prime S$$\triangleleft \mathrm{s}$図 6.

$\mathrm{o}s\prime \mathrm{D}$荻 uQdd荻 uQd\’eGumbel モデルの下での$\mathrm{Q}\mathrm{Q}$$u\cdot u[]$plot.

First -Second

Second -Third.
$|1$$8\ovalbox{\tt\small REJECT}_{}$.

.

$\mathrm{E}[] \mathrm{K}[] u\mathrm{u}\mathrm{d}\mathrm{Q}_{[][]\prime}\mathfrak{g}.$$\epsilon \mathrm{K}\mathrm{m}\mathrm{n}\mathrm{u}\mathrm{d}\mathrm{Q}\mathrm{u}\cdot \mathfrak{n}\mathrm{r}\cdot$

図 7.

Gumbel モデルの下での指数確率紙.

52GEV モデル
このモデルの下で, 上位 $r=1,2,3$ 個のデータを用いた場合の最尤推定値は, 次の様になった :$r=1$ ,$\hat{\mu}$,$\hat{\sigma}_{r}$$\hat{\xi}_{t}$2,0.425 0.426
0.035 0.034-0.378 -0.156

3
0.432
0.036
-0.280

最大値の分布は上限のある Weibull 分布と見なせる.
の決定のためにplot を書かせたのが図 8, 9, 10 で, 指数確率紙が図 11 である. これらから, このモデルの下でも Gumbel モデルの場合と同様に上位 2 個までのデータが使えると言推定値から,$r$$\xi<0$ の可能性が強く,$\mathrm{P}\mathrm{P}$える.

従って, この場合の $T$-return level の推定は$\hat{q}_{2}(T)=0.426+0.034\{(-\log(1-1/T))^{0.156}-1\}/(-0.156)$とすればよい. また, 最大孔食深さの上限値は 0646 と推定される.

72$\mathrm{P}\mathrm{P}$pIot. 1st 0. S..

$\mathrm{r}=1$.

1$\mathrm{P}\mathrm{P}$plot 1st$0l0$.

$0\mathrm{J}$or$’\rho$

$0A$

図 9.

plot, $r=1$ .

$\mathrm{P}\mathrm{P}$

p 化 t 1 社屋歌..

$\mathrm{r}=3$I.

$0\mathrm{J}$First and Second.

$\mathrm{u}\cdot u$.

$2\mathrm{n}\mathrm{d}\mathrm{O}.\mathrm{S}" \mathrm{r}--3$$\mathrm{P}\mathrm{P}$$\mathrm{P}\mathrm{P}\mathrm{p}\mathrm{b}\mathrm{t}$.

$\infty 0\mathfrak{n}u\mathrm{n}[]-$

3rd 屋歌..

.

$\mathrm{r}=3$$\mathrm{P}\mathrm{P}$$\mathrm{u}\cdot u$

plot, $r=3$ .

First and Second.

$\ulcorner-3$.

Second and Third.

$\ulcorner-3$$\xi$$.\xi$$.\xi$

plot, $r=2$ .

$.\cdot\xi$

u.d\’e$\ulcorner-2$$\mathrm{O}.\mathrm{S}..r=\mathit{2}$鋪I
GEV モデノレの下での

pIot 2nd$\prime D$

GEV モデ) の下での

憶憶 pIot$u_{9}u$

図 11.

.

$0l0*$荻$uu$

GEV モデルの下での

図 10.

$\mathrm{P}\mathrm{P}$$.\cdot\xi$$\mathrm{u}\mathfrak{g}\mathrm{d}d$$\mathrm{P}\mathrm{P}$.

$.\cdot\xi$$0\mathrm{D}\mathrm{o}z$

図 8.

$\mathrm{O}.\mathrm{S}..\mathrm{r}=2$.

$\mathrm{E}[][]\prime \mathrm{u}[] \mathrm{Q}\mathrm{u}[] \mathrm{U}\cdot$

GEV モデルの下での指数確率紙, $r=2,$ $r=3$ .

0

’

2

$\epsilon[][][] u\mathrm{Q}[[]\prime u\cdot\cdot$

$

.

.

73 73
付録
Al. Gumbel モデルの場合の Fisher 情報量, Smith (1986)
サイズ

の上位

$n$

$r$

個の同時分布の Fisher 情報行列は(Al.l)$I_{F}( \theta)=\frac{n}{\sigma^{2}}(\begin{array}{ll}r -B_{r}-B_{r} C_{r}\end{array})$となる. ただし,$\theta=(\mu, \sigma)$,$B_{r}=r\psi(r+1)$,$C_{r}=r\{\psi^{2}(r+1)+\psi’(r+1)+1\}$

.

(AL2)

逆行列は
$\frac{\sigma^{2}}{n(rC_{r}-B_{r}^{2})}(\begin{array}{ll}C_{r} B_{r}B_{r} r\end{array})$と表される.
ここで,$T$

これが, 最尤推定量$\hat{\theta}_{r}=(\hat{\mu}_{r},\hat{\sigma}_{r})$(A1.3)の漸近分散行列である.

-return level $q(T)$ の推定は,$\hat{q}_{r}(T)=\hat{\mu}_{r}+\hat{\sigma}_{r}g(T)$とすればよかった. 推定量$AV(\hat{q}_{r}(T))$$\hat{q}_{r}(T)$,$g(T)=-\log(-\log(1-1/T))$の漸近分散は$=$$AV(\hat{\mu}_{r})+(g(T))^{2}AV(\hat{\sigma}_{r})+2g(T)\mathrm{A}\mathrm{C}\mathrm{o}\mathrm{v}(\hat{\mu}_{r},\hat{\sigma}_{r})$$=$$\frac{\sigma^{2}}{n(rC_{r}-B_{r}^{2})}[C_{r}+(g(T))^{2}r+2g(T)B_{r}]$(A1.4)となる.

A2. GEV モデルの場合の Fisher 情報量? Tawn (1988)
GE 好皀妊襪両豺腓 Fisher 情報行列はかなり複雑になる.
サイズの上位 個の同時分布の Fisher 情報行列を
$n$

$r$

$I_{F}(\theta)=n\{\begin{array}{l}E_{\theta}\{\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\mu\partial\sigma}\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\mu\partial\xi}\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\mu\partial\sigma}\}E_{\theta}\{\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\sigma\partial\xi}\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\mu\partial\xi}\}E_{\theta}\{-\frac{\partial^{2}l}{\partial\sigma\partial\xi}\}E_{\theta}\{\}\end{array}\}$とする,

ここで $\theta=(\mu, \sigma, \xi)$. このとき,$E_{\theta} \{-\frac{\partial^{2}l}{\partial\mu^{2}}\}=\frac{(1+\xi)^{2}}{\sigma^{2}(1+2\xi)}G(2)$,$E_{\theta} \{-\frac{\partial^{2}l}{\partial\mu\partial\sigma}\}=\frac{1}{\sigma^{2}\xi(1+2\xi)}\{(1+2\xi)G(1)-(1+\xi)^{2}G(2)\}$,(A2.1)

74$E_{\theta} \{-\frac{\partial^{2}l}{\partial\mu\partial\xi}\}=\frac{1}{\sigma\xi^{2}(1+2\xi)}[(1+\xi)^{2}G(2)-(1+2\xi)G(1)\{\xi\psi(r+\xi+1)+\frac{1+\xi+\xi^{2}}{1+\xi}\}]$,,$E \theta\{-\frac{\partial^{2}l}{\partial\sigma^{2}}\}=\frac{1}{\sigma^{2}\xi^{2}(1+2\xi)}\{r(1+2\xi)-2(1+2\xi)G(1)+(1+\xi)^{2}G(2)\}$$E_{\theta} \{-\frac{\partial^{2}l}{\partial\sigma\partial\xi}\}=\frac{1}{\sigma\xi^{3}(1+2\xi)}[(1+2\xi)G(1)\{\xi\psi(r+\xi+1)+\frac{1+(1+\xi)^{2}}{1+\xi}\}-r\xi(1+2\xi)\psi(r+1)$$-(1+\xi)^{2}G(2)-r(1+2\xi)]$ ,$E \theta\{-\frac{\partial^{2}l}{\partial\xi^{2}}\}=\frac{1}{\xi^{4}(1+2\xi)}[(1+\xi)^{2}G(2)-2(1+2\xi)G(1)\{\xi\psi(r+\xi+1)+\frac{1+\xi+\xi^{2}}{1+\xi}\}$$+r(1+2\xi)\{1+2\xi\psi(r+1)+\xi^{2}[1+\psi’(r+1)+\psi^{2}(r+1)]\}]$ ,ただし, $G(j)=G(j;r, \xi)=\Gamma(r+j\xi+1)/\Gamma(r),$ $j=1,2$ .

参考文献
Hosking, J. R. M., Wallis, J. R. and Wood, E. F. (1985). Estimation of the generalized extremevalue disribution by the method of probability-weighted moments. Technometrics, 27,
251-261.
Nagaraja, H. N. (1982). Record values and extreme value distributions. J. Appl. Probab., 19,233-239.
Scarf, P. A., Cottis, R. A. and Laycock, P. J. (1992). Extrapolation of extreme pit depths inspace and time using the deepest pit depths. J. Electrochem. Soc., 139, 2621-2627.
$r$

Smith, R. L. (1986). Extreme value theory based on the86, 27-43.

$r$largest annual events. J. Hydrol.,Tawn, J. A. (1988). An extreme value theory model for dependent observations. J. Hydrol.,101, 227-250.

Weissman, I. (1978). Estimation of parameters and large quantiles based on theobservations. J. Arner. Statist. Assoc. 73, 812-815.

$k$largest