制限付き最尤法 REML; Restricted Maximum Likelihood

Presentation on theme: "制限付き最尤法 REML; Restricted Maximum Likelihood"— Presentation transcript:

1 制限付き最尤法 REML; Restricted Maximum Likelihood
増田　豊

2 正規分布と最尤法 Normal Distribution and Maximum Likelihood

3 正規分布 【正規分布の形状を示す】 1. 横軸･･･とり得る値、縦軸･･･出現頻度
2. この分布のグラフは f(y) のようになる。これを確率密度関数という。 　この部分の面積が1となるように考えられている。 　高さ（頻度）を求めるにはこの関数に値を代入する。 3. 分布の形状は m と s によって決定できる。正規分布を N(m,s)と表す。 4. 左右対称である。 5. 身長・体重など、生物の観察値は正規分布に従うことが知られている。

4 正規分布と母数 分布の特性を表す値を母数という 正規分布の母数は m と s2 である ある変数 Y が正規分布に従うとき
Y ~ N( m, s2 ) と表記する 【正規分布の母数を示す】

5 m と分布のイメージ m は頻度が最大になる位置（平均）を表す

6 s2 と分布のイメージ s2 は分布の広がり（分散）を表す

7 母集団と標本 【母集団と標本の関係】 正規母集団からの標本も正規分布をなす 標本（データ）から母数を推定することができる → 統計的推測

8 統計的推測の例 ある中学校の生徒30人の身長を測定した 母集団（この学校の全生徒）の身長の平均と分散を推定したい。
身長測定値は正規分布に従うとする。

9 頻度分布と基礎統計量 【基礎統計量を求める】 この結果から、母集団の平均は165、分散は10程度であると考えられる。
では、このデータに実際に正規分布を当てはめてみる。

10 正規分布の当てはめ N ( 165, 15 ) N ( 165, 10 ) N ( 165, 5 ) 【実際に正規分布を当てはめる】
このグラフから当てはまりの程度を判断することは難しい。 そこで、仮定した分布にデータがどのくらいよく当てはまっているかを示す目安が必要である。 その目安となる数値が尤度である。

11 尤度 （l） データが、特定の分布にどの程度当てはまっているかを示す数値 尤度が大きいほど、よく当てはまっている
尤度の自然対数がよく用いられる （対数尤度 ... log l）

12 尤度の計算 尤度と対数尤度は以下のように計算する

13 尤度の計算

14 対数尤度の計算

15 例データの尤度

16 尤度が最大になる点 ただし s2 は不偏分散、 S2 は標本分散

17 最尤法と最尤推定量 データを固定して母数を動かし、尤度を最大にする母数を探す方法を最尤法という。 このときの値を最尤推定量という。
この場合は、以下のようになる。 および

18 最尤法と尤度関数 尤度 l の値はmとs2の値によって変動する 尤度を母数 m と s2 の関数と見なすとき、これを尤度関数という。
最尤法は、尤度関数を最大にする母数を求める方法である。 結局は関数の最大値を求める問題である。

19 尤度関数の最大値を探す 母数を適当に動かし、そのつど尤度を計算する方法が、直接検索法である
直接検索法は、求める母数が多くなると、ちょっと（かなり）面倒である 二分法やシンプレックス法などがある 直接検索法を用いることは derivative-free（DF）であると言われる

20 最尤法の問題点 記録数が少ないときに偏りがある 尤度関数を最大にする点を探すのが困難 平均と分散の両方を推定することによるもの
この例の場合では、不偏分散ではない 尤度関数を最大にする点を探すのが困難 この例の場合は、簡単に求められる 一般には、面倒な数値計算が必要である REMLでもこの問題は解消されない

21 制限付き最尤法の導出 Derivation of REML

22 モデル 以下のモデルを仮定する y は多変量正規分布に従う 期待値と分散を以下のように仮定する y = Xb + e = 1m + e
y:観察値, m:集団平均, e:残差（変量効果） y は多変量正規分布に従う 期待値と分散を以下のように仮定する E(y) = Xb = 1m Var(y) = Var(e) = Is2

23 対数尤度関数

24 最尤推定量 (ML) log L を最小にする m と s2 の値 log L を m と s2 で偏微分し、0とする
（関数の最小値を求める） この場合は、以下のようになる および

25 偏微分した式に注目する s2 の式に m が含まれている いま m が未知なので、 で代用する の値による変動はどの程度かを調べる
推定量の中に推定量を含んでいる この結果、　　の値によって　　が変化する 　　の値による変動はどの程度かを調べる

26 式の変形

27 条件付き期待値 標本　　と分散　　を仮定したときの期待値 分散は未知なので、適当な値を仮定する

28 条件付き期待値の値 標本平均　　の分布は次のようになる 分散の定義から、条件付き期待値の値は

29 推定式 適当な値　　を設定し、　　を計算する 求めた　　を　　に代入し、計算を繰り返す 値が動かなくなったとき「収束」したという

30 収束値 この収束値をREML推定量という この場合は、不偏分散である 分散についてはMLよりもREMLのほうが好ましい推定量であるようだ

31 EMアルゴリズム 「期待値をとる→代入」を繰り返す方法 計算式が単純である 解が必ず母数空間内に収束する 収束が非常に遅い
（収束値に近づくほど、速度が遅くなる）

32 REML推定量を導く方法 ML推定量を修正する 尤度関数そのものを修正する 今回の例のように式変形して期待値をとる
制限付き尤度関数を最大にする母数を求める

33 尤度関数の比較

34 REML推定値の求めかた Computational Methods of REML Estimates

