上位 r 個の観測値に基づく確率点の推定

2004 c

上位 r 個の観測値に基づく確率点の推定

高橋 倫也

¹

²

（受付

2003

7

31

2004

1

27

要 旨

いくつかの単位領域（または単位期間）のそれぞれで測定したデータを利用して，上側微小確 率点を推定する．単位領域ごとの最大値データのみを用いるのが古典的極値解析の手法であっ た．ここでは，最大値だけでなく上位

r

r

キーワード： 漸近相対効率，漸近分散，一般極値分布，Gumbel分布．

1.

与えられた領域（または期間）での最大値を推定するために，いくつかの単位領域（または単 位期間）ごとの最大値（極値）データを利用するのが古典的極値解析の手法である．その推定精 度を上げる方法として「上位

r

r

上位

r

Tawn

1988

Weissman

r

定を初めて議論した．彼は母集団分布が

Gumbel

Gumbel

を詳しく調べた．Smith（1986）は

Gumbel

r

N

10

Dupuis

Smith

Tawn

1988

Smith

1986

GEV

Lowestoft

この手法の他の分野への応用には次のような研究がある．

Robinson and Tawn

Smith

（1997）は

1993

3000 m

1

5–1–1

2

2–19–1

た可能性があるか，すなわち王の記録はそれまでの競技者の能力をはるかに超えたものかどう かの議論をしている．

Robinson and Tawn

1995

5

GEV

1500 m

5

タは

90％ 信頼区間におさまり，かなり稀少なデータだが否認できないという結論を述べた．

Smith

タは明らかに外れ値であると断定した．

Strand and Boes

1998

5

Gumbel

al.

1992

GEV

Scarf and Laycock

Coles and Walshaw

Palutikof et al.

上位

r

Nagaraja

1982

Scarf

1993

以下，2節で極値理論から導かれる上位

r

3

r

4

r

r

PWM

2.

一般極値（

generalized extreme value

(2.1) G _ξ (z) = exp[−(1 + ξz) ^−1/ξ ] , 1 + ξz > 0 (ξ ∈ R)

とする．ここで

G _ξ

ξ < 0

Weibull

ξ = 0

G ₀ (z) = lim _ξ→0 G _ξ (z) = exp( − e ^−z )

Gumbel

0

Fr´ echet

分布

F

Y ₁ , Y ₂ , . . . , Y _n

Y _1:n ≥ Y _2:n ≥ · · · ≥ Y _n:n

F

G _ξ

a n > 0，b n ∈ R，

n = 1, 2, . . .

n→∞ lim P

Y _1:n − b _n a _n ≤ z

= G _ξ (z) , ∀ z ∈ R .

r

P

Y _1:n − b _n

a n ≤ z ₁ , Y _2:n − b _n

a n ≤ z ₂ , . . . , Y _r:n − b _n a n ≤ z r

, z ₁ ≥ z ₂ ≥ · · · ≥ z r

は同時密度関数

g _{ξ, 12···r} (z ₁ , z ₂ , . . . , z _r ) = g _ξ (z ₁ ) · · · g _ξ (z _r−1 )

G _ξ (z ₁ ) · · · G _ξ (z _r−1 ) g _ξ (z _r ) , g _ξ (z) = dG _ξ (z)/dz

G _{ξ, 12···r}

David and Nagaraja

2003

10.6

ここで，

(2.2) g _{ξ, 12···r} (z ₁ , z ₂ , . . . , z r ) =

exp

−

r

j=1

z _j − e ^{−z r} , ξ = 0

r

j=1

(1 + ξz j ) ^−1/ξ−1 exp[−(1 + ξz r ) ^−1/ξ ] , ξ = 0

確率ベクトル

(Z ₁ , Z ₂ , . . . , Z r )

G ξ, 12···r

j

≥ 1

G ξ,j

は

r

(2.3) G _ξ,j (z) =

j−1

k=0

exp( − kz − e ^−z )/k! , ξ = 0

j−1

k=0

(1 + ξz) ^−k/ξ exp[ − (1 + ξz) ^−1/ξ ]/k! , ξ = 0

