New アクチュアリー「数学」演習sugiura/2016/act_math2016b2.pdf · 2018. 9. 14. ·...
52
アクチュアリー「数学」演習 杉浦 誠 最終変更日: 2018 年 9 月 14 日 目次 1 回帰分析 1 1.1 回帰直線 (単回帰) ......................................... 1 1.2 重回帰 ............................................... 2 1.3 非線形回帰 ............................................. 3 1.4 確率分布の前提を置いた回帰モデルの分析 ............................ 4 1.5 統計の復習 1 正規母集団と二項母集団 ............................. 7 2 時系列解析 8 2.1 時系列に現れる確率過程と用語の定義 ............................... 8 2.2 AR(p)(p 次の自己回帰モデル, Auto-regressive Model) ..................... 9 2.3 MA(q)(q 次の移動平均モデル, Moving-average Model) .................... 11 2.4 ARMA(p, q) (Autoregressive Moving-average Model) ..................... 12 2.5 時系列モデルに基づく予測 .................................... 13 2.6 統計の復習 2 順序統計量 .................................... 14 3 確率過程 17 3.1 マルコフ連鎖とマルチンゲール .................................. 17 3.2 ポアソン過程 ............................................ 20 3.3 ブラウン運動 ............................................ 21 4 シミュレーション 22 4.1 確率変数を生成する技法 ...................................... 22 4.2 分散減少法 ............................................. 26 4.3 統計の復習 3 適合度、独立性の検定 .............................. 29 5 損保数理に関する確率統計の話題から 32 5.1 最尤推定量の漸近挙動 ....................................... 32 5.2 極値問題 .............................................. 38 5.3 安定分布 .............................................. 46 1
Transcript of New アクチュアリー「数学」演習sugiura/2016/act_math2016b2.pdf · 2018. 9. 14. ·...
アクチュアリー「数学」演習
杉浦 誠
最終変更日: 2018年 9月 14日
目次
1 回帰分析 1
1.1 回帰直線 (単回帰) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
1.2 重回帰 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2
1.3 非線形回帰 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3
1.4 確率分布の前提を置いた回帰モデルの分析 . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4
1.5 統計の復習 1 正規母集団と二項母集団 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2 時系列解析 8
2.1 時系列に現れる確率過程と用語の定義 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 8
2.2 AR(p) (p次の自己回帰モデル, Auto-regressive Model) . . . . . . . . . . . . . . . . . . . . . 9
2.3 MA(q) (q 次の移動平均モデル, Moving-average Model) . . . . . . . . . . . . . . . . . . . . 11
2.4 ARMA(p, q) (Autoregressive Moving-average Model) . . . . . . . . . . . . . . . . . . . . . 12
2.5 時系列モデルに基づく予測 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.6 統計の復習 2 順序統計量 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 14
3 確率過程 17
3.1 マルコフ連鎖とマルチンゲール . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
3.2 ポアソン過程 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
3.3 ブラウン運動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
4 シミュレーション 22
4.1 確率変数を生成する技法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
4.2 分散減少法 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
4.3 統計の復習 3 適合度、独立性の検定 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29
5 損保数理に関する確率統計の話題から 32
5.1 最尤推定量の漸近挙動 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 32
5.2 極値問題 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.3 安定分布 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
1
これは 2014年度後期に情報理論 IIとして行うアクチュアリー試験「数学」用の講義ノートです。教科書・
参考書として以下を用いています。
• 日本アクチュアリー会編 モデリング 日本アクチュアリー会• 藤田岳彦 著 確率・統計・モデリング問題集 日本アクチュアリー会• 藤田岳彦 著 弱点克服大学生の確率・統計 東京図書, 2010
• 黒田耕嗣 著 生保年金数理 培風館, 2007
• 岩沢宏和 黒田耕嗣 著 損害保険数理 (アクチュアリー数学シリーズ 4), 日本評論社, 2015
• 国沢清典編 確率統計演習 2 統計 培風館, 1966
• 稲垣宣生 著 数理統計学 改訂版 裳華房, 2003
• 小寺平治 著 明解演習 数理統計 共立出版, 1986
• E.L. Lehmann, George Casella: Theory of Point Estimation, Second Edition, Springer, 1998
• S.I. Resnick: Extreme Values, Regular Variation and Point Processes, Springer, 1987
• 高橋 倫也, 志村 隆彰: 極値統計学 (ISMシリーズ:進化する統計数理), 近代科学社, 2016
• Breiman, L.: Probability, Addison-Wesley, 1968. (Classics in applied mathematics, 7, Society for
Industrial and Applied Mathematics, 1992. Reprint版)
• Durrett, R.: Probablity Theory and Examples, 4th ed., Cambridge University Press, 2010.
教科書・参考書は今後増えていく予定です。
2
•「藤田岳彦 著 確率・統計・モデリング問題集 日本アクチュアリー会」に従って述べていく。
1 回帰分析
1.1 回帰直線 (単回帰)
2種類のデータの観測値 (xi, yi), (i = 1, 2, · · · , n) が与えられているとする。
x =1
n
n∑i=1
xi, (データの平均)
sx2 =
1
n
n∑i=1
(xi − x)2 =1
n
n∑i=1
xi2 − x2 = x2 − x2, (データの分散)
sxy =1
n
n∑i=1
(xi − x)(yi − y) =1
n
n∑i=1
xiyi − x y = xy − x y, (データの共分散)
rxy =sxysxsy
, (データの相関係数) ただし、sx =√sx2
などがデータの性質や関係を表す基本的な量である。以下の性質があった。
• −1 ≤ rxy ≤ 1.
• rxy = 1 (−1) ⇐⇒ ある定数 a > 0 (a < 0) が存在し ∀iに対して yi = axi + b.
• rxy ≒ 1 (−1)のとき、正の相関 (負の相関)が強いという。
• a, b, c, dを定数とし ac > 0のとき、rax+b,cy+d = rxy. (相関係数は単位のとりかたによらない。)
• 最小二乗法xi から予測される値 α + βxi と現実の値 yi との差の二乗
yi
xi
α+ βxi
y = α+ βx
x
y
O
の和 Q =n∑i=1
{yi − (α + βxi)}2 が最小となるように係数
α = α, β = β を定める:
0 =∂Q
∂α= −2
n∑i=1
(yi − (α+ βxi)
)= −2n(y − α− βx)
0 =∂Q
∂β= −2
n∑i=1
xi
(yi − (α+ βxi)
)= −2n(xy − αx− βx2)
これより正規方程式
α+ βx = y
αx+ βx2 = xy行列表示で
(1 x
x x2
)(α
β
)=
(yxy
)(1.1)
これを解いて β =−x y + xy
x2 − x2=sxys2x
= rxysysx
, α = y − βx = y − rxysysxxを得る。
この y = α + βxを xを説明変数、y を被説明変数とする回帰直線という。単に、xから y への回帰直線と
いうこともある。α = y − βx, β = rxysysxより、回帰直線は y − y = β(x− x)あるいは
y − ysy
= rxyx− xsx
つまり y の標準化 = 相関係数× xの標準化 (1.2)
と表されることに注意する。また、回帰直線は (x, y)を通ることを注意する。
問題 1.1 xから y への回帰直線が y = α1 + β1x, y から xへの回帰直線が x = α2 + β2y であるとする。
(1) β1β2 > 0のとき rxy, sy/sx を β1, β2 を用いて表せ。
(2) 更に、β1β2 = 1と仮定する。x, y を α1, α2, β1, β2 を用いて表せ。
1
• 決定係数
yi = α+ βxi を yi の内挿値、ei = yi − yi を残差という。このとき、n∑i=1
ei = 0,
n∑i=1
xiei = 0が成立する。実
際、α = y − βxと β =sxys2xに注意すれば
n∑i=1
ei =
n∑i=1
(yi − yi) = n(y − (α+ βx)
)= 0,
n∑i=1
xiei =
n∑i=1
xi(yi − yi) = n(xy − (αx+ βx2)
)= n
(xy − (y − βx)x− βx2
)= n
(sxy − βs2x
)= 0.
全変動、回帰変動、残差変動について以下の関係式が成り立つ:
全変動 ≡n∑i=1
(yi − y)2 =
n∑i=1
(ei + yi − y)2 =
n∑i=1
e2i + 2
n∑i=1
ei(yi − y) +n∑i=1
(yi − y)2
=
n∑i=1
e2i + 2
n∑i=1
ei(α+ βxi − y) +n∑i=1
(yi − y)2
=
n∑i=1
e2i +
n∑i=1
(yi − y)2 ≡ 残差変動+回帰変動.
決定係数 R2 = 1− 残差変動全変動
=回帰変動全変動
と定める。R2 が 1に近いほど回帰直線がデータによくあてはまっ
ている。また、以下のように R2 = rxy2 が示せる。
全変動 ≡n∑i=1
(yi − y)2 = ns2y
回帰変動 ≡n∑i=1
(yi − y)2 =
n∑i=1
{α+ βxi − (α+ βx)
}2
= β2n∑i=1
(xi − x)2 =(rxy
sysx
)2· ns2x.
問題 1.2 次に対し x, y, s2x, s2y, sxy, β, αと、全変動, 決定係数 R2, 回帰変動, 残差変動を求めよ。
(1) (xi, yi) = (i, i2) (i = 1, 2, . . . , n) ヒント:
n∑i=1
i4 =n(n+ 1)(2n+ 1)(3n2 + 3n− 1)
30(導けるかな?)
(2)
i 1 2 3 4 5 6
xi 1 2 2 4 5 5
yi 5 14 11 21 18 26
(3)
i 1 2 3 4 5 6 7 8
xi 3 6 8 9 6 7 3 5
yi 4 7 8 9 4 5 5 6
((2), (3)は電卓を
用い分数で表せ。)
1.2 重回帰
単回帰では説明変数が 1つだったが、ここでは 2個以上の場合を考える。簡単のため 2個として説明する。
データの観測値 (x1i, x2i, yi), (i = 1, 2, · · · , n) が与えられているとし、
Q =n∑i=1
{yi − (α+ β1x1i + β2x2i)}2 が最小となる α, β1, β2 は、∂Q
∂α=∂Q
∂β1=∂Q
∂β2= 0より正規方程式
1 x1 x2x1 x21 x1x2x2 x1x2 x22
α
β1β2
=
yx1yx2y
を解いたもの
となる。次のように書けることに注意する。(XT は行列 X の転置を表す。)
X =
1 1 · · · 1x11 x12 · · · x1nx21 x22 · · · x2n
とすると XXT = n
1 x1 x2x1 x21 x1x2x2 x1x2 x22
, X
y1y2...yn
= n
yx1yx2y
.
2
また、以下が成り立つ。
回帰式 y = α+ β1x1 + β2x2 は (x1, x2, y)を通る。
残差 ei = yi − yi = yi − (α+ β1x1i + β2x2i)についてn∑i=1
ei = 0,
n∑i=1
x1iei = 0,
n∑i=1
x2iei = 0が成立.
自由度修正決定係数 R2 = 1− 残差変動/(n− k − 1)
総変動/(n− 1).
ここで、残差変動 =
n∑i=1
e2i , 総変動 =
n∑i=1
(yi − y)2, nは観測値の数、k は説明変数の数である。
問題 1.3 五個のデータ (x11, y21, y1), · · · , (x15, y25, y5) が与えられている。ここで、∑x1i = 3,
∑x2i = 2,∑
yi = 5,∑x1ix2i = 4,
∑x1iyi = 12,
∑x2iyi = 8,
∑x21i = 10,
∑x22i = 12,
∑y2i = 16であった。y を
x1, x2 で線形回帰するときの回帰式を求めよ。
• ダミー変数データ (xi, yi), (i = 1, 2, · · · , n) から、奇数時点と偶数時点で定数項 αを変えた回帰式 y = α+ βxを考える。
ダミー変数 di =
{1 i = 奇数
0 i = 偶数
として、データ (xi, di, yi)から、回帰式 y = α+ β1d+ β2xを考えると、
奇数時点では回帰式 y = α+ β1 + β2x, 偶数時点では回帰式 y = α+ β2xと、定数項のみを変えた回帰式を求
めることができる。
同様に、奇数時点と偶数時点で係数 β を変えた回帰式 y = α+ βx は、上記のダミー変数を用い、
データ (xi, dixi, yi)から、回帰式 y = α+ β1x+ β2(dx)を考えると、
奇数時点では回帰式 y = α+ (β1 + β2)x, 偶数時点では回帰式 y = α+ β1xと求めることができる。
問題 1.4 (1) 問題 1.2 (2)のデータについて、定数項ダミーを用いて奇数時点と偶数時点で定数項 α を変え
た回帰式 y = α+ βxを推定せよ。
(2) 問題 1.2 (2) のデータについて、係数ダミーを用いて奇数時点と偶数時点で係数 β を変えた回帰式
y = α+ βxを推定せよ。
(3) 問題 1.2 (3)のデータの前半 4つ上半期、後半 4つは下半期についてであった。適当な定数項ダミー dを
入れることにより、上半期と下半期で定数項 αを変えた回帰式 y = α+ βxを推定せよ。
1.3 非線形回帰
あるタイプの非線形関数で当てはめるべきケースがある。ここでは応用上よく用いられるものを紹介する。
対数線形モデル y = αxβ の両辺の対数をとると、log y = logα+ β log x.
新しい変数として y′ = log y, x′ = log xをとるとよい。
指数関数モデル y = αeβx は、log y = logα+ βxと変形せよ。.
変数 y のみを y′ = log y に変える。
ロジスティック関数モデル y =eα+βx
1 + eα+βx(β > 0) これは微分方程式
dy
dx= βy(1− y) (0 < y < 1) の解
変数 y のみを y′ = logy
1− yに変える。
2項回帰モデル 発生確率 y (0 ≤ y ≤ 1) が説明変数に依存して決まる回帰モデル。
これをある確率分布の分布関数 F を用いて y = F (α+βx)と表すと、y′ = F−1(y)とおくと y′ = α+βx.
F (x) = Φ(x)が N(0, 1)の分布関数のとき、プロビット・モデルという。
F (x) =ex
1 + ex(ロジスティック分布の分布関数) のとき、ロジット・モデルという。これはロジス
ティック関数モデルと同一のものである。
3
問題 1.5 (x, y)のデータが表のとおりに与えられている。このデータから、ロジット・モデル
y =eα+βx
1 + eα+βx(β > 0) を用いた回帰式を求めるとき、α, β の
値を求めよ。また、プロビット・モデルの場合も求めよ。ただし
小数点以下第 2位まで求めよ。
x 1.2 1.4 2.7 3.5 4.8
y 10% 10% 50% 80% 90%
1.4 確率分布の前提を置いた回帰モデルの分析
1.4.1 推定量の分布
説明変数を xi, 誤差項を確率変数 εi とし、被説明変数 Yi を
Yi = α+ βxi + εi, i = 1, 2, · · · , n (1.3)
とし、次を仮定する。
ε1, · · · , εn は独立で εi ∼ N(0, σ2). (1.4)
このとき、最小二乗推定量 α, β および誤差項の分散 σ2 を考える:
β =sxYs2x
=
∑ni=1(xi − x)(Yi − Y )∑n
i=1(xi − x)2, α = Y − βx, σ2 =
1
n− 2
n∑i=1
(Yi − (α+ βxi)
)2. (1.5)
定理 1.1 仮定 (1.4)の下、α, β, σ2 について以下が成立する。
(1)
(α
β
)∼ N
((α
β
), σ2
(1n + x2
ns2x− xns2x
− xns2x
1ns2x
))= N
((α
β
),σ2
n
(1 x
x x2
)−1).
特に、α ∼ N(α, σ2
(1
n+
x2
ns2x
)), β ∼ N
(β,
σ2
ns2x
).
(2)(n− 2)σ2
σ2∼ χ2
n−2. ただし、χ2n−2 は自由度 n− 2のカイ二乗分布を表す。
(3) αと σ2 は独立。また、β と σ2 は独立。
証明: α, β を ε1, · · · , εn の線形結合で表す。ci =xi − x∑n
i=1(xi − x)2とし、
n∑i=1
ci = 0および
n∑i=1
cixi =1
ns2x
( n∑i=1
x2i − xn∑i=1
xi
)=x2 − x2
s2x= 1
に注意すると、
β =
∑ni=1(xi − x)(Yi − Y )∑n
i=1(xi − x)2=
n∑i=1
ci(α+ βxi + εi)−n∑i=1
ciY = β +
n∑i=1
ciεi,
α = Y − βx =1
n
n∑i=1
(α+ βxi + εi)−(β +
n∑i=1
ciεi
)x = α+
n∑i=1
( 1n− cix
)εi.
従って、(1.4)より (α, β)は二次元正規分布に従い、E[α] = α, E[β] = β. また、n∑i=1
c2i =ns2x
(ns2x)2=
1
ns2xに
注意して、
V (β) =
n∑i=1
c2iV (εi) =1
ns2xσ2,
V (α) =
n∑i=1
( 1n− cix
)2V (εi) =
n∑i=1
( 1
n2− 2cix
n+ c2ix
2)σ2 =
( 1n+
x2
ns2x
)σ2,
4
Cov(α, β) =n∑i=1
( 1n− cix
)ciV (εi) = −
x
ns2xσ2
より (1)の ∼は示せる。最後の等号は 1 +x2
s2x=x2
s2xに注意して逆行列を計算せよ。
(2) σ2 は ei = Yi− (α+ βxi)とすると、∑ni=1 ei = 0,
∑ni=1 ciei = 0の二つの制約条件があるため、自由度が
2 つ減って χ2n−2 に従うと説明される。(3) のため厳密な証明の概略を述べる: 1 行目が (1/
√n, · · · , 1/
√n),
2行目が(x1 − x√
ns2x, · · · , xn − x√
ns2x
)で与えらる直行行列を Aとする。このとき、
ε1...
εn
= A
ε1...
εn
と定めるとε1, · · · , εn は独立で εi ∼ N(0, σ2)となる。また、ε1 =
1√n
n∑i=1
εi, ε2 =n∑i=1
xi − x√ns2x
εi,n∑i=1
ε2i =n∑i=1
ε2i より
α− α =1√nε1 −
x√ns2x
ε2, β − β =1√ns2x
ε2,
(n− 2)σ2 =
n∑i=1
{εi − (α− α)− (β − β)xi
}2
=
n∑i=1
{εi −
1√nε1 −
xi − x√ns2x
ε2
}2
= · · · =n∑i=3
ε2i .
よって (2)は明らか。また、α, β は ε1, ε2 の、σ2 は ε3, · · · , εn の関数なので (3)も従う。 □
1.4.2 α, β の区間推定と検定
自由度 nの t分布は独立な Z ∼ N(0, 1)と Y ∼ χ2n を用いて、T =
Z√Y/n
の分布と定義されることに注意
する。また、tn(α)で自由度 nの t分布の上側 α点: T ∼ tn のとき P (T ≥ tn(α)) = αとする。
以下、定理 1.1を引用なしに頻繁に用いる。
• αの信頼区間: Z =α− α√V (α)
∼ N(0, 1), Y =(n− 2)σ2
σ2∼ χ2
n−2 で Z と Y は独立なので、
T =Z√
Y/(n− 2)=
α− α√V (α)σ2/σ2
=α− α√
σ2(1n + x2
ns2x
) ∼ tn−2.従って、信頼度 1− εでの αの信頼区間は
α− tn−2(ε/2)
√σ2( 1n+
x2
ns2x
)≤ α ≤ α+ tn−2(ε/2)
√σ2( 1n+
x2
ns2x
).
• β の信頼区間: Z =β − β√V (β)
∼ N(0, 1), Y =(n− 2)σ2
σ2∼ χ2
n−2 で Z と Y は独立なので、
T =Z√
Y/(n− 2)=
β − β√V (β)σ2/σ2
=β − β√σ2
1
ns2x
∼ tn−2.
従って、信頼度 1− εでの β の信頼区間は
β − tn−2(ε/2)
√σ2
ns2x≤ β ≤ β + tn−2(ε/2)
√σ2
ns2x.
注意 1.1 σ2 を計算する際は σ2 =1
n− 2(1− rxy2)ns2y と計算するのがよい。これは次のように導かれる。
σ2 =1
n− 2
n∑i=1
e2i =1
n− 2(残差変動) =
1
n− 2(全変動−回帰変動) =
1
n− 2(1−決定係数 R2)(全変動)
5
=1
n− 2(1− rxy2)ns2y.
• 検定: 次の手順で有意水準 εの両側検定を行うことができる。
帰無仮説 H0 : β = β0, 対立仮説 H1 : β = β0
H0 のもとで、T =β − β√σ2
ns2x
∼ tn−2. よって、t分布表から tn−2(ε/2)を求め、標本からの実現値 tに対して、
|t| > tn−2(ε/2)なら H0 を棄却、|t| ≤ tn−2(ε/2)なら H0 を採択 すればよい。
同様に、片側検定の場合、
帰無仮説 H0 : β = β0, 対立仮説 H1 : β > β0 のときは、
t > tn−2(ε)なら H0 を棄却、t ≤ tn−2(ε)なら H0 を採択 すればよい。
帰無仮説 H0 : β = β0, 対立仮説 H1 : β < β0 のときは、
t < tn−2(ε)なら H0 を棄却、t ≥ tn−2(ε)なら H0 を採択 すればよい。
問題 1.6 問題 1.2 (2)のデータについて、σ2 の実現値を求め、α, β の 95%信頼区間を求めよ。ただし小数点
以下第 3位まで求めよ。
問題 1.7 問題 1.2 (3)のデータについて、帰無仮説 H0 : β = 0, 対立仮説 H1 : β > 0を、有意水準 5%で検
定せよ。
1.4.3 点予測、区間予測
説明変数 xn+1 が与えられたときの Yn+1 の予測量 Yn+1 は、α, β を用いて、Yn+1 = α + βxn+1 となり。
これは正規分布に従う。(これは α, β が ε1, · · · , εn の線形結合であることによる。)
予測誤差 Yn+1 − Yn+1 = −(α− α)− (β − β)xn+1 + εn+1 について定理 1.1(1)より、
E[Yn+1 − Yn+1] = −E[α− α]− xn+1E[β − β] + E[εn+1] = 0,
V (Yn+1 − Yn+1) = V (εn+1 − (α− α)− (β − β)xn+1) = V (εn+1 − α− βxn+1)
= V (εn+1) + V (−α− βxn+1), (∵ εn+1 は α, β と独立)
= σ2 + V (α) + 2xn+1 Cov(α, β) + x2n+1V (β)
= σ2 + σ2( 1n+
x2
ns2x
)+ 2xn+1
−xσ2
ns2x+ x2n+1
σ2
ns2x
= σ2
{1 +
1
n+
(xn+1 − x)2
ns2x
}従って、予測誤差 Yn+1 − Yn+1 ∼ N
(0, σ2
{1 +
1
n+
(xn+1 − x)2
ns2x
})となる。
これより、σ2 が既知であればこれより区間推定できる。
σ2 が未知の場合、Z =Yn+1 − Yn+1√V (Yn+1 − Yn+1)
∼ N(0, 1), W =(n− 2)σ2
σ2∼ χ2
n−2 で Z とW は独立なので、
T =Z√
W/(n− 2)=
Yn+1 − Yn+1√V (Yn+1 − Yn+1)σ2/σ2
=Yn+1 − Yn+1√
σ2
{1 +
1
n+
(xn+1 − x)2
ns2x
} ∼ tn−2.従って、信頼度 1− εでの Yn+1 の信頼区間は
Yn+1− tn−2(ε2
)√σ2
{1 +
1
n+
(xn+1 − x)2
ns2x
}≤ Yn+1 ≤ Yn+1 + tn−2
(ε2
)√σ2
{1 +
1
n+
(xn+1 − x)2
ns2x
}となる。ただし、上式で Yn+1 と σ2 は実現値を表す。
6
問題 1.8 問題 1.2 (3)のデータに対して推定された回帰式を用いて (問題 1.7も参照のこと)、x8+1 = 4に対
する点予測および信頼係数 95%信頼区間を求めよ。ただし小数点以下第 3位まで求めよ。
1.5 統計の復習 1 正規母集団と二項母集団
定義 1.1 正規母集団の統計において次の分布は特に重要である。
χ2 分布: X1, . . . , Xn が i.i.d.で N(0, 1)に従うとき、X21 + · · ·+X2
n ∼ χ2n (自由度 nの χ2 分布).
t分布: Y, Z は独立で Y ∼ χ2n, Z ∼ N(0, 1)のとき、T =
Z√Y/n
∼ tn (自由度 nの t分布).
F 分布: X,Y は独立で X ∼ χ2m, Y ∼ χ2
n のとき、W =X/m
Y/n∼ Fmn (自由度 (m,n)の F 分布).
次の定理は確率統計学 Iで定理 3.7で示した。
定理 1.2 X1, . . . , Xn が独立で、それぞれ同一の正規分布 N(µ, σ2)に従うとするとき、次が成立する。
(1) 標本平均 X =1
n
n∑i=1
Xi は N(µ,σ2
n
)に従う。
(2) 不偏分散 U2 =1
n− 1
n∑i=1
(Xi −X)2 について、n− 1
σ2U2 =
1
σ2
n∑i=1
(Xi −X)2 ∼ χ2n−1.
(3) X と U2 は独立。
応用例
• 正規母集団の母平均の区間推定、検定において、母分散 σ2 が既知の場合、定理 1.2 (1)を用いて行うことが
できた。例えば、標本平均の実現値が xのとき、信頼度 1− εでの母平均 µのの信頼区間は
x− u(ε/2)√σ2
n≤ µ ≤ x+ u(ε/2)
√σ2
n,
ここで、u(α)は N(0, 1)の上側 α点を表す。母分散が未知であっても、標本数が大きい場合は母分散を不偏
分散の実現値としてこれを用いた。
• 正規母集団の母分散の区間推定、検定において、定理 1.2 (2)を用いて行うことができた。例えば、不偏分散
の実現値が u2 のとき、信頼度 1− εでの母分散 σ2 の信頼区間は
(n− 1)u2
χ2n−1(ε/2)
≤ σ2 ≤ (n− 1)u2
χ2n−1(1− ε/2)
,
ここで、χ2n−1(α)は χ2
n−1 の上側 α点を表す。
•正規母集団の母平均の区間推定、検定において、母分散 σ2が未知の場合、定理 1.2より T =X − µ√U2/n
∼ tn−1
となること (各自証明を試みよ)を用いて行うことができた。例えば、標本平均, 不偏分散の実現値が x, u2 の
とき、信頼度 1− εでの母平均 µのの信頼区間は
x− tn−1(ε/2)√u2
n≤ µ ≤ x+ tn−1(ε/2)
√u2
n,
ここで、tn−1(α)は tn−1 の上側 α点を表す。
• 2つの正規母集団の母数の比較に関する検定を前期の数理統計学 Iの最後の節で取り上げた。これは区間推
定にも用いられる。例えば、
X1, . . . , Xm, Y1, . . . , Ynをそれぞれ正規母集団N(µ1, σ12), N(µ2, σ2
2)からの無作為標本とする。このとき、
標本平均をX, Y ,標本分散を S2X , S2
Y とすると (母分散 σ12と σ2
2は既知であれば)、Z =X − Y − (µ1 − µ2)√
σ12/m+ σ22/nは N(0, 1)に従う。従って、標本平均の実現値を x, yとすると、平均の差 µ1−µ2 の信頼度 1− εの信頼区間は
x− y − u(ε/2)√σ21
m+σ22
n≤ µ1 − µ2 ≤ x− y + u(ε/2)
√σ21
m+σ22
n,
7
ここで、母分散が未知であっても、標本数が大きい場合は母分散 σ21 , σ
22 をそれぞれの不偏分散の実現値 u21, u
22
で置き換えて成立する。
• 上記は大標本での二項母集団の区間推定や検定にも用いることができる。例えば、母比率 p1 の二項母集団からの大きさ mの標本比率を P1, 母比率 p2 の二項母集団からの大きさ n
の標本比率を P2 とすると、二項分布の正規分布近似を考え、
P1 ∼ N(p1,
p1(1− p1)m
), P2 ∼ N
(p2,
p2(1− p2)n
)よりP1 − P2 ∼ N
(p1 − p2,
p1(1− p1)m
+p2(1− p2)
n
).
これより、標準化し、標本比率と根号内の母比率をその実現値 p1, p2 に置き換えることで、母比率の差 p1− p2の信頼度 1− εの信頼区間
p1 − p2 − u(ε/2)√p1(1− p2)
m+p2(1− p2)
n≤ p1 − p2 ≤ p1 − p2 + u(ε/2)
√p1(1− p2)
m+p2(1− p2)
n
を得る。
例題 1.1 ある政策の支持率を予想するために、母集団から男性 900人、女性 800人をそれぞれ無作為に抽出
して調査を行ったところ、男性は 300人、女性は 320人が支持すると回答した。母集団全体の男女比は 5 : 4
であるとして、母集団全体での支持率を近似法を用いて、信頼度 95%で区間推定せよ。
解: 男女の支持率を p1, p2 とし、標本比率を P1, P2 とする。このとき、P1 ∼ N(p1,
p1(1− p1)900
), P2 ∼
N(p2,
p2(1− p2)800
)と近似される。男女比を考慮すると全体の支持率は P =
5
9P1 +
4
9P2 となるから、P1 と
P2 は独立なので、P ∼ N(59p1 +
4
9p2,(59
)2 p1(1− p1)900
+(49
)2 p2(1− p2)800
). これより、標準化し、標本比
率と根号内の母比率をその実現値 p1 =300
900, p2 =
320
800に置き換えることで、
5
9p1 +
4
9p2 ± u(0.025)
√(59
)2 p1(1− p1)900
+(49
)2 p2(1− p2)800
= 0.36296 · · · ± 0.02281 · · · ={0.3857 · · ·0.3401 · · ·
従って、0.340 ≤ p ≤ 0.386. □
問題 1.9 ある都市の市長選挙の結果を予想するために、60才未満の者 120人に意見を求めたところ 48人が
保守系を支持すると言った。一方、60才以上の者 80人について調べたところ、56人が保守系を支持した。
(1) 近似法を用いて、60 才以上の人の支持率と 60 才未満の人の支持率の差の信頼係数 95% 信頼区間を求め
よ。ただし小数点以下第 3位まで求めよ。
(2) 投票率を考慮すると、この都市の 60才未満の人と 60才以上の人の比は 4 : 5である。近似法を用いて、こ
の選挙での保守系の得票率を信頼度 95%で区間推定せよ。ただし小数点以下第 3位まで求めよ。
2 時系列解析
2.1 時系列に現れる確率過程と用語の定義
確率過程 Yt, t = 0,±1,±2, . . .を考える。ここでは、時間パラメータ tは負の値もとることに注意する。
定義 2.1 確率過程 Yt, t = 0,±1,±2, . . .が定常 (stationary)であるとは、
条件 1 E[Yt] = µ (定数)
条件 2 Cov(Yt, Yt−h) = γh, 特に、Cov(Yt, Yt) = V (Yt) = γ0 (定数)
となるときにいう。
8
定義 2.2 γh = Cov(Yt, Yt−h)を時差 hの自己共分散という。
ρh =γhγ0
=Cov(Yt, Yt−h)
V (Yt)を時差 hの自己相関 (コレログラム)という。
定常過程の典型例: Xt, t = 0,±1,±2, . . ., を i.i.d. とし、数列 {an}が∑∞n=0 |an| <∞を満たすとする。
Yt =∞∑i=0
aiXt−i とおくと、これは定常過程。実際、E[Xt] = µ, V (Xt) = σ2 とすると、
E[Yt] =
∞∑i=0
aiE[Xt−i] = µ
∞∑i=0
ai,
Cov(Yt, Yt−h) = Cov( ∞∑i=0
aiXt−i,
∞∑j=0
ajXt−h−j
)=
∞∑i=0
∞∑j=0
aiaj Cov(Xt−i, Xt−h−j)
= σ2∞∑i=0
ah+jaj , (∵ Cov(Xt−i, Xt−h−j) = 0 ⇐⇒ i = h+ j)
定義 2.3 {yt}nt=1 を {Yt}nt=1 の実現値とする。このとき y = 1n
∑ni=1 yt とする。
γh =1
n
∑nt=h+1(yt − y)(yt−h − y)を時差 hの標本自己共分散という。
ρh =γhγ0
=
∑nt=h+1(yt − y)(yt−h − y)∑n
t=1(yt − y)2を時差 hの標本自己相関という。
2.2 AR(p) (p次の自己回帰モデル, Auto-regressive Model)
{Yt}が p次の自己回帰モデル AR(p)に従うとは、
Yt = ϕ0 + ϕ1Yt−1 + ϕ2Yt−2 + · · ·+ ϕpYt−p + εt (2.1)
となることをいう。ここで、ϕ0, ϕ1, . . . , ϕp は定数であり、誤差項 εt は Yt−1, Yt−2, . . ., εt−1, εt−2, . . .と独立
で E[εt] = 0, V (εt) = σ2 とし、同一の分布に従う。{Yt} ∼ AR(p)と表す。
命題 2.1 (定常性) 特性方程式 ϕ(x) = 1− (ϕ1x+ ϕ2x2 + . . .+ ϕpx
p) = 0の解*1の絶対値がすべて 1より大
きいとき、AR(p)は定常性をもつ。
p = 1, つまり、Yt = ϕ0+ϕ1Yt−1+εtのとき、ϕ(x) = 1−ϕ1xより、定常 ⇐⇒ |1/ϕ1| > 1 ⇐⇒ |ϕ1| < 1.
このとき、(2.1)の tを t− 1, t− 2, . . . として代入していくと、
Yt = ϕ0 + ϕ1Yt−1 + εt = ϕ0 + ϕ1(ϕ0 + ϕ1Yt−2 + εt−1) + εt
= ϕ0(1 + ϕ1) + ϕ21Yt−2 + (εt + ϕ1εt−1) = ϕ0(1 + ϕ1) + ϕ21(ϕ0 + ϕ1Yt−3 + εt−3) + (εt + ϕ1εt−1)
= ϕ0(1 + ϕ1 + ϕ21) + ϕ31Yt−3 + (εt + ϕ1εt−1 + ϕ21εt−2)
= · · ·= ϕ0(1 + ϕ1 + · · ·+ ϕn1 ) + ϕn+1
1 Yt−n−1 + (εt + ϕ1εt−1 + · · ·+ ϕn1 εt−n).
よって、|ϕ1| < 1に注意して n→∞とすると次を得る。
Yt =ϕ0
1− ϕ1+
∞∑i=0
ϕi1εt−i. これは AR(1)のMA(∞)表現である。 (2.2)
{Yt} ∼ AR(p)とし、(2.1)で期待値をとると、µ = E[Yt] = ϕ0 + ϕ1µ+ ϕ2µ+ · · ·+ ϕpµとなるので
Yt − µ = ϕ1(Yt−1 − µ) + · · ·+ ϕp(Yt−p − µ) + εt.
*1 ラグ作用素 L: LYt = Yt−1 を用いると、(2.1) は Yt = (ϕ1L+ ϕ2L2 + · · ·+ ϕpLp)Yt + ϕ0 + εt, 即ち ϕ(L)Yt = ϕ0 + εt と表せる。
9
これに、Yt − µ, Yt−1 − µ, · · · をかけて期待値をとれば、
γ0 = ϕ1γ1 + ϕ2γ2 + · · ·+ ϕpγp + σ2
γ1 = ϕ1γ0 + ϕ2γ1 + · · ·+ ϕpγp−1
γ2 = ϕ1γ1 + ϕ2γ0 + · · ·+ ϕpγp−2...
γp = ϕ1γp−1 + ϕ2γp−2 + · · ·+ ϕpγ0
γp+1 = ϕ1γp + ϕ2γp−1 + · · ·+ ϕpγ1
γp+t = ϕ1γp+t−1 + ϕ2γp+t−2 + · · ·+ ϕp−1γt+1 + ϕpγt (t ≥ 0)
が成立する。逆に上式の両辺を γ0 で割ることで、µ, γ0, γ1, . . . , γp は、ϕ0 = (1− (ϕ1 + · · ·+ ϕp))µおよびρ1 = ϕ1 + ϕ2ρ1 + · · ·+ ϕpρp−1ρ2 = ϕ1ρ1 + ϕ2 + · · ·+ ϕpρp−2
...ρp = ϕ1ρp−1 + ϕ2ρp−2 + · · ·+ ϕp
行列で
ρ1ρ2...ρp
=
1 ρ1 · · · ρp−1ρ1 1 · · · ρp−2...
.... . .
...ρp−1 ρp−2 · · · 1
ϕ1ϕ2...ϕp
(2.3)
を満たす。これをユール=ウォーカーの方程式という。ρh = γh/γ0 に注意。
• AR(1): Yt = ϕ0 + ϕ1Yt−1 + εt のとき、 µ = ϕ0 + ϕ1µより µ =ϕ0
1− ϕ1.
γ0 = ϕ1γ1 + σ2, γ1 = ϕ1γ0, γ2 = ϕ1γ1, · · · より
γ0 =σ2
1− ϕ21, γ1 =
ϕ1σ2
1− ϕ21. また、γh =
ϕh1σ2
1− ϕ21, 自己相関 ρh = ϕh1 .
• AR(2): Yt = ϕ0 + ϕ1Yt−1 + ϕ2Yt−2 + εt のとき、
特性方程式 ϕ(x) = 1− ϕ1x− ϕ2x2 = 0のすべての解の絶対値が 1より大きい条件を求めると、
(i) 実数解をもつ場合 D = ϕ21 + 4ϕ2 ≥ 0, ϕ(−1) = 1 + ϕ1 − ϕ2 > 0, ϕ(1) = 1− ϕ1 − ϕ2 > 0.
(ii) 虚数解をもつ場合 D = ϕ21 + 4ϕ2 < 0,∣∣∣∣∣−ϕ1 ±√−(ϕ21 + 4ϕ2)i
2ϕ2
∣∣∣∣∣2
=
(−ϕ12ϕ2
)2
+
(√−(ϕ21 + 4ϕ2)
2ϕ2
)2
=ϕ21 + {−(ϕ21 + 4ϕ2)}
4ϕ22=−1ϕ2
> 1
ϕ2 < 0より、ϕ2 > −1となる。以上をまとめると、
AR(2)が定常性をもつ⇐⇒ ϕ2 = 1− |ϕ1|かつ ϕ2 > −1⇐⇒ (ϕ1, ϕ2)が (0, 1), (−2,−1), (2,−1)を頂点とする三角形の内部
特に、放物線 ϕ2 = − 14ϕ
21 の上側が実数解、下側が虚数解に対応する。(ϕ1ϕ2-平面に図示せよ。)
平均: 期待値をとって、µ = ϕ0 + ϕ1µ+ ϕ2µ+ 0 より µ =ϕ0
1− (ϕ1 + ϕ2).
分散、自己共分散: Yt − µ = ϕ1(Yt−1 − µ) + · · ·+ ϕ2(Yt−2 − µ) + εt より、
γ0 = ϕ1γ1 + ϕ2γ2 + σ2,
γ1 = ϕ1γ0 + ϕ2γ1, γ1 =ϕ1
1− ϕ2γ0
γ2 = ϕ1γ1 + ϕ2γ0, γ2 =( ϕ211− ϕ2
+ ϕ2
)γ0
γn = ϕ1γn−1 + ϕ2γn−2 (n ≥ 2) この漸化式を解けば γn が計算できる。
これより γ0 =ϕ21
1− ϕ2γ0 +
( ϕ211− ϕ2
+ ϕ2
)γ0 + σ2, すなわち、
γ0 =σ2
1− ϕ21
1−ϕ2− (
ϕ21
1−ϕ2+ ϕ2)
=1− ϕ21 + ϕ2
σ2
(1− ϕ2)2 − ϕ21> 0.
10
問題 2.1 定常 AR(2)モデル Yt = 3 − 0.3Yt−1 + 0.4Yt−2 + εt, E[εt] = 0, V (εt) = 3において、µ = E[Yt],
自己共分散 γ0, γ1, γ2, γn (n ≥ 3), 自己相関 ρ1, ρ2, ρn (n ≥ 3) を求めよ。
問題 2.2 定常 AR(2)モデル Yt = 1 +2
3Yt−1 +
1
6Yt−2 + εt, で、
(1) 定常性を示せ。(2) µ = E[Yt]を求めよ。(3) σ2 = 2として γ0, γ1, γ2, γ3 を求めよ。
問題 2.3 (1) 定常 AR(2)モデル Yt = 1 + aYt−1 + 0.2Yt−2 + εt で、定常性をもつための aの範囲を求めよ。
(2) 定常 AR(2)モデル Yt = 1− 0.5Yt−1 + bYt−2 + εt で、定常性をもつための bの範囲を求めよ。
2.3 MA(q) (q次の移動平均モデル, Moving-average Model)
{Yt}が q 次の移動平均モデルMA(q)に従うとは、
Yt = θ0 + εt − θ1εt−1 − θ2εt−2 − · · · − θqεt−q (2.4)
となることにいう。ここで、θ0, θ1, . . . , θq は定数であり、誤差項 εt は Yt−1, Yt−2, . . ., εt−1, εt−2, . . .と独立で
E[εt] = 0, V (εt) = σ2 とし、同一の分布に従う。{Yt} ∼ MA(q)と表す。
MA(q)は必ず定常性を満たす。実際、E[Yt] = θ0 であり、
γ0 = V (Yt) = (1 + θ21 + · · ·+ θ2q)σ2,
γ1 = Cov(Yt, Yt−1) = Cov(εt − θ1εt−1 − θ2εt−2 − · · · − θqεt−q, εt−1 − θ1εt−2 − θ2εt−3 − · · · − θqεt−q−1)= (−θ1 + θ1θ2 + θ2θ3 + · · ·+ θq−1θq)σ
2.
同様に γ2 = (−θ2 + θ1θ3 + · · ·+ θq−2θq)σ2, · · · , γq = −θqσ2, γn = 0 (n ≥ q + 1).
自己相関 ρh は ρ0 = 1, ρ1 =−θ1 + θ1θ2 + θ2θ3 + · · ·+ θq−1θq
1 + θ21 + · · ·+ θ2q, ρ2 =
−θ2 + θ1θ3 + · · ·+ θq−2θq1 + θ21 + · · ·+ θ2q
, · · · .
定常な AR(p) {Yt}は、Yt = ξ0 + εt +∑∞i=1 ξiεt−i と表されることより、MA(∞).
実際、AR(p)のMA(∞)表現を、Yt = ξ0 + εt + ξ1εt−1 + ξ2εt−2 + · · · とすると、
Yt = ϕ0 + ϕ1Yt−1 + · · ·ϕpYt−p + εt
= ϕ0 + ϕ1(ξ0 + εt−1 + ξ1εt−2 + ξ2εt−3 + · · · ) + ϕ2(ξ0 + εt−2 + ξ1εt−3 + ξ2εt−4 + · · · ) ++ · · ·+ ϕp(ξ0 + εt−p + ξ1εt−p−1 + ξ2εt−p−2 + · · · ) + εt
= ϕ0 + (ϕ1 + · · ·+ ϕp)ξ0 + εt + ϕ1εt−1 + (ϕ1ξ1 + ϕ2)εt−2 + (ϕ1ξ2 + ϕ2ξ1 + ϕ3)εt−3 + · · · .
よって、ξ0 =ϕ0
1− (ϕ1 + · · ·+ ϕp)(= µ), ξ1 = ϕ1, ξ2 = ϕ1ξ1 + ϕ2 = ϕ21 + ϕ2, ξ3 = ϕ1ξ2 + ϕ2ξ1 + ϕ3 = · · · .
逆に、上で Yt と εt の役割を入れ替えれば、次がわかる。
命題 2.2 (反転可能性) 特性方程式 θ(x) = 1 − (θ1x + θ2x2 + . . . + θqx
q) = 0の解の絶対値がすべて 1より
大きいとき、MA(q)は AR(∞)で表現できる。すなわち、反転可能である。
命題 2.3 (識別可能性) 特性方程式 θ(x) = 0の解の絶対値がすべて 1以上のとき、MA(q)は識別可能である。
すなわち、平均、分散、自己共分散からモデルが一意的に定まる。
例題 2.1 Yt がMA(1), つまり Yt = θ0 + εt − θ1εt−1 について、µ = E[Yt], γ0 = V (Yt), γ1 = Cov(Yt, Yt−1)
のとき、θ0, θ1, σ2 はどうなるか。
11
解: 別の δ0, δ1, ω2 があったとすると、
µ = θ0 = δ0, γ0 = σ2(1 + θ21) = ω2(1 + δ21), γ1 = (−θ1)σ2 = (−δ1)ω2.
γ0, γ1 の式より
θ11 + θ21
=δ1
1 + δ21, ゆえに θ1(1 + δ21)− δ1(1 + θ21) = (θ1 − δ1)(1− θ1δ1) = 0
よって、θ1δ1 = 1であれば、別のパラメータでも同じ µ, γ1, γ1 をとり得る。
ここで、特性方程式 1− θ1x = 0の解は x =1
θ1. よって、
∣∣∣ 1θ1
∣∣∣ ≥ 1, すなわち |θ1| ≤ 1となる方を取ってく
れば、一意的に定まる。 □
問題 2.4 MA(2)モデル Yt = 3+ εt − aεt−1 + aεt−2 が反転可能性をもつ実数 aの範囲を求めよ。また、識別
可能性をもつ実数 aの範囲を求めよ。
問題 2.5 MA(2)モデル Yt = 2+ εt +5
6εt−1 −
1
6εt−2, E[εt] = 0, V (εt) = 2において、µ = E[Yt], 自己共分
散 γ0, γ1, γ2, γ3 を求めよ。
問題 2.6 MA(2)モデル Yt = 1 + εt −1
4εt−1 −
1
4εt−2, ただし、εt ∼ N(0, 32)とする。このとき、自己共分
散 γ0, γ1, γ2, γ3 を求めよ。また、Yt の分布、(Yt, Yt+1)の同時分布を求めよ。
問題 2.7 定常AR(2)モデル Yt = 1+0.8Yt−1+0.3Yt−2+εtのMA(∞)表現 Yt = ξ0+εt+ξ1εt−1+ξ2εt−2+· · ·における、ξ0, ξ1, ξ2, ξ3 を求めよ。
問題 2.8 定常 AR(2)モデル Yt = 3 +5
6Yt−1 −
1
6Yt−2 + εt において、次の二つの方法で Yt のMA(∞)表現
を求めよ。
(1) Yt = µ+∞∑i=0
aiεt−i とおくと、{ai}が ai+2 =5
6ai+1 −
1
6ai, a0 = 1, a1 =
5
6を満たすことを示す。
(2) ラグ作用素 LYt = Yt−1 を用い、ϕ(L)(Yt − µ) = εt とし、1/ϕ(x)を xのべき級数で表す。
2.4 ARMA(p, q) (Autoregressive Moving-average Model)
{Yt}が次数 p, q の自己回帰移動平均モデル ARMA(q, p)に従うとは、
Yt = c+ ϕ1Yt−1 + ϕ2Yt−2 + · · ·+ ϕpYt−p + εt − θ1εt−1 − θ2εt−2 − · · · − θqεt−q (2.5)
となることにいう。ここで、c, ϕ1, . . . , ϕp, θ1, . . . , θq は定数であり、誤差項 εt は Yt−1, Yt−2, . . ., εt−1, εt−2, . . .
と独立で E[εt] = 0, V (εt) = σ2 とし、同一の分布に従う。{Yt} ∼ ARMA(q)と表す。
命題 2.4 (1) ϕ(x) = 1−(ϕ1x+ϕ2x2+. . .+ϕpxp) = 0の解の絶対値がすべて 1より大きいとき、ARMA(p, q)
は定常性をもつ。
(2) θ(x) = 1− (θ1x+ θ2x2 + . . .+ θqx
q) = 0の解の絶対値がすべて 1より大きいとき、ARMA(p, q)は反転
可能である。
(3) θ(x) = 0の解の絶対値がすべて 1以上であり、かつ ϕ(x) = 0と共通解をもたないとき、ARMA(p, q)は
識別可能である。
例題 2.2 定常な ARMA(1, 1)モデル Yt = ϕ1Yt−1 + εt − θ1εt−1 のパラメータ ϕ1, θ1 および σ2(= V (εt))を
用いて自己相関 ρh (h = 1, 2, · · · ) を表せ。
12
解: Yt − ϕ1Yt−1 = εt − θ1εt−1 で
V (Yt − ϕ1Yt−1) = V (Yt)− 2ϕ1 Cov(Yt, Yt−1) + ϕ21V (Yt−1) = γ0(1 + ϕ21)− 2ϕ1γ1,
V (εt − θ1εt−1) = σ2(1 + θ21),
より γ0(1 + ϕ21)− 2ϕ1γ1 = σ2(1 + θ21). また、
γ1 = Cov(Yt, Yt+1) = Cov(Yt, ϕ1Yt + εt+1 − θ1εt)= ϕ1 Cov(Yt, Yt) + Cov(ϕ1Yt−1 + εt − θ1εt−1, εt+1 − θ1εt) = ϕ1γ0 − θ1σ2.
よって、γ0(1 + ϕ21)− 2ϕ1(ϕ1γ0 − θ1σ2) = σ2(1 + θ21)より γ0 =σ2(1 + θ21 − 2ϕ1θ1)
1− ϕ21,
γ1 =σ2(ϕ1 + ϕ1θ
21 − ϕ21θ1 − θ1)
1− ϕ21=σ2(θ1 − ϕ1)(ϕ1θ1 − 1)
1− ϕ21. ∴ ρ1 =
γ1γ0
=(θ1 − ϕ1)(ϕ1θ1 − 1)
1 + θ21 − 2ϕ1θ1.
また、h ≥ 2のとき、
γh = Cov(Yt, Yt+h) = Cov(Yt, ϕ1Yt+h−1 + εt+h − θ1εt+h−1) = ϕ1γh−1.
よって、γh = ϕh−11 γ1 より、ρh = ϕh−11 ρ1 =(θ1 − ϕ1)(ϕ1θ1 − 1)
1 + θ21 − 2ϕ1θ1(h ≥ 2). □
問題 2.9 定常ARMA(1, 1)モデル Yt = 2+1
3Yt−1−
1
2εt−1+εt, E[εt] = 0, V (εt) = 42において、µ = E[Yt],
γh = Cov(Yt, Yt−h)を求めよ。
2.5 時系列モデルに基づく予測
• 偏自己相関: 確率過程の自己相関を ρ1, ρ2, ρ3, · · · とし、Yt = ϕ0 + ϕ1Yt−1 + ϕ2Yt−2 + εt のとき、
ρ1 = ϕ11 ρ1 = ϕ21 + ϕ22ρ1 ρ1 = ϕ31 + ϕ32ρ1 + ϕ33ρ2
ρ2 = ϕ21ρ1 + ϕ22 ρ2 = ϕ31ρ1 + ϕ32 + ϕ33ρ1 · · ·ρ3 = ϕ31ρ2 + ϕ32ρ1 + ϕ33
で定まる ϕhh を Yt と Yt−h の偏自己相関という。(ユール=ウォーカーの方程式 (2.3)で ϕi を ϕhi とした。)
問題 2.10 次の定常 AR(2)モデル, MA(2)モデルにおいて偏自己相関 ϕ11, ϕ22, ϕ33 を求めよ。
(1) 定常 AR(2)モデル Yt = 3− 0.3Yt−1 + 0.4Yt−2 + εt, (cf . 問題 2.1).
(2) 定常 AR(2)モデル Yt = 1 +2
3Yt−1 +
1
6Yt−2 + εt, (cf . 問題 2.2).
(3) 定常MA(2)モデル Yt = 2 + εt +5
6εt−1 −
1
6εt−2, (cf . 問題 2.5).
• 時点 nまでのデータ {yt}が与えられたときの h時点先の Yn+h の予測:
一般に、Y の X による予測値として E[Y |X]をとる。これを (2乗平均)最良予測値という。
これは、E[(Y − g(X))2] を最小にするのは g(X) = E[Y |X] となることによる。また、
E[f(X)|X] = f(X), X と Y が独立なら E[X|Y ] = E[X]
に注意する。(前期に定理 1.7として示した。)
[AR(p)の予測量] Yt = ϕ0 + ϕ1Yt−1 + ϕ2Yt−2 + · · ·+ ϕpYt−p + εt のとき、
Yn+1 の予測値として次の Yn+1 をとる。
Yn+1 = E[Yn+1|Yn, Yn−1, . . .] = E[ϕ0 + ϕ1Yn + ϕ2Yn−1 + · · ·+ ϕpYn−p+1 + εn+1|Yn, Yn−1, . . .]= ϕ0 + ϕ1Yn + ϕ2Yn−1 + · · ·+ ϕpYn−p+1
13
最後の等号は εn+1 が Yn, Yn−1, . . .と独立であることと E[εn+1] = 0による。
Yn+1 = Yn+1 + εn+1 に注意して
Yn+2 = ϕ0 + ϕ1Yn+1 + ϕ2Yn + · · ·+ ϕpYn−p+2 + εn+2
= ϕ0 + ϕ1(Yn+1 + εn+1) + ϕ2Yn−1 + · · ·+ ϕpYn−p+1 + εn+1,
Yn+2 = E[Yn+2|Yn, Yn−1, . . .] = ϕ0 + ϕ1Yn+1 + ϕ2Yn + · · ·+ ϕpYn−p+2
同様に、Yn+3 = ϕ0 + ϕ1Yn+2 + ϕ2Yn+1 + ϕ3Yn + · · ·+ ϕpYn−p+3 となる。
[MA(q)の予測量] Yt = θ0 + εt − θ1εt−1 − θ2εt−2 − · · · − θqεt−q のとき、
Yn+1 = E[Yn+1|εn, εn−1, . . .] = θ0 − θ1εn − · · · − θq−1εn+2−q − θqεn+1−q,
Yn+2 = E[Yn+2|εn, εn−1, . . .] = θ0 − θ2εn − · · · − θqεn+2−q,
...
Yn+q = E[Yn+2|εn, εn−1, . . .] = θ0 − θqεn,Yn+h = 0, h = q + 1, q + 2, . . . .
すなわち、εn+1, . . . , εn+h をゼロとしたものになる。
なお、誤差項 εn, εn−1, . . .は直接観察されるものではないため、パラメータ推定等の結果として得られる残
差 en, en−1, . . .を用いて計算するものとなる。(残差についてはモデリングの教科書 p.2-22を参照のこと。)
問題 2.11 定常 AR(1)モデル Yt = 2 +1
3Yt−1 + εt, εt ∼ N
(0,
1
5
), について、YT = 2.7が与えられたとす
る。
(1) 最良予測値 YT+1 を求め、2乗平均誤差 E[(YT+1 − YT+1)2|YT = 2.7]を求めよ。
また、最良予測値 YT+2 を求め、2乗平均誤差 E[(YT+2 − YT+2)2|YT = 2.7]を求めよ。
(2) YT+1 と YT+2 の信頼度 95%の信頼区間を求めよ。ただし小数点以下第 2位まで求めよ。
問題 2.12 定常MA(2)モデル Yt = 1 + εt −2
3εt−1 −
1
3εt−2, εt ∼ N
(0,
1
4
), について、
εT−2 = −1
3, εT−1 =
1
3, εT =
2
3が与えられたとする。
(1) 最良予測値 YT+1 を求め、2乗平均誤差 E[(YT+1 − YT+1)
2∣∣∣εT−2 = −1
3, εT−1 =
1
3, εT =
2
3
]を求めよ。
また、最良予測値 YT+2 を求め、2乗平均誤差 E[(YT+2− YT+2)
2∣∣∣εT−2 = −1
3, εT−1 =
1
3, εT =
2
3
]を求めよ。
(2) YT+1 と YT+2 の信頼度 95%の信頼区間を求めよ。ただし小数点以下第 3位まで求めよ。
2.6 統計の復習 2 順序統計量
前期に扱った 1.7 順序統計量 を復習する。
X1, . . . , Xn が i.i.d.で連続型とする。X1, . . . , Xn を小さいものから並べたものを
X(1), X(2), . . . , X(n)
とし、X(j) を j 番目の順序統計量という。X1 の (共通の)密度関数を f(x), 分布関数を F (x)とする。
明らかに、X(1) = min{X1, . . . , Xn}, X(n) = max{X1, . . . , Xn}である。
また、n = 2m− 1(奇数)のとき X(m), n = 2m(偶数)のとき1
2{X(m) +X(m+1)}が X1, . . . , Xn の中央値で
あった。
対称性に注意すると、(X(1), . . . , X(n))の同時密度関数 g が
g(t1, t2, · · · , tn) ={n!f(t1)f(t2) · · · f(tn) (t1 < t2 < · · · < tn)
0 (その他)
14
となることがわかる。また、同様に対称性より x < y のとき、∫∫· · ·∫x<t1<t2<···<tk<y
f(t1)f(t2) · · · f(tk) dt1dt2 · · · dtk =1
k!{F (y)− F (x)}k
ここで、F (x)は X1 の分布関数である。よって、X(j) の密度関数は、この周辺密度なので
fX(j)(x) = n!
∫· · ·∫
t1<···<tj−1<x
f(t1) · · · f(tj−1) dt1 · · · dtj−1f(x)∫· · ·∫
x<tj+1<···<tn
f(tj+1) · · · f(tn) dtj+1 · · · dtn
=n!
(j − 1)!1!(n− j)!F (x)j−1f(x)(1− F (x))n−j
となる。また、i < j とすると (X(i), X(j))の同時密度関数は x < y に対しては
f(X(i),X(j))(x, y) =n!
(i− 1)!1!(j − 1− i)!1!(n− j)!F (x)i−1f(x)(F (y)− F (x))j−1−if(y)(1− F (y))n−j ,
x ≥ y に対しては f(X(i),X(j))(x, y) = 0となる。
独立な U1, . . . , Un ∼ U(0, 1) の順序統計量 U(k) がベータ分布 Beta(k, n− k + 1)に従うこと等は前期に示
した (cf . 前期の講義ノート pp.20–21)。また、これは二項母集団の母比率の区間推定 (精密法) に応用された
(cf . 前期の講義ノート pp.27–28, 検定は p.36)。
問題 2.13 ある人の射撃命中率は 70% であるといわれている。この人のある日の成績は 6発中 2発が命中し
たのみであった。この日はとくに命中率が低かったといえるか。5%の有意水準で検定せよ。
ここでは、指数分布に従う母集団について考察する。
指数分布に従う母集団の母平均の区間推定
• 指数分布 Ex(λ), ガンマ分布 Γ(a, λ), カイ二乗分布 χ2n, ベータ分布 Beta(a, b)の復習
X ∼ Γ(a, λ)とは X の密度関数が fX(x) =
λa
Γ(a)xa−1e−λx (x > 0)
0 (x ≤ 0)とする。
(a) 指数分布 Ex(λ) とガンマ分布 Γ(1, λ) は同じ分布、カイ二乗分布 χ2n とガンマ分布 Γ(n2 ,
12 ) は同じ分布。
(b) X,Y が独立で、X ∼ Γ(a, λ), Y ∼ Γ(b, λ)とし、c > 0を定数とすると次が成立:
(i) cX ∼ Γ(a, λ/c), (ii) X + Y ∼ Γ(a+ b, λ), (iii)X
X + Y∼ Beta(a, b),
(iv) a = m/2, b = n/2ならW =Y/n
X/m∼ Fnm であるから、
X
X + Y=
1
1 + Y/X=
1
1 + nmW
と変形すれば、F分布 Fnm とベータ分布 Beta(m2 ,n2 )の関係も従う。
命題 2.5 X1, . . . , Xn は i.i.d.で Ex(1/λ)に従うとき、2
λ(X1 + · · ·+Xn) ∼ χ2
n.
証明は上記 (a)の前半と、(b)の (i), (ii)より容易。
この命題より Sn = X1 + · · ·+Xn とすると、
P(χ22n
(1− 1
2ε)≤ 2
λSn ≤ χ2
2n
(12ε))
= 1− ε
となるから、括弧内の不等式を λについて解き、Sn = X1 + · · ·+Xn の実現値 sに置き換えることで、平均
値 λの信頼係数 1− εの信頼区間が 2s
χ22n(
12ε)≤ λ ≤ 2s
χ22n(1− 1
2ε)となることがわかる。
また、検定も同様にできる。帰無仮説 H0 : λ = λ0, 対立仮説 H1 : λ = λ0 を有意水準 εで検定するには、
実現値を s = x1 + · · ·+ xn とするとき、
2
λ0s < χ2
2n(1− ε/2)または χ22n(ε/2) <
2
λ0sのとき H0 を棄却し、 (2.6)
15
χ22n(1− ε/2) <
2
λ0s < χ2
2n(ε/2)のとき H0 を採択する。 (2.7)
片側検定の場合も同様に考えることができる。(各自考えよ。)
次に、データが n個で打ち切られた場合を考える。
命題 2.6 X1, . . . , XN は i.i.d.で Ex(1/λ)に従うとし、j 番目の順序統計量を Tj = X(j), j = 1, . . . , n(< N),
とする。
(1) (T1, T2, · · · , Tn)の密度関数 f(t1, t2, · · · , tn)は
f(t1, · · · , tn) =
N !
(N − n)!1
λnexp[− 1
λ
{ n∑j=1
tj + (N − n)tn}]
(0 < t1 < t2 < · · · < tn)
0 (その他)
(2.8)
(2) (2.8)を尤度関数とする λの最尤推定量は λ =1
n
{ n∑j=1
Tj + (N − n)Tn}.
(3) (2)の λに対し、2n · λλ∼ χ2
2n.
略証: (1) P (Tk ≥ tn) = e−1λ tn と順序統計量の定義から容易に示せる。
(2) 対数尤度関数 l(λ) = log f(t1, · · · , tn)を λに関し微分せよ。
(3) s < 1/2に対し E[es(2nλ/λ)] = (1− 2s)−n (χ22n 分布の積率母関数)を示せばよい。
E[es(2nλ/λ)] = E[e2sλ {
∑nj=1 Tj+(N−n)Tn}]
=N !
(N − n)!1
λn
∫ ∞0
e−1−2s
λ t1 dt1
∫ ∞t1
e−1−2s
λ t2 dt2 · · ·∫ ∞tn−2
e−1−2s
λ tn−1 dtn−1
∫ ∞tn−1
e−1−2s
λ (N−n+1)tn dtn
= · · · = (1− 2s)−n. (n = 4 < N として各自証明せよ。) □
この命題より S = T1 + · · ·+ Tn + (N − n)Tn とすると、
P(χ22n
(1− 1
2ε)≤ 2
λS ≤ χ2
2n
(12ε))
= 1− ε.
よって、S = T1 + · · ·+ Tn + (N − n)Tn の実現値をとすると、平均値 λの信頼係数 1− εの信頼区間が2s
χ22n(
12ε)≤ λ ≤ 2s
χ22n(1− 1
2ε)(2.9)
となることがわかる。
検定についても同様にできる。帰無仮説 H0 : λ = λ0, 対立仮説 H1 : λ = λ0 を有意水準 ε で検定するに
は、実現値 s = t1 + · · ·+ tn + (N − n)tn に対し、(2.6), (2.7)とすればよい。
データが X で打ち切られた場合 (寿命試験では定時打ち切りという)
大きさ N の標本がその値の X 以下であるもののみが測定され、観測された n個の値は x1, · · · , xn であったとする。このときは s = x1 + · · ·+ xn + (N − n)X をその実現値として、
平均値 λの信頼係数 1− εの信頼区間は (2.9)となる。
帰無仮説 H0 : λ = λ0, 対立仮説 H1 : λ = λ0 を有意水準 εで検定するには、(2.6), (2.7)とすればよい。
問題 2.14 以下の問いに答えよ。(1), (2)は小数点以下第 1位まで求めよ。
(1) ある電気商品の寿命を調べたところ、次の通りであった。平均寿命の 95%信頼区間を求めよ。
(∗) 732, 838, 915, 1211, 1355, 1420, 1638 (単位 時間)
(2) 10個の商品について調べた結果、(1)の (∗)のような標本が得られた。3個の部品はいづれも 1638時間
以上もったものとすれば、平均寿命の 95%信頼区間を求めよ。
(3) 10個の商品について調べた結果、(1)の (∗)のような標本が得られ、残りの 3個の部品はいづれも 1700
時間以上もったものとする。この商品の平均寿命は 1400時間といってよいか。有意水準 5%で検定せよ。
16
3 確率過程
3.1 マルコフ連鎖とマルチンゲール
定義 3.1 {Xt}がマルコフ過程 (連鎖)であるとは、任意の 0 ≤ t1 < t2 < · · · < tn < tと x1, x2, · · · , xn, xに対して、
P (Xt ≤ x|Xt1 = x1, Xt2 = x2, · · · , Xtn = xn) = P (Xt ≤ x|Xtn = xn)
となることをいう。つまり、将来は現在のみに依存し、過去には依存しない確率過程をいう。
このとき、{Xt}は次のようになる:
P (X1 = x1, X2 = x2, · · · , Xn = xn)
= P (Xn = xn|X1 = x1, · · · , Xn−1 = xn−1)P (X1 = x1, · · · , Xn−1 = xn−1)
= P (Xn = xn|Xn−1 = xn−1)P (X1 = x1, · · · , Xn−1 = xn−1)
...
= P (Xn = xn|Xn−1 = xn−1) · · ·P (X3 = x3|X2 = x2)P (X2 = x2|X1 = x1).
推移確率 P (Xt+1 = y|Xt = x) = pxy, 推移確率密度関数 fXt+s|Xt(y|x) =
f(Xt,Xt+s)(x, y)
fXt(x)= p(s, x, y) が
t によらないマルコフ過程は時間的一様 (推移確率が定常) なマルコフ過程という。ここではこの場合のみを
扱う。
• 有限マルコフ連鎖 Xt が S = {1, 2, · · · ,m}に値をとるの場合を考える。(S を状態空間という。)
pij = P (Xt+1 = j|Xt = i)とし、P =
p11 · · · p1m...
...
pm1 · · · pmm
とおく。P を推移確率行列という。
npij = P (Xt+n = j|Xt = i)とすると、Pn =
np11 · · · np1m...
...
npm1 · · · npmm
となる。特に、チャップマン-コロモゴ
ロフの方程式 n1+n2pij =
m∑k=1
n1pik n2pkj が成り立つ。これは
P (Xt+n1+n2 = j|Xt = i) =
m∑k=1
P (Xt+n1 = k|Xt = i)P (Xt+n1+n2 = j|Xt+n1 = k)
から容易にわかる。また、P (X0 = i) = pi とすると、
(P (Xn = 1), · · · , P (Xn = m)) = (P (Xn−1 = 1), · · · , P (Xn−1 = m))P = (p1, · · · , pm)Pn.
特に、π = (π1, · · · , πm)が π = πP を満たすとき、π を定常分布という。有限マルコフ連鎖では、∀i, j に対してある nがあって npij > 0とできる (このときマルコフ連鎖は既約であるという)とき定常分布はただ一つ
であり、特にある n0 があって ∀i, j に対し n0pij > 0を満たせば limn→∞
Pn =
π1 · · · πm
· · ·π1 · · · πm
となる。• 2重マルコフ過程
P (Xt+1 ≤ x|X0 = x0, X1 = x1, · · · , Xt−1 = xt−1, Xt = xt) = P (Xt+1 ≤ x|Xt = xt, Xt−1 = xt−1)
が成立しているとき、Xt を 2重マルコフ過程という。この場合は Yt = (Xt, Xt+1)とペアにして考えると Yt
はマルコフ過程となる。
17
例題 3.1 Xt を {1, 2, 3} 上のマルコフ連鎖でその推移確率を P =
0 1 0
0 1/2 1/2
1/2 0 1/2
とする。P (Xn = i|X0 = j) = npji とするとき、(1) 2p11, 3p11 を求めよ。(2) np11 を求めよ。
解: (1) 2p11 = p11p11 + p12p21 + p13p31 = 0 · 0 + 1 · 0 + 0 · 12= 0,
また、2p21 = 0 · 0 + 1
2· 0 + 1
2· 12=
1
4, 2p31 =
1
2· 0 + 0 · 0 + 1
2· 12=
1
4より、
3p11 = p11 · 2p11 + p12 · 2p21 + p13 · 2p31 = 0 · 0 + 1 · 14+ 0 · 1
4=
1
4.
(2) P の固有値は 1, i/2, −i/2. よって、すべての固有値で重複度 1より P は対角化可能なので、正則行列 U
があって、P = UDU−1 とできる*2。ただし、D は 1, i/2, −i/2を対角成分とする対角行列である。よって、Pn = UDnU−1 で Dn は 1, (i/2)n, (−i/2)n を対角成分とする対角行列なので、定数 a, b, cがあって、
np11 = a+ b( i2
)n+ c(− i2
)nと表せる。ここでは、np11 は実数で、±i = e±iπ/2 より(
± i2
)n=(12
)n (e±iπ/2
)n=(12
)n(cos
nπ
2± i sin nπ
2
)となるので、定数 α, β, γ があって、
np11 = α+(12
)n{β cos
nπ
2+ γ sin
nπ
2
}. (3.1)
と表せる。ここで、P 0 = (0pij)は単位行列、P = (1pij)に注意する。よって、
1 = 0p11 = α+ β, 0 = 1p11 = α+1
2γ, 0 = 2p11 = α− 1
4β
となるので α = 1/5, β = 4/5, γ = −2/5. よって、
P1(Xn = 1) = p(n)11 =
1
5+(12
)n{45cos
nπ
2− 2
5sin
nπ
2
}. □
問題 3.1 例 3.1のマルコフ連鎖 Xt に対して (1) 2p22, 2p23, 3p23 を求めよ。(2) np23 を求めよ。
ヒント: np23 も (3.1)の右辺のように表せる。
例題 3.2 次の推移行列をもつ {1, 2, 3, 4} 上のマルコフ連鎖の定常分布を求めよ。((b), (c) は問題 3.2 と
する。)
(a)
0.4 0.6 0 0
0 0.7 0.3 0
0.1 0 0.4 0.5
0.5 0 0 0.5
(b)
0.7 0.3 0 0
0.2 0.5 0.3 0
0 0.3 0.6 0.1
0 0 0.2 0.8
(c)
0.7 0 0.3 0
0.2 0.5 0.3 0
0.1 0.2 0.4 0.3
0 0.4 0 0.6
解: 定常分布を π = (π1, π2, π3, π4)とする。πP = π を書き直して、
0.4π1 + 0.1π3 + 0.5π4 = π1, 0.6π1 + 0.7π2 = π2, 0.3π2 + 0.4π3 = π3, 0.5π3 + 0.5π4 = π4,
最初の 3つの式より π1 = π2/2 = π3 = π4を得る。従って、π1+π2+π3+π4 = 1より、π = (1/5, 2/5, 1/5, 1/5)
となる。 □
問題 3.2 例 3.2の (b), (c)を推移行列とするマルコフ連鎖の定常分布をそれぞれ求めよ。
*2 ここでは対角化やジョルダン標準形における正則行列 U の存在だけ用い、その具体形を必要としない計算法を紹介する。
18
定義 3.2 {Xt}がマルチンゲールであるとは、すべての tについて
E[Xt+1|Xt, Xt−1, · · · , X2, X1] = Xt, a.s. (3.2)
となることをいう。
マルチンゲールを扱うのに次の条件つき期待値についての定理が有用である (cf . 前期の定理 1.7)。
定理 3.1 条件つき期待値について以下が成立する。ただし、X は多次元確率変数であってもよい。
(1) E[E[Y |X]] = E[Y ],
(2) E[aY + bZ|X] = aE[Y |X] + bE[Z|X], a.s. (a, bは定数。)
(3) E[g(X)h(Y )|X] = g(X)E[h(Y )|X], a.s. 特に E[g(X)|X] = g(X), a.s.
(4) X,Y が独立なら E[Y |X] = E[Y ], a.s.
(5) E[E[Z|X,Y ]|X] = E[Z|X], a.s.
{Xt}がマルチンゲールのとき、t > sに対して E[Xt|Xs, Xs−1, · · · , X1] = Xs となる。実際、
E[Xs+2|Xs, Xs−1, · · · , X1] = E[E[Xs+2|Xs+1, Xs, Xs−1, · · · , X1]|Xs, Xs−1, · · · , X1]
= E[Xs+1|Xs, Xs−1, · · · , X1] = Xs
を繰り返して証明される。ここで、最初の等号では定理 3.1 (5)を次の等号は (3.2)を t = s+ 1とし次の等号
は t = sとし用いた。また、その両辺の期待値をとれば、定理 3.1(1)より E[Xt] = E[Xs]が従う。
マルチンゲールの意味: 時刻 s までの情報のもとで未来 (t > s) の期待財産と時刻 s での財産が一致する、
すなわち、sから tまでの条件付き期待財産増分が 0であるときをいう。たとえば、公平な賭けを行っている
ギャンブラーの財産過程がその典型例である。
• 1次元対称ランダムウォークaを整数、ξ1, ξ2, · · · ,を i.i.d.で P (ξi = 1) = P (ξi = −1) = 1
2 とする。このとき、時間パラメータ tを 0以上
の整数として、Zat = a+ ξ1 + ξ2 + · · ·+ ξt を aを出発する 1次元対称ランダムウォークという。特に a = 0
のとき、Zt = Z0t と表し、(1次元)対称ランダムウォークという。
Zt は Z上のマルコフ連鎖で、その推移確率は P (Zt+1 = y|Zt = x) =
{12 y = x, x+ 1
0 その他となる。また、
E[Zt+1|Zt, Zt−1, · · · , Z1] = E[Zt|Zt, Zt−1, · · · , Z1] + E[ξt+1|Zt, Zt−1, · · · , Z1]
= Zt + E[ξt+1] = Zt
より対称ランダムウォークはマルチンゲールとなる。ここで、最初の等号は定理 3.1(1)を次の等号ではそれぞ
れ項に対して (3), (4)を用いた。
問題 3.3 Zt を対称ランダムウォークとするとき、次を求めよ。(0 < s < t, a, bを自然数とする。)
(1) E[Zt − Zs], (2) V (Zt − Zs), (3) P (Za+b = a− b), (4) P (Z5 = 3, Z8 = 2)
ヒント: (4)は Zt の独立増分性を用いよ。
問題 3.4 ξ1, ξ2, · · · ,を i.i.d.で P (ξi = 1) = P (ξi = −1) = 12 とし、Zt = ξ1 + ξ2 + · · ·+ ξt とするとき、
E[Z2t+1 − (t+ 1)|ξt, ξt−1, · · · , ξ1] = Z2
t − t.
を示せ。(これを Z2t − tが ξt, ξt−1, · · · , ξ1 に関してマルチンゲールであるという。)
「藤田岳彦 著 確率・統計・モデリング問題集」p.36– の練習問題 3.1–9, 3.20–も上記の知識で解けます。特
に、練習問題 3.25–は漸化式に関する問題なので是非挑戦ください。
19
3.2 ポアソン過程
定義 3.3 Nt, t ≥ 0がパラメータ (強度)λのポアソン過程であるとは、
(1) N0 = 0で Nt ∼ Po(λt) (t ≥ 0), つまり P (Nt = k) =(λt)k
k!e−λt, k = 0, 1, 2, . . ..
(2) (独立増分性) 0 < t1 < t2 < · · · < tn として、Nt1 , Nt2 −Nt1 , · · · , Ntn −Ntn−1 は独立。
(3) (定常増分性) t > s > 0として、Nt −Ns ∼ Nt−s.が成立するときにいう。
モデリングの教科書 p.3-7にあるように (1)の代わりに、
(1’) (微笑時間でジャンプの起こる確率は小さい) P (Nt+h−Nt ≥ 2) = o(h),すなわち、limh→0
Nt+h −Nth
= 0
を仮定してもよい。((1’), (2), (3)から (1)の導出は「藤田岳彦 著 確率・統計・モデリング問題集」p.33–を参
照のこと。)
例題 3.3 Nt をパラメータ λ のポアソン過程とする。次を求めよ。ただし、0 < s < t, k, l = 0, 1, 2, · · · とする。
(1) P (Ns = k,Nt = k + l), (2) E[NsNt], (3) Cov(Ns, Nt).
解: Ns と Nt −Ns が独立であることと Ns ∼ Po(λs), Nt −Ns ∼ Nt−s ∼ Po(λ(t− s))を用いる。
(1) P (Ns = k,Nt = k+l) = P (Ns = k,Nt−Ns = l) = P (Ns = k)P (Nt−Ns = l) =(λs)k
k!
(λ(t− s))l
l!e−λt.
(2) N ∼ Po(λ)のとき E[N ] = V (N) = λに注意する。
E[NtNs] = E[(Nt −Ns)Ns] + E[N2s ] = E[Nt −Ns]E[Ns] + V (Ns) + (E[Ns])
2
= λ(t− s)λs+ λs+ (λs)2 = λ2st+ λs.
(3) Cov(Ns, Nt) = E[NsNt]− E[Ns]E[Nt] = λ2st+ λs− λsλt = λs. □
例 3.4 Nt をパラメータ λのポアソン過程とすると、Nt − λtはマルチンゲールとなる。
証明: s < tとする。独立増分性より、Nt −Ns と Nu, 0 ≤ u ≤ sは独立. よって、
E[Nt − λt|Nu, 0 ≤ u ≤ s] = E[Nt −Ns +Ns − λt|Nu, 0 ≤ u ≤ s]= E[Nt −Ns] +Ns − λt = λ(t− s) +Ns − λt = Ns − λs □
問題 3.5 Nt をパラメータ λのポアソン過程、0 < s < tとするとき次を求めよ。
(a) E[Nt(Nt + 1)], (b) E[N2sNt], (c) E[NsN
2t ], (d) E[Nt|Ns = 1], (e) E[Ns|Nt = 2]
• 計数過程としてのポアソン過程Ti = i番目のイベントが発生する時刻とし、nt を時刻 tまでに起こったイベントの発生回数とする。このとき
nt を計数過程 (counting process)といい、特にそのイベントの発生間隔 T1, T2 − T1, T3 − T2, · · · が独立で同分布のときの nt を再生過程という。
ポアソン過程は発生の間隔 T1, T2 − T1, T3 − T2, · · · が独立で指数分布 Ex(λ)の場合である。実際、nt ≤ kとは時刻 t までに k + 1 回目のイベントが起こっていない、すなわち、Tk+1 > t となるので、Tn ∼ Γ(n, λ)
より、
P (nt ≤ k) = P (Tk+1 > t) =
∫ ∞t
λk+1
k!xke−λx dx = (部分積分を繰り返す) =
k∑j=0
(λt)j
j!e−λt. (3.3)
問題 3.6 T1, T2 − T1, T3 − T2 が独立で指数分布 Ex(λ)に従うとき、P (3T1 < T2), P (3T1 < 2T2 < T3),
P (2T2 < T3)を求めよ。
20
Poisson分布とカイ二乗分布の関係式
Tk+1 ∼ Γ(k + 1, λ)のとき Y2 = 2λTk+1 ∼ χ22(k+1) より、X ∼ Po(λ)のとき (3.3)で t = 1とすると
P (X ≤ k) = P (Tk+1 > 1) = P (Y2 > 2λ)
となる。また、P (X ≥ k) = 1− P (X ≤ k − 1) = 1− P (Y1 > 2λ) = P (Y1 ≤ 2λ),
ここで、Y2 ∼ χ22k となる。これより、Poisson 分布に従う母集団の母平均の区間推定 (精密法) の式、
1
2nχ22s(1− 1
2ε) ≤ λ ≤1
2nχ22(s+1)(
12ε)を得る。ここで、sは標本を x1, . . . , xn, その総和が s =
n∑k=1
xk を表す。
検定は、帰無仮説 H0 : λ = λ0, 対立仮説 H1 : λ = λ0 を有意水準 εで検定するには、
2nλ0 < χ22s(1− ε/2)または 2nλ0 > χ2
2(s+1)(ε/2)のとき H0 を棄却する。
χ22s(1− ε/2) < 2nλ0 < χ2
2(s+1)(ε/2)のとき H0 を採択する。
3.3 ブラウン運動
Wt, t ≥ 0が (標準)ブラウン運動であるとは、
(1) W0 = 0でWt ∼ N(0, t) (t ≥ 0).
(2) (独立増分性) 0 < t1 < t2 < · · · < tn として、Wt1 , Wt2 −Wt1 , · · · , Wtn −Wtn−1 は独立。
(3) (定常増分性) t > s > 0として、Wt −Ws ∼Wt−s.
となるときにいう。このとき、t 7→Wt が (確率 1で)連続となることが知られている。
例題 3.5 Wt を標準ブラウン運動とする。次を求めよ。ただし、0 < s < tとする。
(1) Wt の密度関数 fWt(x), (2) (Ws,Wt)の同時密度関数 f(Ws,Wt)(x, y), (3) E[WsWt],
(4) E[Wt|Wu, 0 ≤ u ≤ s] =Ws, すなわち、Wt はマルチンゲールになる。
解: Ws とWt −Ws が独立であることとWs ∼ N(0, s), Wt −Ws ∼Wt−s ∼ N(0, t− s)を用いる。(1) fWt(x) =
1√2πt
e−x2
2t . (2) f(Ws,Wt−Ws)(x, z) =1√2πs
e−x2
2s1√
2π(t− s)e−
z2
2(t−s) より、
f(Ws,Wt)(x, y) =1√2πs
e−x2
2s1√
2π(t− s)e−
(y−x)2
2(t−s) .
(3) E[WsWt] = E[Ws(Wt −Ws) +W 2s ] = E[Ws]E[Wt −Ws] + E[W 2
s ] = s.
(4) 独立増分性より、Wt −Ws とWu, 0 ≤ u ≤ sは独立. よって、
E[Wt −Ws +Ws|Wu, 0 ≤ u ≤ s] = E[Wt −Ws|Wu, 0 ≤ u ≤ s] +Ws = E[Wt −Ws] +Ws =Ws □
問題 3.7 Wt を標準ブラウン運動とする。次を求めよ。ただし、0 < s < t, α ∈ Rとする。
(a) E[W 3t ], (b) E[W 4
t ], (c) E[eαWt ], (d) E[WteαWt ], (e) E[W 2
sW2t ],
(f) E[W 2t |Ws], (g) E[eαWt |Ws], (h) E[Wse
αWt ], (i) P (Wt > 2Ws > 0)
21
4 シミュレーション
確率変数 X の密度関数を f(x)とすると、
E[g(X)] =
∫ ∞−∞
g(x)f(x) dx
であるが、この定積分を解析的に求めるのではなく、モンテカルロ法によって数値解析的に求める。
[モンテカルロ法] X1, X2, . . .を i.i.d.で Xi ∼ X とすると、g(X1), g(X2), . . .も i.i.d.で g(Xi) ∼ g(X)で
あり、大数の法則により
limn→∞
1
n
n∑i=1
g(Xi) = E[g(X)]
となる。これより、X の分布をもつ乱数 x1, x2, . . .を十分多く (N 個として)発生できれば、∫ ∞−∞
g(x)f(x) dx ≒ 1
N
N∑i=1
g(xi)
とできる。実際、U ∼ U(0, 1)のとき E[eU2
] =
∫ 1
0
ex2
dxを 10個の U(0, 1)に従う乱数
0.0623, 0.3246, 0.4265, 0.2796, 0.0424, 0.5218, 0.6877, 0.9943, 0.2771, 0.1668
を用いて (順に u1, . . . , uN , N = 10と表す)、
E[eU2
] =
∫ 1
0
ex2
dx ≒ 1
N
N∑i=1
eu2i ≒ 1.311
と求めることができる。別の 10個の乱数で求めると 1.5859, 同様に 100個の乱数で求めると 1.43676 となる
(以上 Excelを用いた)。この積分の Rで数値計算の値は 1.46265 である。
誤差評価や収束するスピードについて、中心極限定理より1n
∑ni=1 g(Xi)− E[g(X)]√
V (g(X))n
は N(0, 1) で近似でき
ることに注意する。これより、1
n
n∑i=1
g(Xi)− E[g(X)] =√
V (g(X))n N(0, 1) + o
( 1√n
)であり、もし V (g(X))
を小さくできれば、この誤差評価が良くなることがわかる (cf . 分散減少法)。
4.1 確率変数を生成する技法
4.1.1 逆関数法
一様分布に従う確率変数 U から F (x)を分布関数とする確率変数 X を生成するには、
X = F←(U), ただし、F←(u) = inf{x ; F (x) ≥ u }
とすればよい。具体的には、連続型ならX のとり得る値となる x上で F (x)は連続で (狭義)単調増加なので、
F−1(u) = F←(u)となり、X = F−1(U)の分布関数は F (x)となる。実際、
P (X ≤ x) = P (F−1(U) ≤ x) = P (U ≤ F (x)) = F (x).
また、X が離散型で P (Xk = xk) = pk, k = 1, 2, · · · , N とするとき、以下のように X を定めればよい。
0 ≤ U ≤ p1 ←→ X = x1p1 < U ≤ p1 + p2 ←→ X = x2
p1 + p2 < U ≤ p1 + p2 + p3 ←→ X = x3...
...p1 + p2 + · · ·+ pN−1 < U ≤ p1 + p2 + · · ·+ pN = 1 ←→ X = xN
(4.1)
22
例 4.1 X ∼ Ex(1)に対し、これを一様分布 U ∼ U(0, 1)から逆関数法を用いて生成し、E[X]とE[(1+X)−1]
をシミュレートする。
X の分布関数は F (x) = 1− e−x (x > 0)であった。よって、F−1(u) = − log(1− u)に注意して、X = − log(1− U)とおく。
次に左のように一様乱数 U が与えられたと
する。U から逆関数法で X を生成し、これか
ら (1 +X)−1 を計算したのが下の 2行である。
U 0.84 0.30 0.76 0.35 0.16
X 1.8326 0.3566 1.4271 0.4307 0.1743
(1 +X)−1 0.3530 0.7371 0.4120 0.6990 0.8516
例えば、U = 0.84に対しては、アクチュアリー資格試験 数学 付表 V 自然対数表を用いて計算すると、
X = − log(1− 0.84) = − log1.6
10= 2.3026− 0.4700 = 1.8326,
(1 +X)−1 = (1 + 1.8326)−1 = 0.35303 · · ·
同様にして、他の U に対しても上記の表のような X, (1 +X)−1 を得る。
これより、E[X]をシミュレートすると、1
5(1.8326 + 0.3566 + 1.4271 + 0.4307 + 0.1743) = 0.84426、
E[(1 +X)−1]については、1
5(0.3530 + 0.7371 + 0.4120 + 0.6990 + 0.8516) = 0.61054となる。 □
問題 4.1 一様分布に従う確率変数 U から次の分布に従う確率変数 X を逆関数法により生成する。
右表の一様乱数 U から生成される X の実現値をそれぞれ求め、
それにより平均値 E[X]をシミュレートせよ。
1 2 3 4 5
U 0.05 0.30 0.52 0.35 0.16
(1) 平均 2の指数分布に従う確率変数 X.
(2) 右の確率分布表で与えられる確率変数 X.
X 1 2 3 4 5
P 0.23 0.17 0.22 0.18 0.20
問題 4.2 0 < p < 1とする。U ∼ U(0, 1)のとき、X = int( logU
log(1− p)
)とすると X ∼ Ge(p)を示せ。ここ
で、int(a)は aの整数部分を表す。次に、p = 0.2とし、問題 4.1で与えられた一様乱数 U からこの式によっ
て生成されるX の実現値をそれぞれ求め、それにより平均値 E[X]をシミュレートせよ。(注意 モデリングの
教科書 pp.4-22–には、効率的な離散型確率変数に従う乱数の発生方法が述べられている。)
4.1.2 棄却法
密度関数 g(y)である確率変数 Y から、密度関数 f(x)である確率変 X を生成したいとする。この際、
ある定数 cがあって、f(y)/g(y) ≤ c for all y (4.2)
であるならば、次の手順により、g に従う Y を用いて f に従う確率変数 X を求めることができる。
[手順 1] 密度関数が g である確率変数 Y を生成する。(例えば、逆関数法で生成する。)
[手順 2] (0, 1)上の一様分布に従う確率変数 U を生成する。
[手順 3] U ≤ f(Y )
cg(Y )ならば、X = Y と (X として Y を採用)し、U >
f(Y )
cg(Y )ならば手順 1にもどる。
ここで cは (4.2)を満たす最小のものを選択することにより、計算手順数が効率化できることに注意する。
定理 4.1 棄却法で生成された X の密度関数は f .
証明: ∀y に対して 0 ≤ f(y)cg(y) ≤ 1 に注意する。
P (X ≤ x) = P(Y ≤ x
∣∣∣U ≤ f(Y )
cg(Y )
)=P(Y ≤ x,U ≤ f(Y )
cg(Y )
)P(U ≤ f(Y )
cg(Y )
) =
∫ x
−∞
∫ f(y)cg(y)
0
g(y) dudy∫ ∞−∞
∫ f(y)cg(y)
0
g(y) dudy
23
=
∫ x
−∞
f(y)
cg(y)g(y) dy∫ ∞
−∞
f(y)
cg(y)g(y) dy
=
1
c
∫ x
−∞f(y) dy
1
c
∫ ∞−∞
f(y) dy
=
∫ x
−∞f(y) dy. □
系 4.2 棄却法で、定理 4.1の証明より X = Y として採用される確率は P
(U ≤ f(Y )
cg(Y )
)=
1
cである。ここ
で、T で一つ採用されるまでの手順 1–3の繰り返し回数を表すと、P (T = k) =(1− 1
c
)k−1 1c, k = 1, 2, . . .,
となり、平均繰り返し回数は
E[T ] =
∞∑k=1
k(1− 1
c
)k−1 1c=
1c
{1− (1− 1c )}2
= c
となる。
注意 離散型の場合に棄却法で確率変数を生成する方法についてはモデリングの教科書 p.4-9を参照のこと。
例 4.2 Z ∼ N(0, 1)に対し、X = |Z|を Y ∼ Ex(1)から棄却法を用いて生成する。
X の密度関数は f(x) =2√2πe−x
2/2 (x > 0)であった。よって、Y の密度関数 g(x) = e−x (x > 0)で
f(x)
g(x)=
√2
πex−x
2/2 =
√2
πe−[(x−1)
2−1]/2 ≤ f(1)
g(1)=
√2e
π.
よって、c =√2e/π ととり、 f(x)
cg(x) = exp(x− x2
2 −12
)= exp
{− (x−1)2
2
}. 従って、
[手順 1] 独立な U1, U2 ∼ U(0, 1)を準備する。
[手順 2] (Y の逆関数法による生成)
FY (y) = 1− e−y より、F−1Y (u) = − log(1− u)に注意して、Y = − log(1− U1)とおく。
[手順 3] U2 ≤ exp{−(Y − 1)2/2
}なら X = Y を採用。そうでなければ手順 1にもどる。
次に左のように一様乱数 U1, U2 が与えられたとする。
U1 から逆関数法で Y を生成し、さらに棄却法で X
を生成する。
FY (y) = 1− e−y より F−1Y (u) = − log(1− u).まず、Y = − log(1− U1)を求め、記入する。
U1 0.84 0.30 0.76 0.35 0.16
U2 0.08 0.77 0.69 0.57 0.77
Y 1.8326 0.3566 1.4271 0.4307 0.1743
Y ′ 0.7046 0.8105 0.9139 0.8521 0.7117
次に、Y ′ = e−(Y−1)2/2 を求めた。例えば、U1 = 0.84に対しては、
Y = − log(1− 0.84) = − log1.6
10= 2.3026− 0.4700 = 1.8326,
Y ′ = e−(1.8326−1)2/2 = e−0.347··· = e−0.35 = (e−0.10)3e−0.05 = (0.9048)3 · 0.9512 = 0.7046
と、アクチュアリー資格試験 数学 付表 V 自然対数表および VI 指数関数表を用いて計算した*3。
同様にして、他の U1 に対しても上記の表のような Y, Y ′ を得る。最後に、Y ′ と U2 を比較して棄却される
データは最後の一つとなり、E[X]をシミュレートすると、1
4(1.8326 + 0.3566 + 1.4271 + 0.4307) = 1.01175
となる。 □
問題 4.3 X ∼ Γ( 32 , 1)に対し、Y ∼ Ex(λ) から棄却法を用いて生成する。
(1) 繰り返し回数が最小となる λを求めよ。また、1つ採用するための繰り返し回数の期待値を求めよ。
(2) (1)で求めた λに対して、例 4.2の一様乱数 U1, U2 から、逆関数法により Y を生成し、この Y から棄却
法で X を生成する。これにより、X の期待値をシミュレートせよ。
*3 e0.01, e−0.01 の値より Y ′ の値の誤差は 1% 程度となる。
24
問題 4.4 f(x) = 30x(1− x)4, 0 < x < 1を密度関数にもつ確率変数を Y ∼ U(0, 1)から棄却法で生成する。
(1) 繰り返し回数を最小にする定数 cを定めよ。
(2) X を生成する手順を述べ、右表の Y および U ∼ U(0, 1)の
シミュレーション結果から生成される X の平均値を求めよ。
Y 0.05 0.30 0.52 0.35 0.16
U 0.08 0.77 0.69 0.57 0.77
4.1.3 合成法
F1, · · · , Fn を分布関数で各々に対応する確率変数を生成する手段を持っているとする。このとき、F がその重み付け平均
F (x) =
n∑i=1
hiFi(x),
n∑i=1
hi = 1, hi > 0
を分布関数とする確率変数 Y は次のようにシミュレーションする。(簡単のため n = 3で述べる。)
独立な U1, U2 ∼ U(0, 1)を準備し、
0 ≤ U1 ≤ h1 ←→ Y = F−11 (U2)h1 < U1 ≤ h1 + h2 ←→ Y = F−12 (U2)
h1 + h2 < U1 ≤ h1 + h2 + h3 = 1 ←→ Y = F−13 (U2)
と Y を定める。ただし、Fi が離散型の場合は (4.1)のように定義するものとする。
証明: 0 ≤ h < k ≤ 1に対し P (h ≤ U1 ≤ k, F−1i (U2) ≤ x) = P (h ≤ U1 ≤ k)P (U2 ≤ Fi(x)) = (k − h)Fi(x)より
P (Y ≤ x) = P (0 ≤ U1 ≤ h1, F−11 (U2) ≤ x) + P (h1 < U1 ≤ h1 + h2, F−12 (U2) ≤ x)
+P (h1 + h2 < U1 ≤ 1, F−13 (U2) ≤ x)= h1F1(x) + h2F2(x) + h3F3(x). □
例題 4.3 fX(x) =
{3x3 + 1
2x (0 < x < 1)
0 (その他)となる確率変数 X を合成法を用いて作れ。
また、右の表の U1, U2 をシミュレーション結果として生成さ
れる X の平均値を求めよ。
U1 0.23 0.94 0.78 0.55 0.20
U2 0.94 0.56 0.62 0.92 0.91
解: F1(x) = x4, F2(x) = x2 とすると、FX(x) = 34F1(x) +
14F2(x) (0 < x < 1) となることに注意する。
よって、0 ≤ U1 ≤ 3/4なら X = U1/42 , 3/4 < U1 ≤ 1なら X = U
1/22 ととればよい。
U1 = 0.23について、0 ≤ U1 ≤ 1/4より、X = F−11 (U2) = (0.94)1/4 = 0.9846 · · · ,U1 = 0.94について、1/4 < U1 ≤ 1より、X = F−12 (U2) = (0.56)1/2 = 0.7483 · · · ,
同様にすると順に、X = (0.62)1/2 = 0.7874 · · · , X = (0.92)1/4 = 0.9793 · · · , X = (0.91)1/4 = 0.9766 · · ·となり、E[X]をシミュレートすると、
1
5(0.985 + 0.748 + 0.787 + 0.979 + 0.9766) = 0.853となる。 □
問題 4.5 合成法に関する以下の問に答えよ。ただし、X の実現値を求める際は、fX(x)の右辺の左から順に、
U1 の範囲を定めるものとし、U2 により逆関数法を用いて X の実現値を決定せよ。
(1) fX(x) =
{15x−1/2 + 4
5x+ 15 (0 < x < 1)
0 (その他)となる確率変数 X を合成法を用いて作れ。また、例題 4.3
の表の U1, U2 をシミュレーション結果として生成される X の平均値を小数点以下第 3位まで求めよ。
(2) g1(x) =
{2(x+ 1)−3 (x > 0)
0 (その他)とし、g2(x)を [0, 5]上の一様分布の密度関数とする。
fX(x) = 23g1(x) +
13g2(x)となる確率変数 X を合成法を用いて作れ。また、例題 4.3の表の U1, U2 をシミュ
レーション結果として生成される X の平均値を小数点以下第 2位まで求めよ。
(3) g1(x) =
{0.5e−0.5x (x > 0)
0 (その他), g2(x) =
{ex (x < 0)
0 (その他)とする。fX(x) = 2
3g1(x) +13g2(x)となる
25
確率変数X を合成法を用いて作れ。また、(U1, U2) = (0.32, 0.75), (0.81, 0.35)としてX の実現値をそれぞれ
求めよ。
4.2 分散減少法
この章の最初に述べたが、X1, X2, . . .を i.i.d.とし、g(Xi)の平均を θ, 標準偏差を σ とする。このとき、そ
の標本平均 θn =1
n
n∑i=1
g(Xi)は、E[θn] = E[g(X1)] = θ であり、
V (θn) =1
n2
n∑i=1
V (g(Xi)) =σ2
n, すなわち E[(θn − θ)2] =
σ2
n.
これより、θn と θ の標準誤差 (二乗平均誤差の平方根) は c/√nという形となる。これは nをシナリオ数とす
ると精度 c/√nとなる。即ち、1000本のシナリオを用いたモンテカルロ法により得られた値の精度をもう一
ケタ上げるためには、その 100倍である 10万本のシナリオが必要となる。
ここではこの cを減少させることにより精度をよくする方法である分散減少法を紹介する。(詳しくはモデ
リングの教科書 pp.4-31–を参照のこと。)
4.2.1 負の相関法
定理 4.3 f, g がともに単調増加もしくはともに単調減少なら、Cov(f(X), g(X)) ≥ 0.
証明: x ≤ y なら (f(x) − f(y))(g(x) − g(y)) ≥ 0, x ≥ y でも (f(x) − f(y))(g(x) − g(y)) ≥ 0 より、常に
(f(X)− f(Y ))(g(X)− g(Y )) ≥ 0は成立する。ここで、X,Y が独立で同じ分布に従うとすると、
0 ≤ E[(f(X)− f(Y ))(g(X)− g(Y ))] = E[f(X)g(X)]− E[f(X)g(Y )]− E[f(Y )g(X)] + E[f(Y )g(Y )]
= E[f(X)g(X)]− E[f(X)]E[g(Y )]− E[f(Y )]E[g(X)] + E[f(Y )g(Y )] (∵ X,Y は独立)
= E[f(X)g(X)]− E[f(X)]E[g(X)]− E[f(X)]E[g(X)] + E[f(X)g(X)] (∵ X,Y は同分布)
= 2(E[f(X)g(X)]− E[f(X)]E[g(X)]) = 2Cov(f(X), g(X)). □
これより、U ∼ U(0, 1)とし、h(x)が単調増加 (減少)なら、−h(1− x)も単調増加 (減少)なので、
Cov(h(U), h(1− U)) = −Cov(h(U),−h(1− U)) ≤ 0
となる。
X の分布関数を F とし、その逆関数 F−1 が存在するとする。Ui, i = 1, 2, . . ., を U(0, 1)に従う i.i.d.とす
ると、1− Ui ∼ U(0, 1)なので、
Xi = F−1(Ui), Yi = F−1(1− Ui)
とおくと、分布 F に従う確率変数 Xi, Yi を生成できる。ここで、
θ(1)2M =
1
2M
M∑i=1
{g(Xi) + g(Yi)} =1
2M
M∑i=1
{g(F−1(Ui)) + g(F−1(1− Ui))}
とすると、これは E[θ(1)2M ] = θを満たす。ただし、θは g(X)の平均であった。ここで、もし g(x)が単調なら、
g(F−1(u))も単調なので、Cov(g(F−1(Ui)), g(F−1(1− Ui))) ≤ 0となることに注意すると、
V (θ(1)2M ) =
1
(2M)2
M∑i=1
V [g(F−1(Ui)) + g(F−1(1− Ui))]
=1
(2M)2
M∑i=1
{V [g(F−1(Ui))] + 2Cov(g(F−1(Ui)), g(F
−1(1− Ui))) + V [g(F−1(1− Ui))]}
26
≤ 1
(2M)2
M∑i=1
{V [g(F−1(Ui))] + V [g(F−1(1− Ui))]
}=
1
2MV (g(X)) = V (θ2M ).
ここで、3行目の最初の等号は F−1(Ui), F−1(1 − Ui)ともに X と同じ分布に従うことを用いた。以上より、
θ(1)2M は θ の不偏推定量であり、標本平均 θ2M より分散が小さくなることがわかった。また、一様分布に従う
確率変数を 2M 個生成する必要はなく、M 個のみ生成すればよいので、一様分布に従う確率変数を生成する
時間も節約できる。
例題 4.4
∫ 1
0
sinπ
2x dxを、(a) 2M 個の一様乱数 U1, U2, · · · , U2M を用いてモンテカルロシミュレーション
をした場合の誤差の分散 V (θ2M )を求めよ。
(b) 2M 個の一様乱数 U1, U2, · · · , UM , 1− U1, 1− U2, · · · , 1− UM を用いた誤差の分散 V (θ(1)2M )を求めよ。
(c) 問題 4.1の一様乱数 U を用いて、θM と θ(1)2M , M = 5, の実現値を求めよ。
解: 以下、独立性を考えない式では Ui を U と表す。g(x) = sinπ
2xとおく。θ2M =
1
2M
2M∑i=1
g(Ui)について、
E[g(Ui)] =
∫ 1
0
sinπ
2x dx =
2
π, E[g(Ui)
2] =
∫ 1
0
sin2π
2x dx =
1
2より
V (θ2M ) =1
(2M)2
2M∑i=1
V (g(Ui)) =1
2MV (g(U)) =
1
2M
(12− 4
π2
)≒ 1
2M0.09471565 · · · .
θ(1)2M =
1
2M
M∑i=1
{g(Ui) + g(1− Ui)}について、
E[{g(Ui) + g(1− Ui)}2] = E[g(Ui)2] + E[g(1− Ui)2] + 2E[g(Ui)g(1− Ui)]
=1
2+
1
2+
∫ 1
0
2 sinπ
2x sin
π
2(1− x) dx = 1 +
∫ 1
0
2 sinπ
2x cos
π
2x dx = 1 +
2
πより
V (θ(1)2M ) =
1
(2M)2
M∑i=1
V (g(Ui) + g(1− Ui)) =1
4MV (g(U) + g(1− U))
=1
4M
[E[{g(U) + g(1− U)}2]− (E[g(U) + g(1− U)])2
]=
1
4M
(1 +
2
π− 16
π2
)≒ 1
2M0.007740417054 · · · .
右表のより、θ5 =2.084
5= 0.417,
θ(1)2·5 =
2.084 + 4.362
10= 0.645. □
U 0.05 0.30 0.54 0.35 0.18 計
sin π2U 0.078 0.454 0.750 0.522 0.279 2.084
sin π2 (1− U) 0.997 0.891 0.661 0.853 0.960 4.362
問題 4.6
∫ 1
0
√x dxを、(a) 2M 個の一様乱数 U1, U2, · · · , U2M を用いてモンテカルロシミュレーションをし
た場合の誤差の分散 V (θ2M )を求めよ。
(b) 2M 個の一様乱数 U1, U2, · · · , UM , 1− U1, 1− U2, · · · , 1− UM を用いた誤差の分散 V (θ(1)2M )を求めよ。
(c) 問題 4.1の一様乱数 U を用いて、θM と θ(1)2M , M = 5, の実現値を求めよ。
27
4.2.2 制御変量法
θ(3)n =
1
n
n∑i=1
{g(Xi)+c(h(Xi)−θh)}とし、E[h(Xi)] = θhは既知と仮定すると (この hを制御変量という)、
V (θ(3)n ) =1
nV [g(X) + c(h(X)− θh)] =
1
n
{V (g(X)) + 2cCov(g(X), h(X)) + c2V (h(X))
}=
1
n
[V (h(X))
{c+
Cov(g(X), h(X))
V (h(X))
}2
+ V (g(X))− (Cov(g(X), h(X)))2
V (h(X))
](4.3)
より、c = c∗ = −Cov(g(X), h(X))
V (h(X))ととると、分散は最小となり、このとき
V (θ(3)n ) =1
n
{V (g(X))− (Cov(g(X), h(X)))2
V (h(X))
}≤ 1
nV (g(X)) = V (θn) となる。
ここで、V (θ
(3)n )
V (θn)= 1− (Cov(g(X), h(X)))2
V (h(X))V (g(X))= 1− (Corr(g(X), h(X)))2,
つまり分散の減少率は g(X)と h(X)の相関係数 Corr(g(X), h(X))の 2乗となる。
一方、標本値だけでは、
Cov[g(X), h(X)] =1
n− 1
n∑i=1
(g(Xi)− g(X))(h(Xi)− h(X)),
Var[h(X)] =1
n− 1
n∑i=1
(h(Xi)− h(X))2, ここで、 g(X) =1
n
n∑i=1
g(Xi), h(X) =1
n
n∑i=1
h(Xi).
cを c∗ = − Cov[g(X), h(X)]
Var[h(X)]で推定する。
以上より、θ(3)n = g(X)− Cov[g(X), h(X)]
Var[h(X)](h(X)− θh)がその実現値となる。
注意 4.1 (回帰分析との関係) θ(3)n = g(X)− Cov[g(X), h(X)]
Var[h(X)](h(X)− θh)は、回帰分析の式 y = α+βxの
式と関係づけると覚えやすい。ここで、
β =sxys2x
, α = y − βx
であった。yi = g(Xi), xi = h(Xi) とし y = θ(3)n , x = θh
と関連付けると、
y = g(X) =1
n
n∑j=1
g(Xj), sxy = Cov[g(X), h(X)] =1
n− 1
n∑j=1
{h(Xj)− h(X)}{g(Xj)− g(X)}
x = h(X) =1
n
n∑j=1
h(Xj), s2x = Var[h(X)] =1
n− 1
n∑j=1
{h(Xj)− h(X)}2
と対応させることができる。これより、
β =sxys2x
=Cov[g(X), h(X)]
Var[h(X)]= −c∗, α = y − βx = g(X) + c∗h(X)
に注意すると、
θ(3)n = α+ βθh = g(X) + c∗h(X)− c∗θh = g(X)− Cov[g(X), h(X)]
Var[h(X)](h(X)− θh)
を得る。
28
例題 4.5
∫ 1
0
ex dxを、(a) n個の一様乱数 U1, U2, · · · , Un を用いてモンテカルロシミュレーションをした場
合の誤差の分散 V (θn) を求めよ。(b) 制御変量 h(U) = U をとったときの誤差の分散 V (θ(3)n ) を計算せよ。
(c) 問題 4.1の一様乱数 U の標本値を用いて、θn と θ(3)n , n = 5, の実現値を求めよ。
解: g(U) = eU とおく。θn =1
n
n∑i=1
g(Ui)について、
E[g(Ui)] =
∫ 1
0
eu du = e− 1, E[g(Ui)2] =
∫ 1
0
e2u du =e2 − 1
2より
V (θn) =1
nV (g(U)) =
1
n
(e2 − 1
2− (e− 1)2
)=
1
n
(e− 1)(3− e)2
≒ 1
n0.242035607 · · · .
θ(3)n =
1
n
n∑i=1
{g(Ui) + c(h(U1)− θh)}について、(4.3)より c = −Cov(g(U), h(U))
V (h(U))とすれば、
V (h(U)) = E[U2]− (E[U ])2 =1
12, Cov(g(U), h(U)) =
∫ 1
0
xex dx− 1
2
∫ 1
0
ex dx = 1− e− 1
2より
V (θ(3)n ) =1
n
[V (g(U))− (Cov(g(U), h(U)))2
V (h(U))
]=
1
n
[e2 − 1
2− (e− 1)2 − 12
(1− e− 1
2
)2]=
1
n
(3− e)(7e− 19)
2=
1
n0.0394022292 · · · .
標本値では、Var(U) = 0.03393, Cov[eU , U ] = 0.046032より、c∗ = −1.3567, θh = 0.5 (既知) より次の表を
得る。よって、θ5 =6.733
5= 1.3466,
θ(3)5 = θ5 + c∗(h(U)− θh)= 1.3466− 1.3567(0.284− 0.5) = 1.6396. □
U 0.05 0.30 0.54 0.35 0.18 計
eU 1.051 1.350 1.716 1.419 1.197 6.733
問題 4.7
∫ 1
0
√1− x2 dxを (a) n個の一様乱数 U1, U2, · · · , Un を用いてモンテカルロシミュレーションをし
た場合の誤差の分散 V (θn)を求めよ。
(b) 制御変量 h(U) = U をとったときの誤差の分散 V (θ(3)n )を計算せよ。
(c) 制御変量 h(U) = U2 をとったときの誤差の分散 V (θ(3)n )を計算せよ。
(d) 問題 4.1の一様乱数 U を用いて、θ5 と (b)の θ(3)5 の実現値を求めよ。
4.3 統計の復習 3 適合度、独立性の検定
4.3.1 適合度の検定
n個のデータが d個の階級 A1, . . . , Ad に分類できたとき、Aj が出現する母比率を pj とし、次を考える。
帰無仮説 H0:(p1, · · · , pd) = (p01, · · · , p0d), 対立仮説 H1:(p1, · · · , pd) = (p01, · · · , p0d).
このとき、Aj の出現回数を nj とし、統計量
Tn =
d∑j=1
(nj − np0j )2
np0j=∑ (観察度数−期待度数)2
期待度数(4.4)
を考える。このとき、H0 のもと Tn は自由度 d− 1の χ2 分布に法則収束する (確率統計学 IIで示す)から、も
し有意水準 αで検定するのであれば、(4.4)に実現値を代入した値を tとするとき、
もし t < χ2d−1(α)であれば H0 を採択し、
もし t ≥ χ2d−1(α)であれば H0 を棄却する。
29
ただし、np0j < 5となる階級があるときは、隣接の階級と合併しすべての階級で np0j ≥ 5となるように分けな
おす。
例題 4.6 ある県の成人を母集団とし、無作為抽出された 1000人の血液型を調べたところ下の表のような観測値を得た。この県の血液型の分布は日本人の血液型の分布と同じであ
るといえるか。有意水準 5%で検定せよ。ただし、日本人の血液型の
分布は A : B : O : AB = 38.0 : 21.8 : 30.8 : 9.4である。
血液型 A B O AB
観測値 360 216 330 94
解: A, B, O, AB の比率を順に p1, p2, p3, p4 とし、次の仮説を設定する。
H0 : (p1, p2, p3, p4) = (0.380, 0.218, 0.308, 0.094), H1 は H0 の否定.
有意水準 5%であるから、棄却域は自由度が 4− 1であることに注意して
T ≥ χ23(0.05) = 7.8147
である。期待度数は A, B, O, ABの順に 380, 218, 308, 94であらから、実現値を代入して、
t =(360− 380)2
380+
(216− 218)2
218+
(330− 308)2
308+
(94− 94)2
94= 2.6424 · · · .
この値は棄却域に入らないから、H0 は採択される。従って、この県の血液型の分布は日本人の血液型の分布
に同じといえる。 □
問題 4.8 ある国の成人を母集団とし、無作為抽出された 100人の血液型を調べたところ下の表のような観測値を得た。この県の血液型の分布は日本人の血液型の分布と同じであ
るといえるか。有意水準 5%で検定せよ。ただし、日本人の血液型の
分布は A : B : O : AB = 38.0 : 21.8 : 30.8 : 9.4である。
血液型 A B O AB
観測値 28 23 42 7
問題 4.9 次のデータは 4枚のコインを同時に投げ表の枚数を 100回行い記録したものである。
これは二項分布 B(4, 1/2)に従っているといえるか。有意水準 5%
で検定せよ。
表の枚数 0 1 2 3 4
観測値 5 23 33 30 9
4.3.2 独立性の検定
N 個のデータが 2種類の属性 A,B によるそれぞれの各階級 A1.A2, · · · , Ar および B1, B2, · · · , Bs に分割されて、右の表のような度数氷河できたものとする。これを r × s分割表という。ここで、Ai ∩Bj なる性質をもったデータ数は fij とする。このとき、次を考える。
帰無仮説 H0 : 二つの属性 A,B は独立である
ことを有意水準 α で検定するには、すべての度数 fij ≥ 5 かつ
fi•f•j/N ≥ 5であるとき、
χ2 =
r∑i=1
s∑j=1
(fij −
fi•f•jN
)2fi•f•jN
> χ2ϕ(α) ならば H0 を棄却し、
χ2 < χ2ϕ(α) ならば H0 を採択することにすればよい。ここで
ϕ = (r − 1)(s− 1)である。
AB
B1 B2 · · · Bs 計
A1 f11 f12 · · · f1s f1•
A2 f21 f22 · · · f2s f2•...
......
. . ....
...
Ar fr1 fr2 · · · frs fr•
計 f•1 f•2 · · · f•s N
注意 4.2 (1) H0 が正しいとして分割表を作れば、Ai ∩ Bj のクラスに入るデータ数の期待値は fi•f•j/N で
ある。従って、統計量 χ2 は
χ2 =∑i
∑j
(実現値−期待度数)2
期待度数
30
と (4.4)と同様となる。ただし、自由度は ϕ = (r − 1)(s− 1)である。
(2) 2× 2分割表の簡便計算法: 2× 2分割表のときには、上の χ2 の値は次のようにして計算すると便利であ
る (第 2の等式が成立することを確かめよ)。
χ2 =
r∑i=1
s∑j=1
(fij −
fi•f•jN
)2fi•f•jN
=(f11f22 − f12f21)2N
f1•f2•f•1f•2. (4.5)
(3) 2× 2分割表におけるYatesの修正法: もし、(2)の χ2 の計算のときに、fij の中に 5より小さな値をと
るものがあったなら
χ2 =
(|f11f22 − f12f21| − N
2
)2N
f1•f2•f•1f•2,
とすればよい。このときも自由度はもちろん (2− 1)× (2− 1) = 1である。
(4) Fisher の直接確率計算法: (3) の検定において、度数 fij が小さすぎるとき χ2 分布は利用できない。
H0: A,B は独立である という仮説の下で、全ての周辺度数 fi•, f•j が一定のときの r × s分割表のように各
Ai ∩Bj の度数が fij となる確率が
∏ri=1 fi•!
∏sj=1 f•j !
N !∏ri=1
∏sj=1 fij !
であることを用い、実現値以上に偏った度数分布にな
るすべての場合の確率の和 pを計算し、有意水準と比較する。
例題 4.7 ある都市で有権者の A内閣に対する支持率を調べた。
有権者から男性 150人、女性 100人を抽出し、支持する、しないを調
べたらその人数は右の表のようであった。A内閣の支持率は男性と女
性とで、違いがあるとみてよいか、有意水準 5%で検定せよ。
男 女 計
支持する 75 60 135
支持しない 75 40 115
計 150 100 250
解: 次の仮説を設定する。H0 : 男性女性と支持するしないは関係ない, H1 は H0 の否定.
有意水準 5%であるから、棄却域は自由度が (2− 1)(2− 1) = 1であることに注意して
T ≥ χ21(0.05) = 3.8415
である。2× 2分割表の簡便計算法を用いて実現値を代入すると、
χ2 =(75 · 40− 60 · 75)2 · 250135 · 115 · 150 · 100
= 2.415 · · · .
これは棄却域に入らないから、H0 は採択される。従って、男性女性と支持不支持には関係がないといえる。□
問題 4.10 (1) ある病気の予防注射の効果を調べるために、300人を調査したところ、下の表 (1)の結果を得
た。予防注射の効果があるといえるか。有意水準 5%で検定せよ。
(2) 下の表 (2)は 300人の自動車所有者を年齢と過去 2年間に起こした事故数に応じて分類したものである。
年齢と事故数の間に関係があるかどうかを有意水準 5%で検定せよ。
表 (1)
発病 発病なし
予防注射あり 38 142
予防注射なし 58 62
表 (2)
事故数 0 1 ∼ 2 3
21歳以下 8 23 14
22 ∼ 26 21 42 12
27歳以上 71 90 19
31
5 損保数理に関する確率統計の話題から
5.1 最尤推定量の漸近挙動
ここでは、前期に学んだ最尤推定量に関する、極限定理を扱う。次を参照した。
[LC] E.L. Lehmann, G. Casella: Theory of Point Estimation, Second Edition, Springer, 1998
確率統計学 IIで学ぶ次の定義と定理を述べておく。この節ではこれが基本となる。
定義 5.1 (1) ∀δ > 0 に対して limn→∞
P (|Xn − X| < δ) = 1 となるとき、Xn は X に確率収束するといい、
Xn → X in prob. と表す。
(2) 確率変数 Y の分布関数を FY (y) = P (Y ≤ y)と表すとする。Xn が X に法則収束するとは、FX の任意
の連続点 cに対して、 limn→∞
FXn(c) = FX(c)となるときにいう。このとき、Xn → X in law と表す。
定理 5.1 X1, X2, . . .を i.i.d. とし、Sn =n∑i=1
Xi とする。
(大数の法則) limx→∞
xP (|X1| > x) = 0 ならば、数列 an があって1
nSn − an → 0 in prob. となる。特に、
E[|X1|] <∞であれば、1
nSn → E[X1] in prob. となる。(正確には大数の弱法則という。)
(中心極限定理) E[X12] < ∞ とする。µ = E[X1], σ
2 = V (X1) とすると、Sn − nµ√
n→ N(0, σ2) in law と
なる。
θ は未知母数とし Θをその母数空間とする。X1, X2, . . . , Xn は i.i.d. でその密度関数を f(x|θ)あるいは確率関数を p(x|θ) = Pθ(X = x)とする。以下、連続型の場合を考え、その確率を Pθ, 期待値を Eθ で表すもの
とする。(離散型の場合も同様に示せる。)さらに次を仮定する。
(A1) θ = θ′ なら密度関数として f(x|θ) = f(x|θ′).(A2) A = {x | f(x|θ) > 0 } は θ によらない。
(A3) 母数空間 Θは開集合とし、真のパラメータ θ0はその内点であるとする。
このとき、尤度関数を Ln(θ|x) =∏ni=1 f(xi|θ), 対数尤度を ln(θ|x) =
∑ni=1 log f(xi|θ) と表す。ただし、
x = (x1, x2, . . .)とする。この尤度関数もしくは対数尤度を最大にする θn = θn(x1, . . . , xn):
Ln(θn|x) = maxθ∈Θ
Ln(θ|x) または ln(θn|x) = maxθ∈Θ
ln(θ|x)
に対して、θn = θn(X1, . . . , Xn)が最尤推定量であった。
定理 5.2 (A1)–(A3)の条件のもと、∀θ = θ0 に対して、
Pθ0(Ln(θ0|X) > Ln(θ|X))→ 1, as n→∞,
ただしX = (X1, X2, . . .).
証明: Ln(θ0|X) > Ln(θ|X) ⇐⇒ 1
n
n∑i=1
logf(Xi|θ)f(Xi|θ0)
< 0 に注意する。ここで、大数の法則により
1
n
n∑i=1
logf(Xi|θ)f(Xi|θ0)
→ Eθ0
[log
f(X|θ)f(X|θ0)
]in prob.
32
とできる。ここで、Jensenの不等式を (log x)′′ < 0に注意して用いると、
Eθ0
[log
f(X|θ)f(X|θ0)
]< logEθ0
[ f(X|θ)f(X|θ0)
]= log
∫ ∞−∞
f(x|θ)f(x|θ0)
f(x|θ0) dx = 0
となり主張を得る。 □
命題 5.3 (Jensenの不等式) f(x)が下に凸であれば
E[f(X)] ≥ f(E[X]).
特に、f ′′(x) > 0であれば等号成立は X が定数のときに限る。
証明: µ = E[X] とする。f(x) は下に凸であるから c ∈ R があって f(x) ≥ c(x − µ) + f(µ) とできる。
従って、E[f(X)] ≥ E[c(X − µ) + f(µ)] = c(E[X]− µ) + f(µ) = f(µ).
f ′′(x) > 0であれば f(x) > f ′(µ)(x− µ) + µ, x = µ, となるので、P (X = µ) = 1でなければ等号は成立し
ない。 □
定理 5.4 (最尤推定量の一致性) (A1)–(A3) の条件のもと、∀x に対し Θ ∋ θ 7→ f(x|θ) は微分可能で θ につ
いての偏微分を f ′(x|θ)と表すとき、尤度方程式
l′(θ|x) =n∑i=1
f ′(xi|θ)f(xi|θ)
= 0 (5.1)
は根 θn = θn(x1, . . . , xn)で θn(X1, . . . , Xn)→ θ0 in prob. となるものをもつ。特に、最尤方程式 (5.1)がた
だ一つの解 θn をもつ、すなわち、θn が最尤推定量であれば、θn は真のパラメータ θ0 に確率収束する。すな
わち、最尤推定量は一致推定量である。
証明: a > 0を [θ0 − a, θ0 + a] ⊂ Θとなるようにとり、
Sn = {x ; ln(θ0|x) > ln(θ0 − a|x) and ln(θ0|x) > ln(θ0 + a|x) }
とする。定理 5.2より Pθ0(X ∈ Sn)→ 1. ここで、∀x ∈ Sn に対して、θn ∈ (θ0 − a, θ0 + a)をそこで ln(θ|x)が極大となるようにとれる。ここで l′(θn|x) = 0に注意する。よって、
Pθ0(|θn − θ0| < a) ≥ Pθ0(X ∈ Sn)→ 1.
ここで、上記の θn が aに依存していることに注意する。これを解決するためには、ln(θ|x)が極大となる θで
|θ − θ0|を最小とするものを選び θn とすればよい。(x 7→ θn の可測性が問題となるが、その証明は省略する。
Lehmann-Casella [LC] p.504を参照のこと。) □
漸近有効性のために次の命題を準備する。
命題 5.5 確率変数 Xn が定数 a > 0に確率収束し、確率変数 Zn が Z に法則収束するとき、Zn/Xn は Z/a
に法則収束する。
証明: cが FZ(z)が連続点 ⇐⇒ c/aが FZ/a(z/a)が連続点であり、FZ(z) = FZ/a(z/a)より、FZ の任意の
連続点 cに対して、 limn→∞
P (Zn/Xn ≤ c/a) = P (Z ≤ c)を示せばよい。ε > 0を任意にとっておく。δ > 0を
|x− c| < δ =⇒ |FZ(x)− FZ(c)| <ε
4
33
ととる。FZ の不連続点は高々可算個なので、c1, c2 を FZ の連続点で c− δ < c1 < c < c2 < c+ δ となるよう
にできる。次に limn→∞
FZn(ci) = FZ(ci) (i = 1, 2) と Xn → a in prob. より、N ∈ Nを
n ≥ N =⇒ |FZn(ci)− FZ(ci)| <ε
4(i = 1, 2), P (An
c) <ε
4
となるように選ぶ。ただし、An ={|Xn − a| <
a
|c|+ 1δ′}, δ′ = min{c− c1, c2 − c, 1}とした。ここで、An
上で c1 < c− δ′ < c
aXn < c+ δ′ < c2 かつ Xn > 0となることに注意する。よって、
P
(ZnXn≤ c
a
)≤ P
({ZnXn≤ c
a
}∩An
)+ P (An
c) = P({Zn ≤
c
aXn
}∩An
)+ P (An
c)
≤ P (Zn ≤ c+ δ′) + P (Anc) ≤ P (Zn ≤ c2) +
ε
4≤ P (Z ≤ c2) +
ε
2,
< P (Z ≤ c) + ε,
P
(ZnXn≤ c
a
)≥ P
({ZnXn≤ c
a
}∩An
)≥ P
({Zn ≤
c
aXn
}∩An
)≥ P ({Zn ≤ c− δ′} ∩An)
≥ P (Zn ≤ c− δ′)− P (Anc) ≥ P (Zn ≤ c1)−ε
4≥ P (Z ≤ c1)−
ε
2> P (Z ≤ c)− ε.
以上より n ≥ N ならば |P (Zn/Xn ≤ c/a)− FZ(c)| < εとなるので主張を得る。 □
命題 5.6 Xn → 0 in prob.であり、確率変数の族 {Yn}は tight, すなわち、∀ε > 0に対してM > 0があっ
て infnP (|Yn| ≤M) ≥ 1− εとすると、XnYn → 0 in prob. となる。
証明: ε > 0を任意にとっておく。条件より、M > 0を P (|Yn| > M) <ε
2とでき、さらに ∀δ > 0に対して
N ∈ Nを n ≥ N ならば P(|Xn| ≥
δ
M
)<ε
2とできる。このとき、n ≥ N なら
P (|XnYn| ≥ δ) = P (|XnYn| ≥ δ, |Yn| ≤M) + P (|Yn| > M)
≤ P (M |Xn| ≥ δ, |Yn| ≤M) + P (|Yn| > M)
≤ P(|Xn| ≥
δ
M
)+ P (|Yn| > M) < ε. □
以下、(A1)–(A3)に加え次を仮定する。X は Xi と同じ分布をもつ確率変数とする。
(A4) f(x|θ)は θ について C3 級.
(A5) Eθ0
[ ∂∂θ
log f(X|θ0)]= 0, I(θ0) := Eθ0
[− ∂2
∂θ2log f(X|θ0)
]= Eθ0
[{ ∂
∂θlog f(X|θ0)
}2]∈ (0,∞).
(A6) ある c > 0とM(x)が存在して、x ∈ A (cf . (A2))と |θ − θ0| < cに対して∣∣∣ ∂3∂θ3
log f(x|θ)∣∣∣ ≤M(x)
であり Eθ0 [M(X)] <∞.
もし∫f(x|θ) dxが θ について 2回微分可能なら (A5)は成立すると期待できる。実際、
∫ ∞−∞
f(x|θ) dx = 1
より ∫ ∞−∞
∂
∂θf(x|θ) dx =
∂
∂θ
(∫ ∞−∞
f(x|θ) dx)
= 0, 同様に∫ ∞−∞
∂2
∂θ2f(x|θ) dx = 0. (5.2)
よって、∂
∂θlog f(x|θ) = 1
f(x|θ)∂f
∂θ(x|θ), ∂2
∂θ2log f(x|θ) = 1
f(x|θ)∂2f
∂θ2(x|θ)−
( 1
f(x|θ)∂f
∂θ(x|θ)
)2より、
Eθ0
[ ∂∂θ
log f(x|θ0)]= Eθ0
[ 1
f(X|θ0)∂f
∂θ(X|θ0)
]=
∫ ∞−∞
∂
∂θf(x|θ0) dx = 0, (5.3)
Eθ0
[ ∂2∂θ2
log f(x|θ0)]= Eθ0
[ 1
f(X|θ0)∂2
∂θ2f(X|θ0)
]− Eθ0
[{ 1
f(X|θ0)∂f
∂θ(X|θ0)
}2](5.4)
34
=
∫ ∞−∞
∂2
∂θ2f(x|θ0) dx− Eθ0
[{ 1
f(X|θ0)∂f
∂θ(X|θ0)
}2]= −Eθ0
[{ ∂
∂θlog f(X|θ0)
}2]= −I(θ0)
となる。ここで、I(θ)は Fisher情報量 (前期の Cramer-Raoの定理のところで学んだ)であった。
定理 5.7 (最尤推定量の漸近有効性) X1, X2, . . . は i.i.d. で (A1)–(A6) の条件を満たすとする。尤度方程式
(5.1)の根 θn = θn(x1, . . . , xn)は
√n(θn(X)− θ0)→ N(0, 1/I(θ0)) in law (5.5)
を満たす。特に、最尤方程式 (5.1)がただ一つの解 θn をもつ、すなわち、θn が最尤推定量であれば、(5.5)が
成り立つ。ただし θn(X) = θn(X1, · · · , Xn)とした。
証明: x = (x1, · · · , xn)を固定し、l′(θn|x)について θ0 のまわりで Taylorの定理を適用すると
l′(θn|x) = l′(θ0|x) + (θn − θ0)l′′(θ0|x) +1
2(θn − θ0)2l′′′(θ∗n|x)
となる。ここで、θ∗n = θ0 + h(θn − θ0), 0 < h < 1, である。仮定より (左辺)= 0なので、
√n(θn(X)− θ0) =
1√nl′(θ0|X)
− 1n l′′(θ0|X)− 1
2n (θn(X)− θ0)l′′′(θ∗n|X).
ここで、(A5)より中心極限定理から
1√nl′(θ0|X) =
1√n
n∑i=1
∂
∂θlog f(Xi|θ0)→ N(0, I(θ0)) in law.
また、(A5)と大数の法則より、
1
nl′′(θ0|X) =
1
n
n∑i=1
∂2
∂θ2log f(Xi|θ0)→ Eθ0
[ ∂2∂θ2
log f(X|θ0)]= −I(θ0) in prob.
次に定理 5.4より θn → θ0 in prob. なので、∀ε > 0に対し N を n ≥ N ならば Pθ0(|θn − θ0| ≥ c) < ε/2と
とると、|θ∗n − θ0| ≤ |θn − θ0|より (A6)と大数の法則より
∣∣∣ 1nl′′′(θ∗n|X)1{|θn−θ0|<c}
∣∣∣ ≤ 1
n
n∑i=1
∣∣∣ ∂3∂θ3
log f(Xi|θ∗n)∣∣∣1{|θn−θ0|<c} ≤ 1
n
n∑i=1
M(Xi)→ Eθ0 [M(X)] in prob.
よって、M0 = Eθ0 [M(X)]として、十分大きい nに対して
Pθ0
(∣∣∣ 1nl′′′(θ∗n|X)
∣∣∣ ≥M0 + 1)≤ Pθ0(|θn − θ0| ≥ c) + Pθ0
(∣∣∣ 1nl′′′(θ∗n|X)1{|θn−θ0|<c}
∣∣∣ ≥M0 + 1)< ε
とできるので、命題 5.5, 5.6より、Z ∼ N(0, I(θ0))として
√n(θn(X)− θ0)→
1
I(θ0)Z in law
となるが、1
I(θ0)Z ∼ N
(0,
1
I(θ0)
)より結論を得る。 □
この定理は最尤推定量 θn(X)が漸近的に平均 θ0, 分散 1nI(θ0)
−1 の正規分布に従うことを示している。
例 5.1 X1, X2, . . .が i.i.d.で各 Xi は指数分布 Ex(λ), λ > 0, に従うとき、λの最尤推定量 λとその漸近分布
を求めよ。
35
解: f(x|λ) = λe−λx とする。尤度方程式
l′n(λ|x) =∂
∂λ
n∑i=1
log f(xi|λ) =n
λ− (x1 + · · ·+ xn) = 0
より、最尤推定量は λn =1
Xn
, ただし Xn =1
n(X1 + · · ·+Xn). Fisher情報量は
I(λ) = −E[ ∂2∂λ2
log f(X|λ)]=
1
λ2.
よって、λは漸近的に N(λ,λ2
n
)に従う。 □
例 5.2 ある保険契約の各契約者の 1 年間のクレーム件数の分布は、互いに独立に同一の幾何分布 Ge(p),
0 < p < 1, に従うものとする。ある年度のクレーム件数を調べたところ、全契約者 10,000人のクレーム件数
は 500件であった。このとき、この母集団の母数 pの 95% 信頼区間を、最尤推定量の漸近有効性を利用して
求めよ。*4
解: f(x|p) = (1− p)xp, x = 0, 1, 2, . . ., とする。尤度方程式
l′n(p|x) =∂
∂p
n∑i=1
log f(xi|p) = −x1 + · · ·+ xn
1− p+n
p= 0
より、最尤推定量は pn =1
1 +Xn
, ただし Xn =1
n(X1 + · · ·+Xn). Fisher情報量は
I(p) = −E[ ∂2∂p2
log f(X|p)]= E
[ X
(1− p)2+
1
p2
]=
1
(1− p)21− pp
+1
p2=
1
(1− p)p2.
よって、p =1
1 +Xは漸近的に N
(p,
(1− p)p2
n
)に従う。ここで、漸近分散
(1− p)p2
nは
(1− p)p2
nに十分
近いと考え、実現値を代入して p = 1/(1 + 50010000 ) = 0.95238 · · · で
p± 1.960
√(1− p)p2
n= 0.95238± 1.960
√0.04762 · 0.952382
10000= 0.95238± 0.00407 =
{0.956450.94831
となるので、95%信頼区間は 0.9483 ≤ p ≤ 0.9565を得る。 □
問題 5.1 a > 0 を定数、θ > 0 を未知母数とする。X1, X2, . . . が i.i.d. で Xi の密度関数を f(x|θ) =
aθaxa−1e−(θx)a
, x > 0, とする。このとき、最尤推定量 θn と、θn の漸近分布を求めよ。
未知母数が k 次元のとき、未知母数 θ = (θ1, . . . , θk)に対し、Fisher情報量 I(θ)を次で定義する。
I(θ) =
I11(θ) · · · I1k(θ)...
...Ik1(θ) · · · Ikk(θ)
,
Iij(θ) = Eθ
[{ ∂
∂θilog f(X|θ)
}{ ∂
∂θjlog f(X|θ)
}]= −Eθ
[ ∂2
∂θi∂θjlog f(X|θ)
], i, j = 1, . . . , k.
上式の変形は (5.4)と同様に示せる。
(A1)–(A6)と同様の条件のもと定理 5.8が成立する。証明は省略する。ただし、(A5)では I(θ0) ∈ (0,∞)
のかわりに、I(θ)が正定値行列であることを仮定する。
*4 [IK] 岩沢宏和 黒田耕嗣 著 損害保険数理 (アクチュアリー数学シリーズ 4), 日本評論社, 2015
36
定理 5.8 (最尤推定量の分散 (損保数理 p.2–2)) X1, X2, . . . は i.i.d. で上記の条件を満たすとする。尤度方
程式
∂
∂θjln(θ|x) =
n∑i=1
∂
∂θjlog f(xi|θ) = 0, j = 1, . . . , k (5.6)
の根 θn = θn(x1, . . . , xn)は
√n(θn(X)− θ0)→ N(0, I(θ0)
−1) in law (5.7)
を満たす。特に、最尤方程式 (5.6)がただ一つの解 θn をもつ、すなわち、θn が最尤推定量であれば、(5.7)が
成り立つ。
この定理は最尤推定量 θn(X)が漸近的に平均ベクトル θ0, 共分散行列1
nI(θ0)
−1 の k 次元正規分布に従う
ことを示している。
例 5.3 α, β > 0 を未知母数とする。X1, X2, . . . が i.i.d. で各 Xi は分布関数が F (x|α, β) = 1 −( β
x+ β
)α(x ≥ 0) で与えられる確率変数であるとする。このとき、α, β の最尤推定量 αn, βn の漸近分布を求めよ。
解: 密度関数は f(x|α, β) = αβα
(x+ β)α+1. よって
∂
∂αlog f(x|α, β) = 1
α+ log β − log(x+ β),
∂
∂βlog f(x|α, β) = α
β− α+ 1
x+ β,
∂2
∂α2log f(x|α, β) = − 1
α2,
∂2
∂α∂βlog f(x|α, β) = 1
β− 1
x+ β,
∂2
∂β2log f(x|α, β) = − α
β2+
α+ 1
(x+ β)2.
ここで、
E[(X + β)−k] =
∫ ∞0
αβα
(x+ β)α+1+kdx =
α
(α+ k)βk, k > −α,
より、Fisher情報量は I(α, β) =
1
α2− 1
(α+ 1)β
− 1
(α+ 1)β
α
(α+ 2)β2
. よって、I(α, β)−1 を計算して、(αn, βn)は
漸近的に N
((α
β
),1
n
α2(α+ 1)2 α(α+ 1)(α+ 2)β
α(α+ 1)(α+ 2)β(α+ 1)2(α+ 2)β2
α
) に従う。 □
問題 5.2 X1, X2, . . . が i.i.d. で µ, σ2 を未知母数とする対数正規分布に従う、すなわち、その密度関数が
f(x|µ, σ2) =1√
2πσ2xexp{− (log x− µ)2
2σ2
}, x > 0, であるとする。このとき、µ, σ2 の最尤推定量 µn, σ2
n
の漸近分布を求めよ。
定理 5.9 (最尤推定量の関数の分散 (損保数理 p.2–4)) k 次元確率変数の列 Xn = (X1n, · · · , Xkn) が漸近的
に平均 θ = (θ1, · · · , θk), 共分散行列1
nΣの正規分布に従い、θ と Σはいずれも nに依存しないものと仮定す
る。gを全微分可能な k変数関数であるとし、Gn = g(X1n, · · · , Xkn)とすると、Gn は漸近的に平均 g(θ), 分
散1
n(∂g)′Σ(∂g)の正規分布に従う。ここで、(∂g) =
(∂g
∂θ1(θ) · · · ∂g
∂θk(θ)
)′(縦ベクトル)である。
この定理は、定理 5.8 の仮定の下、最尤推定量 θn の関数 g(θn) は漸近的に、平均 g(θ0), 分散1
n(∂g)′I(θ0)
−1(∂g)の正規分布に従うことを示している。
37
略証: 簡単のため k = 2 として示す。まず a, b が定数で、(X,Y ) が N((µ1
µ2
),
(σ11 σ12
σ12 σ22
))に従うとき、
aX + bY は平均 aµ1 + bµ2 で次の分散の正規分布に従うことに注意する:
V (aX + bY ) = a2V (X) + 2abCov(X,Y ) + b2V (Y ) =(a b
)(σ11 σ12σ12 σ22
)(ab
). (5.8)
g(x1, x2)は全微分可能なので、多変数の平均値の定理より
Gn − g(θ1, θ2) =∂g
∂θ1(X∗1n, X
∗2n)(X1n − θ1) +
∂g
∂θ2(X∗1n, X
∗2n)(X2n − θ2) (5.9)
となる。ここで、(X∗1n, X∗2n) = (θ1, θ2) + h(X1n − θ1, X2n − θ2), 0 < h < 1, である。仮定より
X1n → θ1 in prob. かつ X2n → θ2 in prob.
が証明できるので、
∂g
∂θ1(X∗1n, X
∗2n)→
∂g
∂θ1(θ1, θ2) in prob. かつ
∂g
∂θ2(X∗1n, X
∗2n)→
∂g
∂θ1(θ1, θ2) in prob.
を得る。さらに、仮定より (√n(X1n − θ1),
√n(X2n − θ2))は平均 (0, 0), 共分散行列 Σの 2次元正規分布に
法則収束するので、その極限分布に従う確率変数を (Y1, Y2)とすると、(5.9)の両辺に√nを掛け n → ∞と
すると、√n(Gn − g(θ1, θ2)) →
∂g
∂θ1(θ1, θ2)Y1 +
∂g
∂θ2(θ1, θ2)Y2 in law
となるので、(5.8)に注意して主張を得る。 □
例 5.4 例 5.2で、この母集団の母平均 µ =1− ppの 95% 信頼区間を、最尤推定量 pn の漸近有効性を利用し
て求めよ。
解: 例 5.2より、最尤推定量 pn =1
1 +Xで、√n(pn − p)→ N(0, (1− p)p2)であった。母平均を pn を用い
て表すと g(p) =1− ppとするとき、g(pn)と表せるから、g(pn)は漸近的に平均 g(p), 分散
1
ng′(p)(1− p)p2g′(p) = 1
n
(− 1
p2
)2(1− p)p2 =
1− pnp2
の正規分布に従うと考えられる。従って、95% 信頼区間は
g(pn)± u(0.025)
√1− pnnp2n
= x± u(0.025)√x(1 + x)
n= 0.05± 1.960
√0.05(1 + 0.05)
10000= 0.05± 0.00449
となるので、95%信頼区間は 0.0455 ≤ µ ≤ 0.0545を得る。 □
問題 5.3 問題 5.2の対数正規分布に従う母集団で、標本数を nとするとき、(µn, σ2n)の漸近有効性を用いて、
次の問いに答えよ。
(1) σ =
√σ2n の漸近分布を求めよ。
(2) 母集団分布の平均の漸近分布を求めよ。ヒント: 対数正規分布での平均は eµn+12 σ
2n であった。
5.2 極値問題
ここでは、損保数理の重要なテーマである極値問題を扱う。基本的な参考図書として次がある。
[R] S.I. Resnick: Extreme Values, Regular Variation and Point Processes, Springer, 1987
[TS] 高橋 倫也, 志村 隆彰: 極値統計学 (ISMシリーズ:進化する統計数理), 近代科学社, 2016
38
定理 5.10 (型収束定理 (損保数理 p.10–55)) 確率変数列 {Xn}が、定数 an > 0, αn > 0, bn, βn ∈ Rと a.s.
には定数でない確率変数 X,Y があって、
Xn − bnan
→ X in law,Xn − βnαn
→ Y in law, (5.10)
であるための必要十分条件は、
limn→∞
anαn
= a > 0, limn→∞
bn − βnαn
= b (5.11)
となることである。このとき、a, bは一意的に定まり、Y と aX + bは同じ分布に従う。
ここで、確率変数X が a.s.には定数ではないというのは、P (X = c) = 1となる定数 cが存在しないという
意味である。また、ある定数 a > 0と b ∈ Rがあって、aX + bと Y が同じ分布に従うとき、すなわち、
FY (x) = P (aX + b ≤ x) = FX
(x− ba
)となるとき、X と Y の分布は同じ型に属するという。
証明は難しいので省略する (cf . [R] pp.7–8)。直感的には、(5.10)より
Y ≍ Xn − βnαn
=anαn
Xn − bnan
+bn − βnαn
≍ anαn
X +bn − βnαn
であるが、X,Y がともに a.s.には定数ではないことより (5.11)と必要十分となると理解できる。
定理 5.11 (Fisher-Tippett(損保数理 p.10–5)) X1, X2, . . .を i.i.d.とし、その最大値Mn = max{X1, · · · , Xn}を考える。このとき、a.s.には定数でない確率変数 Y と定数 cn > 0, dn ∈ Rが存在して、
Mn − dncn
→ Y in law (5.12)
とすれば、Y の分布は次のいずれかの分布関数のもつ分布と同じ型に属する。
(Frechet分布) Φα(x) =
{exp (−x−α) , x > 0,
0, x ≤ 0.
(Weibull分布) Ψα(x) =
{exp (−(−x)α) , x ≤ 0,
1, x > 0
(Gumbel分布) Λ(x) = exp (−e−x).ここで α > 0は定数である。
−4 −2 0 2 4
0.0
0.2
0.4
0.6
0.8
定理の Φα(x), Ψα(x), Λ(x)で与えられ分布
を総称して極値分布という。その分布の密度関
数のグラフは右図のようになる。(α = 2 とし
た。) どのグラフがどの分布に対応するかは、
密度関数が正となる範囲を考えればすぐわかる
であろう。
定理 5.11の証明: Xi, Y の分布関数を
F (x) = P (Xi ≤ x), H(x) = P (Y ≤ x)とすると、Mn の分布関数は (F (x))n であり、
(5.12) より H(x) の任意の連続点 x ∈ R に対
して
P(Mn − dn
cn≤ x
)=(F (cnx+ dn)
)n → H(x), n→∞. (5.13)
39
(5.13)は nを [nt]に置き換えても成り立つから、
P(M[nt] − d[nt]
c[nt]≤ x
)=(F (c[nt]x+ d[nt])
)[nt] → H(x), n→∞. (5.14)
ここで [a]は aの整数部分を表すものとする。同様に、M[nt] の分布関数は (F (x))[nt] であるから、
P(M[nt] − dn
cn≤ x
)=(F (cnx+ dn)
)[nt]={(F (cnx+ dn)
)n}[nt]/n
→ H(x)t (5.15)
を得る。(5.14), (5.15) より、{H(x)}t を分布関数にもつ、したがって a.s. に定数ではない確率変数 Yt が
あって、M[nt] − d[nt]
c[nt]→ Y in law,
M[nt] − dncn
→ Yt in law,
となる。よって、型収束定理より、
limn→∞
cnc[nt]
= γ(t) > 0, limn→∞
dn − d[nt]c[nt]
= δ(t) (5.16)
となる γ(t), δ(t) (t > 0) が存在し、γ(t)Yt + δ(t)と Y が同じ分布に従う:{H(x)
}t= P (Yt ≤ x) = P (Y ≤ γ(t)x+ δ(t)) = H(γ(t)x+ δ(t)). (5.17)
ここで、tのかわりに stを代入して {H(x)
}st= H(γ(st)x+ δ(st)).
一方、 {H(x)
}st={(H(x)
)s}t={H(γ(s)x+ δ(s))
}t= H
(γ(t){γ(s)x+ δ(s)}+ δ(t)
).
上の二式が任意の xに対して成立するので、後で述べる補題 5.12より次の等式を得る。
γ(st) = γ(s)γ(t), δ(st) = γ(t)δ(s) + δ(t), s, t > 0 (5.18)
次に、γ(e) = eθ, θ ∈ R, とすると、(5.18)より s, t > 0と n ∈ Nに対して γ(sn) = γ(s)n であるから、t = sn
として γ(t1/n) = γ(t)1/n となり、したがって、∀r ∈ Qに対して γ(er) = γ(e)r = (eθ)r = (er)θ を得る。こ
こで、(5.16)より γ(t)は可測関数なので次を得る*5。
γ(t) = tθ, t > 0. (5.19)
(a) θ = 0のとき: γ(t) ≡ 1なので (5.18)より δ(st) = δ(s) + δ(t). 従って、上と同様に δ(e) = cとすると、
∀r ∈ Qに対して δ(er) = rδ(e) = rcとなるので、δ(t)の可測性より
δ(t) = δ(elog t) = c log t, t > 0,
となり、(5.17)より {H(x)
}t= H(x+ c log t), t > 0. (5.20)
もし c = 0なら、H(x) = 0, 1となり、これは Y が a.s.に定数であることを意味するので、c = 0となる。次
に、0 ≤ H(x) ≤ 1より{H(x)
}tは tについて非増加なので、c < 0を得る。
次にある x0 ∈ Rがあって H(x0) = 1とすると、(5.20)より H(x0 + c log t) = 1, ∀t > 0. これは ∀x ∈ R
に対し H(x) = 1 を意味しており、H(x) が分布関数であることに矛盾する。したがって ∀x ∈ R に対し
*5 証明はかなり面倒なので省略する。詳しくは、例えば次の [BGT] pp.4–5 を見よ。[BGT] N.H. Bingham, C.M. Goldie and J.L. Teugels: Regular Variation, Cambridge University Press, 1987.
40
H(x) < 1. 同様に、x0 ∈ Rがあって H(x0) = 0とすると、∀x ∈ Rに対し H(x) = 0となり矛盾する。した
がって ∀x ∈ Rに対し H(x) > 0. (5.20)に x = 0を代入して、
{H(0)}t = H(c log t).
ここで、c′ = −1/c > 0, d ∈ Rを exp{−e−d} = H(0) ∈ (0, 1)とし x = c log tを代入すると、t > 0のとき
−∞ < x <∞でH(x) =
{H(0)
}ex/c
= exp{−e−de−c′x} = Λ(c′x+ d).
(b) θ > 0のとき: (5.18)の第 2式より、
γ(t)δ(s) + δ(t) = γ(s)δ(t) + δ(s)
より、t = 1, s = 1のときδ(s)
1− γ(s)=
δ(t)
1− γ(t),
すなわち、δ(s)
1− γ(s)は定数関数なのでそれを cとする。よって、t = 1のとき
δ(t) = δ(s)1− γ(t)1− γ(s)
= c(1− tθ)
であり、
{H(x)}t = H(tθx+ c(1− tθ)) = H(tθ(x− c) + c) (5.21)
特に x = cとして H(c)t = H(c)となるので、H(c) = 0または 1. もし H(c) = 0なら、ある x > cがあって
H(x) > 0となるので、右連続性より
{H(x)}t = H(tθ(x− c) + c)→ H(c) = 0 as t→ +0.
一方、a = H(x) > 0なら at → 1 as t→ +0となるので矛盾する。これよりH(x) ≡ 0となり、これはH(x)
が分布関数であることに矛盾する。したがって、H(c) = 1となる。
次に、0 < H(c− 1) < 1を示す。H(c− 1) = 1とすると、(5.21)よりH(−tθ+ c) = H(tθ(c− 1− c)+ c) ={H(c−1)}t = 1, t > 0, となるので、H ≡ 1となりH が分布関数であることに矛盾する。また、H(c−1) = 0
とすると、全く同様にH ≡ 0となり矛盾する。したがって、p > 0があってH(c− 1) = exp{−pα}とできる。ただし、α = 1/θ とした。このとき、x < cに対して、x− c = −tθ とすると t = (c− x)α に注意して (5.21)
を用い
H(x) = H(x− c+ c) = H(−tθ + c) = {H(c− 1)}t = exp{−(p(c− x))α} = Ψα(p(x− c)).
(c) θ < 0のとき: (b)のときと同様に cを定めると、(5.21)が成り立つ。よって、H(c) = 0, 1となるが、も
し H(c) = 1ならある x < cがあって H(x) < 1となるが、
H(c− 0) = limt→∞
H(tθ(x− c) + c) = limt→∞{H(x)}t = 0
となり、Y = c a.s. となり矛盾する。よって H(c) = 0.
さらに、(b)とまったく同様に 0 < H(1 + c) < 1が示されるので、p > 0があってH(1 + c) = exp{−p−α}とできる。ただし、α = −1/θとした。このとき、x > cに対して、x− c = tθ とすると t = (x− c)−α に注意して (5.21)を用い
H(x) = H(x− c+ c) = H(tθ + c) = {H(1 + c)}t = exp{−(p(x− c))−α} = Φα(p(x− c)). □
41
補題 5.12 Xを a.s.には定数ではない確率変数としF (x)をその分布関数とする。このとき、定数 a > 0, c > 0,
b, d ∈ Rが、∀x ∈ Rに対して F (ax+ b) = F (cx+ d)を満たせば、a = c, b = dとなる。
証明: a ≥ cとして一般性を失わない。G(x) = F (cx+ d)とすると、
F (ax+ b) = F(c(acx+
b− dc
)+ d)= G
(acx+
b− dc
)となることに注意し、α =
a
c≥ 1, β =
b− dcと表す。仮定より、G(x) = G(αx+ β)であるが、これを繰り返
し代入して
G(x) = G(αx+ β) = G(α(αx+ β) + β) = G(α2x+ (α+ 1)β)
= G(α{α2x+ (α+ 1)β}+ β) = · · · = G(αnx+ (αn−1 + · · ·+ 1)β) (5.22)
を得る。よって、α > 1であれば
G(x) = G(αnx+
αn − 1
α− 1β)= G
(αn(x+
β
α− 1
)− β
α− 1
)であるが、F (x)は a.s.には定数ではない確率変数の分布関数だから、G(x)もそうなるので、ある 0 < p < 1
と x0 = −β
α− 1があって、G(x0) = pとできる。このとき、n→∞とすると
p = G(x0) = G(αn(x0 +
β
α− 1
)− β
α− 1
)→
{1, x0 +
βα−1 > 0,
0, x0 +βα−1 < 0,
となるので、いずれにせよ矛盾する。よって α = 1. ここで、もし β = 0なら、(5.22)で n→∞とすると
p = G(x0) = G(x0 + nβ)→{
1, β > 0,0, β < 0,
となり矛盾する。したがって、β = 0となる。 □
例 5.5 X1, X2, . . . が i.i.d.でその分布と定数 cn > 0, dn ∈ Rが次の (1)–(3)で与えられるとき、定理 5.11の
極限分布 Y の分布関数を求め、どの分布と同じ型か調べよ (Φα,Ψα,Λで表せ)。
(1) Parate分布に従う: fX(x) =
{βx−β−1, x ≥ 1,
0, x < 1,で、cn = n1/β , dn = 0. ただし β > 0は定数。
(2) ベータ分布 Beta(1, b)に従う: fX(x) =
{b(1− x)b−1, 0 < x < 1,
0, その他,で、cn = n−1/b, dn = 1.
(3) 指数分布 Ex(λ)に従い、cn =1
λ, dn =
1
λlog n.
解: Φα(x),Ψα(x),Λ(x)はすべて −∞ < x < ∞ で連続なので (5.13)より limn→∞
{F (cnx + dn)}n をそれで表せばよい。
(1) Xi の分布関数は F (x) =
{1− x−β , x ≥ 1,
0, x < 1.
x ≤ 0のとき F (cn + dn) = 0. x > 0のとき nが十分大なら n1/βx > 1となることに注意すると、
{F (cnx+ dn)}n = {1− (n1/βx)−β}n =(1− x−β
n
)n→ e−x
−β
= Φβ(x).
以上より、極限分布の分布関数は Φβ(x)となる。(2), (3)は演習問題とする。 □
問題 5.4 例 5.5 (2), (3)を解け。
42
例題 5.6 独立に平均 9の指数分布に従うX1, X2, . . . , Xn の最大値をMn とする。このとき、M81 が 81以上
である確率を、以下の 2通りの方法で求めよ。
(1) 最大値の分布関数から直接計算する方法。
(2) 最大値を正規化した分布を極値分布により近似 (cf . 例 5.5 (3))して計算する方法。
解: (1) x ≥ 0に対して、FXi(x) = 1− e−x/9 であるから、
P (M81 ≥ 81) = 1− P (M81 < 81) = 1− (1− e−81/9)81 = 0.00994700 · · · .
(2) 例 5.5 (3)より P
(Mn − 9 log n
9≤ x
)→ e−e
−x
, n→∞, となるから、
P (M81 ≥ 81) = 1− P(M81 − 9 log 81
9≤ 81− 9 log 81
9
)≒ 1− exp{−e−9+log 81}
= 1− exp{−81e−9} = 0.00994639 · · · . □
問題 5.5 X1, X2, . . . , Xn は i.i.d.で例 5.5 (1)で β =9
7とした Parate分布に従うとし、その最大値をMn と
する。このとき、M128 が 128以上である確率を、以下の 2通りの方法で求めよ。(ヒント: 128 = 27.)
(1) 最大値の分布関数から直接計算する方法。
(2) 最大値を正規化した分布を極値分布により近似 (cf . 例 5.5 (1))して計算する方法。ただし e−1 = 0.36788
として計算せよ。
定義 5.2 (一般化極値分布) σ > 0, µ ∈ Rのとき、分布関数が
Hξ,µ,σ(x) = exp{−(1 + ξ
x− µσ
)−1/ξ}, 1 + ξ
x− µσ
> 0, (5.23)
である分布のことを一般化極値分布あるいは GEV 分布 (generalized extreme value distribution) という。
ただし、ξ = 0のときは
H0,µ,σ(x) = limξ→0
Hξ,µ,σ(x) = exp{−e−
x−µσ
}, −∞ < x <∞, (5.24)
とする。ξ, µ, σ はそれぞれ形状パラメータ、位置パラメータ。尺度パラメータと呼ばれ、µ = 0, σ = 1とした
とき、特に Hξ(x)と記す。
GEVと Φα(x),Ψα(x),Λ(x)は次のように関係づけられる。
• Frechet分布 Φα(x) = H1/α
(x− 1
1/α
), ξ = 1/α > 0
• Gumbel分布 Λ(x) = H0(x), ξ = 0
• Weibull分布 Ψα(x) = H−1/α
(x+ 1
1/α
), ξ = −1/α < 0
命題 5.13 α > 0, Z ∼ Ex(1)とすると、次が成立する。
(1) X = − logZ は Gumbel分布 Λに従い、MX(t) = E[etX ] = Γ(1− t). 特に、E[X] = γ, V (X) = π2/6.
ここで、γ = limn→∞
(1 +
1
2+ · · ·+ 1
n− log n
)= 0.577215664 · · · (Euler の定数) を表す。
(2) X = Z−1/α は Frechet分布 Φα に従い、E[Xn] = Γ(1− n
α
), 1− n
α> 0.
(3) X = −Z1/α はWeibull分布 Ψα に従い、E[Xn] = (−1)nΓ(1 +
n
α
), 1 +
n
α> 0.
証明: (1)のみ示し、残りは演習問題とする。X の分布とその積率母関数については
P (X ≤ x) = P (Z ≥ e−x) =∫ ∞e−x
e−z dz = e−e−x
= Λ(x),
43
MX(t) = E[etX ] = E[Z−t] =
∫ ∞0
z−te−z dz = Γ(1− t)
より従う。次に、X のキミュラント母関数 ψX(t) = logMX(t)に関して、ψ′X(0) = E[X], ψ′′X(0) = V (X)と
なることに注意する。これを計算するため、ガンマ関数のWeierstrassの公式
1
Γ(x)= xeγx
∞∏p=1
(1 +
x
p
)e−x/p
を用いる。(証明は特殊関数の教科書を調べてください。) これより、
ψX(t) = log Γ(1− t) = − log(1− t)− γ(1− t)−∞∑p=1
{log
(1 +
1− tp
)− 1− t
p
},
ψ′X(t) =1
1− t+ γ −
∞∑p=1
{ −1/p1 + (1− t)/p
+1
p
}=
1
1− t+ γ −
∞∑p=1
{1p− 1
p+ 1− t
},
ψ′′X(t) =1
(1− t)2+
∞∑p=1
1
(p+ 1− t)2.
44
以上より、
E[X] = ψ′X(0) = 1 + γ −∞∑p=1
{1p− 1
p+ 1
}= 1 + γ − 1 = γ,
V (X) = ψ′′X(0) = 1 +
∞∑p=1
1
(p+ 1)2=
∞∑k=1
1
k2=π2
6. □
問題 5.6 命題 5.13 (2), (3)を示せ。
定義 5.3 (最大吸引域) 定理 5.11で Xi を F、(5.12)を満たす Y の分布関数を H とするとき、定数 cn > 0,
dn ∈ Rが存在して、ある ξ に対して H(x) = Hξ(x)できる。このとき、F は一般化極値分布 Hξ の最大吸引
域 (maximum domain of attraction)に属するといい、F ∈MDA(Hξ)と表す。ただし、Hξ のかわりに Φα,
Λ, Ψα を用いることも多い。
多くの確率分布が、上記のどれかの最大吸引域に属することが知られてる。そのための必要十分条件や、も
う少しチェックしやすい十分条件も知られている。ただし、連続型確率分布であっても一般化極値分布の最大
吸引域に属さない例が知られているし、また、Poisson分布や幾何分布も一般化極値分布の最大吸引域に属さ
ないことが知られている。詳しくは [R], [TS]を参照のこと。
次の表は三つのタイプの極限分布に関するまとめである ([TS], pp.35–36より)。ただし、F (x) = 1− F (x)は末尾分布関数、緩慢関数の定義は定理 5.16の次にある。
• Frechet分布 の最大吸引域
Frechet分布 Φα(x) = exp(−x−α), x > 0, α > 0
F ∈MDA(Φα) ωF =∞, F (x) = x−αl(x), l(x)は緩慢関数
基準化定数 cn = F←(1− 1/n), dn = 0
極限への収束 Mn/cn → Φα in law
例: 裾の厚い (fat-tailed)分布
Cauchy分布 f(x) = 1/[π(1 + x2)], x ∈ R, cn = n/π, α = 1
Parate分布 F (x) ∼ Kx−α, K,α > 0, cn = (Kn)1/α
対数ガンマ分布 f(x) = αβ
Γ(β) (log x)β−1x−α−1, x > 1, α, β > 0,
cn = [(log n)β−1n/Γ(β)]1/α
• Weibull分布 の最大吸引域
Weibull分布 Ψα(x) = exp(−(−x)α), x < 0, α > 0
F ∈MDA(Ψα) ωF <∞, F (ωF − x−1) = x−αl(x), l(x)は緩慢関数
基準化定数 cn = ωF − F←(1− 1/n), dn = ωF
極限への収束 (Mn − ωF )/cn → Ψα in law
例: 有限な上限をもつ分布
一様分布 f(x) = 1, 0 < x < 1, cn = 1/n, dn = 1, α = 1
ωF でベキ F (x) = K(ωF − x)α, ωF −K−1/α ≤ x ≤ ωF , K,α > 0,
cn = (Kn)−1/α, dn = ωF
ベータ分布 f(x) = 1B(a,b)x
a−1(1− x)b−1, 0 < x < 1, a, b > 0,
Beta(a, b) cn =[Γ(a)Γ(b+1)nΓ(a+b)
]1/b, dn = 1, α = b
45
• Gumbel分布 の最大吸引域
Gumbel分布 Λ(x) = exp(−e−x), x ∈ R
F ∈MDA(Λ) ωF ≤ ∞, F (x) = c(x) exp(−∫ xαF
g(t)a(t) dt
), αF < x < ωF ,
ただし、x→ ωF のとき c(x)→ c > 0, g(x)→ 1, a′(x)→ 0.
基準化定数 dn = F←(1− 1/n), cn = a(dn)
極限への収束 (Mn − dn)/cn → Λ in law
例: 統計学でよく用いられる分布
ガンマ分布 f(x) = 1Γ(λ)x
λ−1e−x, x > 0, λ > 0
Γ(λ, 1) cn = 1, dn = log n+ (λ− 1) log log n− log Γ(λ)
正規分布 f(x) = 1√2π
exp(−x2
2 ), x ∈ R
N(0, 1) cn = (2 log n)−1/2, dn =√2 log n− [log log n+ log(4π)]/(2
√2 log n)
対数正規分布 f(x) = 1√2πσx
exp(− (log x−µ)22σ2 ), x > 0, µ ∈ R, σ > 0
LN(µ, σ2) cn = σ(2 log n)−1/2dn, dn = exp{µ+ σ
[√2 log n− log logn+log(4π)
2√2 logn
]}逆ガウス分布 f(x) =
(λ
2πx3
)exp[−λ(x−µ)
2
2µ2x
], x > 0, µ, λ > 0,
cn = 2µ2/λ, dn = µ2{2 log n− 3 log log n+ log λ2
4πµ2 + 2λµ
}/λ
ωF <∞で F (x) = K exp(− aωF−x
), x < ωF , a,K > 0,
指数挙動 cn = a/[log(Kn)]2, dn = ωF − a/ log(Kn)
では、データが与えられたときそれがどの一般化極値分布の最大吸引域かが問題となるが、それには通常、
最尤法が用いられる。具体的には、いくつかのブロックに分けて考えるブロック最大値モデル、ブロックの最
大値だけでなく、ある閾値を超過したデータを対象とした閾値超過モデルがある。詳しくは損保数理の教科書、
さらにはそこに出展されている次の書籍の第 7章を参照されたい。
[MFE] A.J. McNeil, R. Frey, P. Embrechts: 定量的リスク管理 -基礎概念と数理技法-, 共立出版, 2008.
5.3 安定分布
極値問題では、独立同分布をもつ確率変数の最大値Mn について、その「型」(cf . 定理 5.10)を調べた。こ
の節ではその和 Sn = X1 + · · ·+Xn について考える。
中心極限定理によると、もし σ2 := V (X1) <∞であれば、
1√n(Sn − nE[X1])→ N(0, σ2) in law
となる。一方、各 Xk が Cauchy分布に従う、すなわちその密度関数が f(x) =1
π(1 + x2)であるとき、その
特性関数は ϕXk(t) = e−|t| であったから、Yn =
1
nSn に対し、
ϕYn(t) = E[eitYn ] =
n∏k=1
E[eitnXk ] =
n∏k=1
ϕXk(t/n) = (e−|t/n|)n = e−|t|
となり、Yn も Cauchy分布に従う、特に、「Yn → Cauchy分布 in law 」を得る。では、一般に Sn について
どんな「型」がありうるのだろうか。一般に次が知られている。
定義 5.4 (安定分布) 確率変数 X が安定分布に従うとは、任意の k ∈ N と X と同分布に従う独立な
X1, . . . , Xk に対して、ある定数 ak > 0, bk があって、X1 + · · ·+Xk ∼ akX + bk となるときにいう。
46
定理 5.14 確率変数 Z が、ある定数 an > 0, bn に対してSn − bnan
→ Z in law となるための必要十分条件は、
Z が安定分布に従うことである。
定理 5.15 確率変数 Z が安定分布に従うとき、Z は正規分布であるか、指数 α ∈ (0, 2)と定数 β,m1,m2 が
あって、その特性関数 ϕZ(t)は次で与えられる。
ϕZ(t) = exp
{iβt+m1
∫ ∞0
(eitx − 1− itx
1 + x2
) dx
x1+α+m2
∫ 0
−∞
(eitx − 1− itx
1 + x2
) dx
|x|1+α
}.(5.25)
さらに、指数 αが 0 < α < 2であれば、c ∈ R, d > 0と |θ| ≤ 1を用いて
ϕZ(t) =
exp{ict− d|t|α
(1 + iθ
t
|t|tan
πα
2
)}, α = 1,
exp{ict− d|t|
(1 + iθ
t
|t|2
αlog |t|
)}, α = 1,
と変形できる。
証明は、例えば次の第 9章にある。
[B] L. Breiman, Probability, Addison-Wesley, 1968. (Classics in applied mathematics, 7, Society for
Industrial and Applied Mathematics, 1992. Reprint版)
ここでは、次の第 3章 7節に述べられている次の定理を紹介するにとどめる。
[D] Durrett, R.: Probablity Theory and Examples, 4th ed., Cambridge University Press, 2010.
定理 5.16 X1, X2, . . .は i.i.d.で、次を満たすとする。
(i) limx→∞
P (X1 > x)
P (|X1| > x)= θ ∈ [0, 1],
(ii) 0 < α < 2と緩慢関数 L(x)があって、P (|X1| > x) = x−αL(x).
このとき、Sn = X1 + · · ·+Xn に対して、
an = inf{x;P (|X1| > x) ≤ 1/n}, bn = nE[X11{|X1|≤an}]
とすると、a.s.には定数でない確率変数 Z があって、Sn − bnan
→ Z in law となる。ここで、Z の特性関数は
(5.32)あるいは (5.25)でm1 = αθ, m2 = α(1− θ), β は (5.33)として与えられる。
ここで、L(x)が緩慢関数であるとは、∀t > 0に対して limx→∞
L(tx)
L(x)= 1となるときにいう。例えば、log x
は緩慢関数である。一方、xα は α = 0のとき緩慢関数ではない。
補題 5.17 X が定理 5.16 (ii) を X = X1 として満たすとすると、次が成立する。
(1) ∀δ > 0に対して、定数 Cδ > 0と xδ > 0が存在して、y ≥ xδ に対して P (|X| > y) ≥ Cδy−(α+δ).特に、E[X2] =∞.
(2) 定数 C > 0と x0 > 0が存在して、x ≥ x0 ならば∫ x
0
yP (|X| > y) dy ≤ Cx2P (|X| > x).
証明: (1) 仮定より、ある xδ > 0が存在して、x ≥ xδ ならば
P (|X| > 2x)
P (|X| > x)= 2−α
L(2x)
L(x)≥ 2−(α+δ)
とできる。よって、n ≥ 1に対して、P (|X| > 2nxδ) ≥ P (|X| > xδ)2−(α+δ)n となる。ここで、∀y ≥ xδ に対
して nを 2n−1 ≤ y/xδ < 2n ととれば、
P (|X| > y) ≥ P (|X| > 2nxδ) ≥ P (|X| > xδ)2−(α+δ)n ≥ Cδy−(α+δ)
47
を得る。ただし Cδ = (xδ/2)α+δP (|X| > xδ)とした。次に、δ > 0を α+ δ < 2と選べば、
E[X2] =
∫ ∞0
2yP (|X| > y) dy ≥∫ ∞xδ
2yCδy−(α+δ) dy =
[2Cδ
2− (α+ δ)y2−(α+δ)
]∞xδ
=∞.
(2) I(x) =
∫ x
0
yP (|X| > y) dy とおくと、t > s > 0に対して
I(tx)− I(sx) =∫ tx
sx
yP (|X| > y) dy = x2∫ t
s
uP (|X| > ux) du. (5.26)
ε > 0を任意とするとある T > 0があって、u ≥ T ならば P (|X| > ux)
P (|X| > u)≥ (1− ε)x−α とできるので、
I(tx)− I(sx) ≥ (1− ε)x2−α∫ t
s
uP (|X| > u) du = (1− ε)x2−α{I(t)− I(s)}, t > s > T.
ここで、 limt→∞
I(t) =
∫ ∞0
yP (|X| > y) dy = E[|X|]/2 =∞より、両辺を I(t)で割って t→∞とすると、
lim inft→∞
I(tx)
I(t)≥ (1− ε)x2−α, 左辺は ε > 0によらないので lim inf
t→∞
I(tx)
I(t)≥ x2−α. (5.27)
(5.26)で s = 1として両辺を I(x)で割って変形すると、
I(tx)
I(x)= 1 +
x2P (|X| > x)
I(x)
∫ t
1
uP (|X| > ux)
P (|X| > x)du ≤ 1 +
x2P (|X| > x)
I(x)
t2 − 1
2.
ここで、x→∞として、(5.27)より t > 1のとき、
t2−α ≤ lim infx→∞
I(tx)
I(x)= 1 +
t2
2lim infx→∞
x2P (|X| > x)
I(x), すなわち、lim inf
x→∞
x2P (|X| > x)
I(x)≥ 2(t2−α − 1)
t2.
よって、x0 > 0があって、x ≥ x0 ならばx2P (|X| > x)
I(x)≥ t2−α − 1
t2となり、証明は完了した。 □
定理 5.16の証明: まず、 limn→∞
nP (|X1| > an) = 1に注意する。これは、定義より nP (|X1| > an) ≤ 1であ
り、また、∀ε > 0に対して条件 (ii)と an の取り方より
(1 + 2ε)−α = limn→∞
P (|X1| > (1+2ε)an1+ε )
P (|X1| > an1+ε )
≤ lim infn→∞
P (|X1| > an)
1/n
となることから従う。よって、(ii)より、x > 0に対して、
nP (|X1| > xan) =P (|X1| > xan)
P (|X1| > an)nP (|X1| > an)→ x−α, (5.28)
さらに、(i)と組み合わせて、x > 0に対して、
nP (X1 > xan) =P (X1 > xan)
P (|X1| > xan)nP (|X1| > xan)→ θx−α, (5.29)
同様に、x < 0に対して nP (X1 < xan)→ (1− θ)|x|−α を得る。ε >を任意とし、In(ε) = {k ∈ N; 1 ≤ k ≤ n, |Xk| > εan}とし、更に、
µ1(ε) = E[X11{εan<|X1|≤an}], µ2(ε) = E[X11{|X1|≤εan}],
Sn(ε) =∑
k∈In(ε)
Xk, Tn(ε) = Sn − bn − (Sn(ε)− nµ1(ε)) =
n∑k=1
{Xk1{|Xk|≤εan} − µ2(ε)},
48
とおく。
Tn(ε)について、Xεk = Xk1{|Xk|≤εan} とかくと、
E[(Tn(ε)/an)2] = V (Tn(ε)/an) = nV (Xε
1/an) ≤ nE[(Xε1/an)
2] = n
∫ ∞0
2xP (|Xε1/an| > x) dx
≤ 2n
∫ ε
0
xP (|X1| > xan) dx =2n
a2n
∫ εan
0
yP (|X1| > y) dy
であるが、εan →∞に注意して補題 5.17 (2)と (5.28)を用いて、次を得る。
E[(Tn(ε)/an)2] ≤ 2n
a2nC(εan)
2P (|X1| > εan)→ 2Cε2−α. (5.30)
Sn(ε)について、#In(ε) = k の条件のもと、Sn(ε)/an の分布は F εn(x) = P (X1/an ≤ x||X1|/an > ε)を
共通の分布関数とする k 個の独立な確率変数の和として表されることに注意する。また、#In(ε)は二項分布
B(n, pεn), pεn = P (|X1| > εan), に従うので、F εn の特性関数を ψεn(t)とすると、
E[eitSn(ε)/an ] = E[E[eitSn(ε)/an |#Iεn]
]=
n∑k=0
(ψεn(t))k
(n
k
)(pεn)
k(1− pεn)n−k = {1 + (ψεn(t)− 1)pεn}n.
ここで、
F εn(x) =
1− P (X1 > xan)
P (|X1| > εan)→ 1− θ(x/ε)−α, x > ε,
P (X1 ≤ xan)P (|X1| > εan)
→ θ(|x|/ε)−α, x < −ε,
より、
ψεn(t)→ ψε(t) := εα∫ ∞ε
eitxθαdx
xα+1+ εα
∫ −ε−∞
eitx(1− θ)α dx
|x|α+1
となることと、npεn → ε−α および∫ ∞ε
αx−α−1 dx = ε−α に注意して、
E[eitSn(ε)/an ]→ eε−α(ψε(t)−1) = exp
(∫ ∞ε
(eitx − 1)θαdx
xα+1+
∫ −ε−∞
(eitx − 1)(1− θ)α dx
|x|α+1
)を示すことができる。また、(5.29)より、
nµ1(ε)/an = nE[(X1/an)1{εan<|X1|≤an}]→∫ 1
ε
xθαdx
xα+1+
∫ −ε−1
x(1− θ)α dx
|x|α+1
を得る。従って、E[exp{it(Sn(ε)− nµ1(ε))/an}]→
exp
{∫ ∞1
(eitx − 1)θαdx
xα+1+
∫ 1
ε
(eitx − 1− itx)θα dx
xα+1
+
∫ −1−∞
(eitx − 1)(1− θ)α dx
|x|α+1+
∫ −ε−1
(eitx − 1− itx)(1− θ)α dx
|x|α+1
}(5.31)
となる。ここで、eitx − 1− itx ∼ −t2x2/2 (x→ 0) より、0 < α < 2のとき、∫ 1
0
(eitx − 1− itx) dx
xα+1,
∫ 0
−1(eitx − 1− itx) dx
|x|α+1
はともに収束する。よって、ε→ 0として (5.31)は
exp
{iβt+
∫ ∞0
(eitx − 1− itx
1 + x2
)θα
dx
x1+α+
∫ 0
−∞
(eitx − 1− itx
1 + x2
)(1− θ)α dx
|x|1+α
}, (5.32)
ただし、 β = (2θ − 1)α{∫ ∞
1
x
1 + x2dx
x1+α+
∫ 1
0
( x
1 + x2− x) dx
x1+α
}(5.33)
に収束する。
定理 5.16の証明を完了させるために次の補題を準備する。
49
補題 5.18 各 ε > 0で hn(ε) → h(ε), gn(ε) → g(ε) (n → ∞) でかつ h(ε) → h(0), g(ε) → g(0) (ε → 0) で
あれば、{εn}を hn(εn)→ h(0)かつ gn(εn)→ g(0)となるように選ぶことができる。
証明: N1 < N2 < . . .を、n ≥ Nm ならば |hn(1/m)− h(1/m)| ≤ 1/mかつ |gn(1/m)− g(1/m)| ≤ 1/mと
なるように選べる。ここで、n < N1 のとき εn = 1, Nm ≤ n < Nm+1 のときは εn = 1/mと選べば、
|hn(εn)− h(0)| ≤ |hn(1/m)− h(1/m)|+ |h(1/m)− h(0)| ≤ 1/m+ |h(1/m)− h(0)|
で、n→∞のときm→∞となるので、hn(εn)→ h(0)となる。gn(εn)→ g(0)も同様に示せる。 □
定理 5.16 の証明の続き: 補題 5.18 を hn(ε) = E[(Tn(ε)/an)2], gn(ε) = E[exp{it(Sn(ε) − nµ1(ε))/an}]
として用いると、{εn} が存在して、(5.30) より hn(εn) → 0, また gn(ε) → (5.32) とできる。これより、
Tn(εn)/an → 0 in L2, 特に、Tn(εn)/an → 0 in prob. であり、一方、(5.32)は t = 0で連続なので確率統計
学 II定理 6.30 (2)より (5.32)はある確率変数 Y の特性関数であり、(Sn(εn)− nµ1(εn))/an → Y in law と
なることが従う。以上より、
Sn − bnan
=Tn(εn)
an+Sn(εn)− nµ1(εn)
an→ Y in law
を得る。(一般に、Xn → 0 in prob.かつ Yn → Y in lawであれば、Xn + Yn → Y in law となることを用い
た。) □
定義 5.5 (吸引域) Xk の分布関数を F とする。ある定数 an > 0, bn が存在して、Sn − bnan
が指数 α ∈ (0, 2)
の安定分布 Z に法則収束するとき、F (x) あるいは Xk の分布は指数 α の安定分布 Z の吸引域に属すると
いう。
定理 5.16は Xi の分布が (i), (ii)を満たせば指数 αの吸引域に属することを意味している。
安定分布の理論はより広い族である無限分解可能分布や対応する確率過程である Levy過程 (独立増分をも
ち時間的一様 Xt+s −Xt ∼ Xs となる、例えば Brown運動や Poisson過程、複合 Poisson過程など)と深く
関連しており、詳しく調べられている。その教科書としては例えば次がある。
[S] 佐藤 健一: 加法過程, 紀伊國屋書店, 1999. (より新しい内容を含む著者による英語版もある。)
50
