モデル選択のベイジアン設定で、反復手順を使用して限界尤度を推定しようとしています。詳細に興味がある場合は、たとえばthisまたはthis paper を参照してください。私はかなり長い間数値の問題と戦っているので、ここで試してみて、運を試してみようと思いました。以下に続く内容を、できるだけ短く、再現可能で、一般的なものにするように努めました。
設定
目標は、フォームの反復プロセスを実行することです
私のモデルの限界尤度yの推定値を取得します(この反復プロセスは、かなり速くyのある値に収束します)。式に表示される他のすべての変数は、すでに計算されているスカラーまたはパラメーターの可能性であるため、固定値です。各反復で変化する唯一の項はyの値です。基本的に、長さ L のベクトル p_l と q_l と、長さ M のベクトルp_mとq_mがあります。
問題
上記の式で使用されている尤度値は、残念ながら対数尤度 (計算してベクトルに保存) ではなく、実際の尤度です。これが私が抱えている数値の問題の原因です。マグニチュードが大きい負の対数尤度を指数化するという古い問題は、おそらくご存知でしょう。私の対数の可能性は非常に負であるため、exp(対数の可能性) はほとんど常に 0 になります。回答が既にオンラインになっている同様の質問がいくつかありますが、それらの解決策のどれも私にとってはうまくいきませんでした。
私が試したこと
exp(log(x)) = x という事実を利用して分数を展開すると、上記の式を次のように書き換えることができると思います。
ここで、Cは任意の定数です。同じ考えに従う証明については、この論文の付録を参照してください。合計の項を扱いやすくするCの適切な値を見つけることができれば、問題は解決されます。ただし、私の場合、 p_l、q_l、p_m 、およびq_mは大きさが大きく異なるため、Cのどの値を減算または加算しようとしても、アンダーフローまたはオーバーフローが再び発生します。したがって、今のところ、ここから先に進む方法が実際にはよくわかりません。コメントやヒントをいただければ幸いです。
コードとデータ
サンプル データ (つまり、対数尤度を含むベクトル) をここで見つけることができます。反復プロセスのコードは次のとおりです。
L <- 1000
M <- 5000
y <- numeric(length=100)
eval_downer <- numeric(length=L)
eval_upper <- numeric(length=M)
y[1] <- 0.5
for(t in 2:100){
for(m in 1:M){
up.m <- q_m[m]
down.m <- (L * q_m[m]) + (M * p_m[m] / y[t-1])
eval_downer[m] <- up.m / down.m
}
for(l in 1:L){
up.l <- p_l[l]
down.l <- (L * q_l[l]) + (M * p_l[l] / y[t-1])
eval_upper[l] <- up.l / down.l
}
upper <- mean(eval_upper)
downer <- mean(eval_downer)
y[t] <- upper / downer
print(t)
}
ありがとうございました！