g _ξ,j

(2.4) g _ξ,j (z) =

exp(−jz − e ^−z )/Γ(j) , ξ = 0

(1 + ξz) ^−j/ξ−1 exp[−(1 + ξz) ^−1/ξ ]/Γ(j) , ξ = 0

定理

1. (Z ₁ , Z ₂ , . . . , Z r )

1 , W ₂ , . . .

S _j =

j

k=1 W _k

（I）

Gumbel (ξ = 0)

（

I.a

{ Z _j , j ≥ 1 } = ^d {− log S _j , j ≥ 1 } .

（I.b）

(Z ₁ , Z ₂ , . . . , Z _r )

Exp

r

j(Z j − Z _j+1 ) , j = 1, 2, 3, . . .

Exp(1)

（I.c）

E ₀ (Z j ) = −ψ(j) , V ₀ (Z j ) = ψ (j) .

（I.d）

Cor ₀ (Z j , Z _j+1 ) =

ψ (j + 1)/ψ (j) .

（

II

GEV (ξ = 0)

（II.a）

{ Z _j , j ≥ 1 } = ^d { (S _j ^−ξ − 1)/ξ, j ≥ 1 } .

（II.b）

(Z ₁ , Z ₂ , . . . , Z r )

ξ

Pareto

r

j

ξ log 1 + ξZ j

1 + ξZ _j+1 , j = 1, 2, 3, . . .

Exp(1)

（II.c）

E _ξ (Z _j ) = Γ(j − ξ) − Γ(j)

ξΓ(j) , V _ξ (Z _j ) = Γ(j − 2ξ)Γ(j) − Γ ² (j − ξ) ξ ² Γ ² (j) .

（II.d）

Cor _ξ (Z j , Z _j+1 ) = 1 j − ξ

Γ(j + 1 − 2ξ)Γ(j + 1) − Γ ² (j + 1 − ξ) Γ(j − 2ξ)Γ(j) − Γ ² (j − ξ) .

証明. （

I

I.a

( − log S ₁ , . . . , − log S _r )

(W ₁ , . . . , W _r )

g _0,12···r

（I.b）

(Z ₁ − Z ₂ , 2(Z ₂ − Z ₃ ), . . . , (r − 1)(Z _r−1 − Z r ))

e ^{−y 1} e ^{−y 2} · · · e ^{−y r−1} , y ₁ , y ₂ , . . . , y _r−1 > 0

となり，

j(Z _j − Z _j−1 )

Exp

1

（I.c） 式（2.4）から求まる．

（I.d）

(Z j , Z _j+1 )

g _0,jj+1 (z _j , z _j+1 ) = 1

Γ(j) exp( − jz _j )g ₀ (z _j+1 ) , z _j ≥ z _j+1

E ₀ (Z j Z _j+1 ) = 1

j ² Γ(j) (jΓ (j + 1) − Γ (j + 1)) = ψ ² (j + 1) + ψ (j + 1) − 1

i ψ(j + 1)

（

II

I

I

（

II.a

P

S _j ^−ξ − 1

ξ ≤ z _j , j = 1, 2, . . . , r = P ( − log S _j ≤ log(1 + ξz _j ) ^1/ξ , j = 1, 2, . . . , r) .

r

j=1

(1 + ξz _j ) ⁻¹ g _0,12···r (log(1 + ξz ₁ ) ^1/ξ , . . . , log(1 + ξz _r ) ^1/ξ )

g _ξ,12···r (z ₁ , . . . , z _r )

（II.b） 形状パラメータ

ξ

Pareto

P (y) = 1 − (1 + ξy) ^−1/ξ , 1 + ξy > 0

（

ξ ≥ 0

y > 0，ξ < 0

0 < y < − 1/ξ

n

V _1:n ≥ V _2:n ≥ · · · ≥ V _n:n

とすると，

V _j:n = ^d U _n−j+1:n ^−ξ − 1

ξ , j = 1, 2, . . . , n

と表される．ただし，1