35 主なアプローチ derivative-freeな方法 EMアルゴリズム 数値計算による方法 直接検索により試行錯誤的に求める
尤度関数を1回だけ微分して期待値をとる 数値計算による方法 収束が非常に速い 尤度関数を2回微分した式が必要になる

36 考え方 尤度関数を最大とする値を得たい = 1次導関数が 0 となる点を探す
式で考えれば簡単だが、コンピュータにとっては困難（具体的な数値しか扱えない） 徐々に解に近づくような計算を繰り返す これを「数値計算」という

37 ニュートン・ラフソン法 g(x) この点で接する接線を求める y=a x+b 傾き a: g’(x0) (x0, g(x0))
切片 b: g(x0)-g’(x0)x0 接線 y= g’(x0) x+g(x0)-g’(x0)x0 (x0, g(x0)) x x0

38 ニュートン・ラフソン法 g(x) x軸と交わる点 x1 を求める 接線 y= g’(x0) x+g(x0)-g’(x0)x0
y=0 を代入して x1 を求める x x1 x0

39 ニュートン・ラフソン法 g(x) これを繰り返して近似値を得る r 回めの反復結果 x x2 x1 x0

40 ニュートン・ラフソン法 r 回めの反復結果 g(x) に相当するのは log Lr の1次導関数である

41 1次および2次導関数 反復式は以下のようになる

42 ニュートン・ラフソン法の特徴 初期値によっては収束しない 収束値付近では速度が上昇する 負の分散が推定されることがある
収束値から離れた値を与えると発散しやすい 収束値付近では速度が上昇する 負の分散が推定されることがある 2次導関数（分母）の式が複雑になる 計算量が多くなりやすい

43 スコアリング法 2次導関数のかわりに、期待値を用いる 反復式は以下のようになる

44 スコアリング法の特徴 非常に収束しやすい 母数効果のみを含むモデルの場合は、反復式が非常に簡単になる
Fisherのスコアリング法とも呼ばれる 尤度関数の2次導関数の負の期待値が、Fisher情報量と呼ばれているため

45 AIアルゴリズム ニュートン・ラフソン法とスコアリング法で用いた分母を「平均」する これを分母として反復式をつくる

46 AIアルゴリズム AI = Average Information である ニュートン・ラフソン法よりは収束しやすい
スコアリング法よりは発散しやすい 混合モデルの場合は、他の方法よりも計算量が少ない

47 アニマルモデルとREML Animal Model and REML

48 モデル 以下のモデルを仮定する y は多変量正規分布に従う 期待値と分散を以下のように仮定する y = Xb + Za +e
y:観察値, b:母数効果, a:育種価, e:残差 X, Z: 計画行列 y は多変量正規分布に従う 期待値と分散を以下のように仮定する E(y) = Xb Var(y) = V = ZGZ’ + R = ZAZ’sa2 + Ise2

49 混合モデル方程式 または

50 Derivative-free (DF) REML
Smith と Graser(1987), Graser ら(1987)が提案した方法 MeyerのDFREML,BoldmanらのMTDFREMLに実装され、爆発的に普及した 尤度の計算を繰り返し、直接検索によって最大となる点を探す

51 Derivative-free (DF) REML
尤度の計算に log|C| と y’Py が必要である C を y’R-1y に掃き出すと y’Py が得られる ピボットの対数の総和が log|C| である コレスキー分解を用いる方法もある

52 Derivative-free (DF) REML
母数の数が少ないとき、計算が高速で、個体数の増加に対応できる 母数の数が増える（多形質を扱う）とき、検索精度・速度が大きく低下する 初期値を変えて、何度も実行する必要がある（コールドスタート）

53 EM REML Harville(1977)が示し、Henderson(1984)が命名した方法 1次導関数を 0 とした式から導く
左辺の期待値が右辺の値に等しい

54 EM REML 以下の反復から求める

55 EM REML 計算式が簡単である 必ず母数空間内に収束する 収束値から遠いときに収束が早く、収束値に近づくほど速度が遅くなる
逆行列の一部 Caa が必要であり、計算量が多くなる

56 AI REML JohnsonとThompson(1995)が最初に応用し、Gilmour(1995)らによって定式化した方法
数値計算に用いられていた NR / MSC 法に基づくが、計算量を大幅に軽減している 現在、主流になっている方法で、パッケージも多い

57 AI REML 数値計算に用いる反復式は以下の通り d は、尤度関数の1次導関数である ニュートン・ラフソン法 Hに2次導関数
Hにこれらの平均を用いるのが AI REML

58 AI REML 少ない反復回数で収束する NR / MSC より計算量が少ない DF / EMよりは計算量が多い
初期値によっては収束しないことも多い 多くのパッケージでは、負の値が推定されたときにEMアルゴリズムを組み合わせる

59 現実問題 Reality

60 実際のプログラムは 数々の高速化技法と組み合わせている コンピュータの物理的な限界を考慮 疎行列(sparse matrix)演算など
メモリ容量の制限 プログラムの内部構造に起因する限界 丸め誤差や有効桁数などへの対策 システムに起因する（本質的でない）問題

61 実際の推定は 推定にかかる時間よりも、データの加工や検証にかける時間が多い（かもしれない）
大きなデータが扱えないとき、細分化した一部（サブセット）を用いることがある まずモデルが適切かどうかを検討すべき 推定値を解釈しなければならない

62 REML以外の推定法 Bayesian approach Method R ベイズ統計学に基づく方法 ギブズサンプリングを応用している
計算量が少なく大きなデータも扱える Method R 計算量がきわめて少なく、巨大データを扱える 統計学的に疑問があり、値が偏ることがある