枝の長さと側枝の数の関係をモデル化する

確率分布の最尤推定の試み

Updated on 12 April 2002

あとは，このパラメータの値の決め方です． 当然，もっともよく観測データと合うように決めたい． となると最尤法の出番でしょうか． あいにく，わたしは今まで最尤法を使った経験はありません．でも， 「困ったときはネットに聞け」とも，「頼みの網（ネット）」とも言います． 久保拓弥さんの 尤度を使った推定法あれこれ のページに掲載されている 「区画」をとらない「率」推定（最初の一歩篇） (PDF, 320KB) では，生態学者に身近な例を題材に，最尤法の初歩が要領よく説明されています． まずはこれを勉強してみることにしました．

最尤法によるパラメータ推定の基本は，データの推定値と観測データとの合致の 度合を表す尤度が最大になるパラメータを探すことです．そのためには， 尤度をパラメータで微分したものがゼロになるところを探せばよいのですが， 解析的には求まらないのがふつうのようです． 久保さんのテキストでは，多次元のNewton-Raphson法を使ってパラメータの 最尤推定をする方法が紹介されています．

そこまで勉強するのは大変そうなので，とりあえずは愚直にパラメータを 動かしながら尤度が最大になるところを探すことにしました． aと bふたつのパラメータのうちまずは a を固定しておき，b を動かして 尤度が最大になるところを探す，次に b をそこに固定しながら a を 動かして尤度が最大になるところを探す，という作業を繰り返します． a,bがじゅうぶんに収束してきたら探索を終了します．

尤度は以下のような手順で計算してみました．

枝の長さを x として，その時に側枝の数が n 本以下である確率 pは p = a / (x^b + a) となります．もし観測データの側枝の数が n以下だったら，確率 pの 事象が起きたことになるし，側枝の数が nより大きかったら 確率 1- pの事象が起きたことになります． この確率を全部の測定データについて求めて掛け合わせたものが尤度です． ただし，1より小さい数をどんどん掛けていくとすごく小さい値になって いくので，まず対数をとってから足し算して，対数尤度を求めるのが ふつうのようです．

対数尤度の値は，当然 aやbの値によって変化します．尤度が最大になる a と bの値を探すことで，観測されたようなデータが得られる可能性が 最大になるような，すなわちもっとも観測データとよく合うパラメータを 見つけることができます．

尤度が最大になるところを探すアルゴリズムを工夫しつつ， Ruby でプログラムを書いてみました（プログラムはきれいにしてから載せる予定）． 下の表が推定したパラメータです．

最大で６本まで側枝を持つものがあるので，P(x, 5)まで求めてあります． ６本の側枝がでる確率は1 - P(x, 5) です． 推定したパラメータを使って，P(x, n)の形を描いてみました．

たとえば長さ40センチの枝の場合，1本も側枝がない確率が10%程度， 1本だけ側枝がある確率が20%，2本の側枝がある確率が一番高くて50%程度， というように読み取れます．

85個のデータを説明するのに， P(x, 0) からP(x, 5) までそれぞれ２個ずつのパラメータ，全部で12個のパラメータ を動員してますから，これで観測データが再現できなかったらがっかりです． 測定された枝長データを使っておそるおそる側枝の本数を計算してみました．

P(x, n)で与えられるのは確率だけですから，実際に何本の側枝を持つかは 乱数で決めます．４回の試行を行った結果が下の図です．

観測データ と，こわいほど似ているような気がするのですが，いかがでしょうか．

今後の課題は，

• パラメータ推定を，ちゃんとした方法でやる．
• P(x, n)の関数形としてもっとよいものがないか検討する．

などです．