> U _1:n ≥ U _2:n ≥ · · · ≥ U _n:n > 0

U(0, 1)

n

V _j:n

II.a

(Z ₁ , Z ₂ , . . . , Z _r )

ξ

Pareto

r

（

I.a

I.b

j log(S _j+1 /S _j )

Exp

1

II.a

j

ξ log 1 + ξZ j

1 + ξZ _j+1

= d j log S _j+1

S _j , j = 1, 2, 3, . . .

（II.c） 式（2.4）から求まる．

（II.d）

(Z j , Z _j+1 )

g _ξ,jj+1 (z j , z _j+1 ) = 1

Γ(j) (1 + ξz j ) ^−j/ξ−1 g _ξ (z _j+1 ) , z j ≥ z _j+1

となる．これから

E ξ (Z j Z _j+1 ) = 1 Γ(j + 1)

1

ξ ² [Γ(j + 1 − 2ξ) − 2 Γ(j + 1 − ξ) + Γ(j + 1)]

+ 1

ξ(j − ξ) [Γ(j + 1 − 2ξ) − Γ(j + 1 − ξ)]

を求め，（

II.c

2

注

1.

I.a

II.a

Nagaraja

1982

I.b

Weissman

1978

II.b

Tawn

1988

注

2.

ξ → 0

注

3. Cor ₀ (Z _j , Z _j+1 )

j

1

Cor ₀ (Z ₁ , Z ₂ ) = 0.626， Cor ₀ (Z ₂ , Z ₃ ) = 0.783， Cor ₀ (Z ₃ , Z ₄ ) = 0.848

= 0

注

4.

ψ(n) = − γ +

n−1

i=1

1

i , ψ (n) = π ² 6 − ⁿ⁻¹

i=1

1

i ² , n = 1, 2, . . .

= 0.57721566...

Euler

図

1

ξ = − 0.4，0，0.4

Z ₁

2

3

図

2

3

4

ξ = − 0.4

0

0.4

(Z ₁ , Z ₂ )

g _ξ,12 (z ₁ , z ₂ ) =

(1 + ξz ₁ ) ^−1/ξ−1 (1 + ξz ₂ ) ^−1/ξ−1 exp{−(1 + ξz ₂ ) ^−1/ξ } , ξ = 0

exp(−z 1 − z ₂ − e ^{−z 2} ) , ξ = 0

z ₁ ≥ z ₂

3.

一般極値分布

G _ξ

µ，σ

= (µ, σ)

(µ, σ, ξ)

T -return level

T

q(T )，

G _ξ

q(T ) − µ σ

= 1 − 1 T ,

(3.1) q(T ) =

µ + σ {− log

− log(1 − 1/T)

} , ξ = 0 µ + σ{

− log(1 − 1/T )

−ξ − 1}/ξ , ξ = 0

図

1. Weibull ( ξ = − 0 . 4)， Gumbel ( ξ = 0)， Fr´ echet ( ξ = 0 . 4)

j (= 1 , 2 , 3)

図

2. Weibull

( ξ = − 0 . 4)

1， 2

3.1

上位

r

(X ₁ , X ₂ , . . . , X r )

1

σ ^r g _0,12···r

x ₁ − µ

σ , . . . , x _r − µ σ

図

3. Gumbel

( ξ = 0)

1， 2

図

4. Fr´ echet

( ξ = 0 . 4)

1， 2

で，GEVモデルの場合は

1 σ ^r g _ξ,12···r

x ₁ − µ

σ , . . . , x r − µ σ

となる．

データは独立な

x ₁ ≥ x ₂ ≥ · · · ≥ x _r

n

x _i1 ≥ x _i2 ≥ · · · ≥ x ir ; i = 1, . . . , n

の

n × r

3.2

3.2 Gumbel

l(µ, σ) = − nr log σ −

n

i=1

r

j=1

x ij − µ σ

+ exp

− x ir − µ σ

となる．パラメータ

= (µ, σ)

∂

∂µ l(µ, σ) = −

n

i=1

− r σ + 1

σ exp

− x ir − µ σ

= 0 ,

∂

∂σ l(µ, σ) = − nr σ +

n

i=1

r

j=1

x ij − µ σ ²

− x ir − µ σ ² exp

− x ir − µ σ

= 0 .

σ

h _r (σ) :=

n

i=1

x ir − x ¯ ¯ r

σ + 1

exp

− x ir − x ¯ ¯ r

σ

= 0 , x ¯ ¯ _r = 1 nr

n

i=1 r

j=1

x _ij

h _r (σ) = 1 σ

n

i=1

x ir − ¯ ¯ x r

σ

2

exp

− x ir − x ¯ ¯ r

σ

を用いて，ニュートン法で

σ _r

µ _r = −

σ _r log

1 nr

n

i=1

exp

− x _ir

σ _r

が求まる．この

r = ( µ

r , σ

r )

r

= (µ, σ)

T -return level q(T )

(3.2) q

r (T ) =

µ r +

σ r {− log(− log(1 − 1/T ))}

とすればよい．また，推定値

µ

_r

σ _r

q

_r (T )

A.1

r

r

{ (x _i1 , x _i2 , . . . , x ir ), i = 1, 2, . . . , n }

G _{0, 12...r}

r

上位

j

X j

G _0,j

(z − µ)/σ

を持つ．したがって，

U _ij = G _0,j

X ij − µ σ

により，X

ij

U(0, 1)

U ij

このことから，次の周辺分布の適合を見る

r

PP plot

r

µ

r

σ r

j(= 1, 2, . . . , r)

{ x _1j , x _2j , . . . , x _nj }

u _ij = G _0,j ((x _ij − µ

_r )/

σ _r ), i = 1, 2, . . . , n

u ij

u _(1)j ≥ u _(2)j ≥ · · · ≥ u _(n)j

1 − i

n + 1 , u _(i)j

, i = 1, 2, . . . , n

をプロットし，

r

PP plot

r

r

PP plot

r

一方，次の決定法も考えられる．

QQ plot： r

r

µ _r

σ _r

j(= 1, 2, . . . , r)

q _(i)j = G ⁻¹ _0,j (1 − i/(n + 1))

j

{x 1j , x _2j , . . . , x nj }

x _(1)j ≥ x _(2)j ≥ · · · ≥ x _(n)j

( µ

_r +

σ _r q _(i)j , x _(i)j ) , i = 1, 2, . . . , n

をプロットし，

r

QQ plot

r

r

QQ plot

r

≥ 2

q _(i)j

これら

PP plot，QQ plot

チェックする方法として，定理

1

指数確率紙：

j = 1, 2, . . .

{ x _ij − x _ij+1 , i = 1, 2, . . . , n }

j(< r)

r

したがって，r として

PP plot，QQ plot

PP plot

QQ plot

Smith

1986

Coles

2001

Smith

補助変量を含む場合

Smith

1986

タが年

(i = 1, 2, . . . , n)

µ _i = α + β

n i , σ _i = σ , i = 1, 2, . . . , n

l(α, β, σ) = −nr log σ −

n

i=1

r

j=1

x ij − α − β i/n σ

+ exp

− x ir − α − β i/n σ

と書ける．これを数値計算で最大化し最尤推定値

( α

r , β

r ,

σ r )

このとき

i

T -return level

α r + β _r

n i +

σ r {− log

− log(1 − 1/T)

}

r

QQ plot

µ r

µ _i = α

_r + β r

n i , i = 1, 2, . . . , n

パラメータが補助変量のもっと複雑な関数の場合も同様にできる．

3.3 GEV

l(µ, σ, ξ) = −nr log σ −

n

i=1

1 ξ + 1

r

j=1

log

1 + ξ

x ij − µ σ

+

1 + ξ

x ir − µ σ

−1/ξ

となる．パラメータ

= (µ, σ, ξ)

尤度方程式は簡単にはならない．ニュートン法で連立非線形の尤度方程式を解かなければな らないが，初期値としては極値データから

PWM

A.3

r = ( µ

r , σ

r , ξ

r )

-return level q(T )

q r (T ) = µ

r +

σ r {(− log(1 − 1/T )) ^{− ξ}

^r − 1}/ ξ

r

とすればよい．また，推定値

µ _r

σ

_r

ξ

_r

q

_r (T )

A.2

Fisher

r

r

{ (x _i1 , x _i2 , . . . , x _ir ), i = 1, 2, . . . , n }

G _{ξ, 12...r}

r

上位

j

X _j

G _ξ,j

(z − µ)/σ

を持つ．したがって，

U ij = G ξ,j

X _ij − µ σ

により，

X _ij

U(0, 1)

U _ij

このことから，次の周辺分布の適合を見る

r

PP plot： r

r

µ

_r

σ

_r

ξ

_r

j(= 1, 2, . . . , r)

{x 1j , x _2j , . . . , x nj }

u ij = G ξ r ,j ((x ij − µ

r )/

σ r ), i = 1, 2, . . . , n

u _ij

u _(1)j ≥ u _(2)j ≥ · · · ≥ u _(n)j

1 − i

n + 1 , u _(i)j

, i = 1, 2, . . . , n

をプロットし，r個の

PP plot

r

r

PP plot

r

また，次の方法も考えられる．

QQ plot： r

r

µ

_r

σ

_r

ξ _r

2, . . . , r)

q _(i)j = G ⁻¹

ξ _r ,j (1 − i/(n + 1))

j

{ x _1j , x _2j , . . . , x _nj }

x _(1)j ≥ x _(2)j ≥ · · · ≥ x _(n)j

( µ

r +

σ r q _(i)j , x _(i)j )

i = 1, 2, . . . , n

をプロットし，

r

QQ plot

r

r

QQ plot

r

≥ 2

q _(i)j

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

1

指数確率紙：

r

r

µ _r

σ

_r

ξ

_r

j(=

1, 2, . . . , r − 1)

j + 1

{ (x _ij , x _ij+1 ), i = 1, 2, . . . , n }

j

ξ _r log σ

r + ξ

r (x ij −

µ r )

σ _r + ξ

_r (x _ij+1 − µ

_r ) , i = 1, 2, . . . , n

を求めて指数確率紙にプロットし，

r − 1

r

r − 1

r

したがって，

r

PP plot

QQ plot

PP plot，QQ plot

Tawn

Tawn

1988

補助変量を含む場合

Tawn

µ

節の

Gumbel

Scarf et al.

(x ij , t i ) , i = 1, . . . , n ; j = 1, . . . , r i ; (x _i1 ≥ · · · ≥ x ir _i ) ,

すなわち，時刻

t _i

r _i

{ (x _i1 , . . . , x _{ir i} )

i = 1, 2, . . . , n }

t

µ t = µt ^β , σ t = σt ^β

l(µ, σ, β, ξ) = −

n

i=1

r _i

log σ + β log t _i

+

1 ξ + 1

r _i

i=1

log

1 + ξ

x _ij t ^−β _i − µ σ

−

1 + ξ

x ir _i t ^−β _i − µ σ

−1/ξ

と書ける．これを数値計算で最大化し最尤推定値

(

µ, σ,

β,

ξ)

r _i (i = 1, 2, . . . , n)

0 = min { r ₁ , . . . , r _n }

r ₀

r

3

µ

_r

σ _r

µ t ^β _i

,

σ t ^β _i

, i = 1, 2, . . . , n .

この方法で

r = r ^∗ ( ≤ r ₀ )

^∗

r _i

1

r _i (i = 1, . . . , n)

4.

ここでは，パラメータ

(µ, σ, ξ)

上位

r( ≥ 2)

T -return level q(T )

)

r

まず，

Gumbel

推定量

q

r (T )

A.1

AV (

q r (T )) = σ ²

n(rC _r − B _r ² ) (C r + 2g(T )B r + (g(T )) ² r) , g(T ) = − log(− log(1 − 1/T ))

q ₁ (T )

q

r (T )

(4.1) e ₀ ( q

r (T ),

q ₁ (T )) = AV ( q

₁ (T ))

AV ( q

r (T )) = (rC _r − B ² _r )(C ₁ + 2g(T )B ₁ + (g(T )) ² )

(C ₁ − B ₁ ² )(C r + 2g(T )B r + (g(T )) ² r)

と表される．この漸近相対効率は

r

T

(4.2) e ₀ (r, T ) = e ₀ ( q

_r (T ) , q

₁ (T ))

r = 2, 3, . . . , 10

T =100，1,000，10,000，100,000

e ₀ (r, T )

1

1

r = 3，T = 10, 000

e ₀ (3, 10, 000) = 1.995

level

3

50

100

等しい漸近分散，を持つと言える．また，漸近相対効率は

r

漸近相対効率

e ₀ (r, T )

命題

1.

T

（I.a）

e ₀ (r, T )

r

∞

（

I.b

g(T ) > 0

e ₀ (r, T )/r

r

（II）

r

0 (r, T)

T

r

r

ψ (r + 1) + 1

ψ (2) + 1

証明. 以下，狭義単調増加（減少）関数を増加（減少）関数と言う．

（I）（I.a） 付録

A.1

Σ _r−1 (

) − Σ r (

)

AV ( q

_r (T )) < AV ( q

_r−1 (T )) .

e ₀ (r, T ) > e ₀ (r − 1, T ) , r = 2, 3, . . .

r ²

r → ∞

分子

= c ₁ (ψ (r + 1) + 1) → c ₁ ,

分母

= c ₂

C _r

r ² + 2g(T ) B _r

r ² + (g(T )) ² r

→ 0 ,

ただし，c

1

2

r→∞ lim e ₀ (r, T ) = ∞ .

表

1. q

1 ( T )

q r ( T )

e 0 ( r, T )．

（I.b） ここで

r

σ ² Σ _r (

) = 1 ψ (r + 1) + 1

ψ ² (r + 1) + ψ (r + 1) + 1 ψ(r + 1)

ψ(r + 1) 1

より，

1

σ ² {rΣ r (

) − (r − 1)Σ _r−1 (

)} =

a b b c

とおくと，

a = ψ ² (r + 1) + ψ (r + 1) + 1

ψ (r + 1) + 1 − ψ ² (r) + ψ (r) + 1 ψ (r) + 1 , b = ψ(r + 1)

ψ (r + 1) + 1 − ψ(r)

ψ (r) + 1 , c = 1

ψ (r + 1) + 1 − 1 ψ (r) + 1

である．ψは増加関数で

ψ

b, c > 0．一方，a > 0

a

よって，

g(T ) > 0

1 g(T )

{rΣ r (

) − (r − 1)Σ _r−1 (

)}

1 g(T )

= σ ² {a + 2b g(T ) + c g ² (T )} > 0 .

r AV ( q

r (T )) > (r − 1)AV ( q

_r−1 (T ))

e ₀ (r, T )

r < e ₀ (r − 1, T )

r − 1 , r = 2, 3, . . .

（II）

r

)

T

g

e(g) = c g ² + 2B ₁ g + C ₁ rg ² + 2B r g + C r

, c

. g

e(g)

e(g)

g

n(g)

n(g) = (B r − rB ₁ )g ² + (C r − rC ₁ )g + (B ₁ C r − B r C ₁ )

となる．ここで，r

≥ 2

g ²

g

B r − rB ₁ = r{ψ(r + 1) − ψ(2)} > 0 ,

C _r − rC ₁ = r { ψ ² (r + 1) + ψ (r + 1) + 1 } − r { ψ ² (2) + ψ (2) + 1 }

= r { ψ ² (r + 1) − ψ ² (2) + ψ (r + 1) − ψ (2) }

= r

ψ(2) +

r

i=2

1 i

2

− ψ ² (2) −

r

i=2

1 i ²

> 0 .

したがって，

g

n(g) > 0

de(g)/dg > 0

e ₀ (r, T )

T

1

r

T ( ≥ 100)

e ₀ (r, T )