5

R でこの事後分布を計算しようとしています。問題は、dbern(p_i, y_i) < 1 の束の積である分子が小さすぎることです。（私のnは約1500です）。したがって、R は 0 を吐き出し、すべての \theta の事後値も 0 です。

明確にするために、各 y_i には独自の p_i があり、これらの p_i を合わせて、n 個の y に対して n 個の要素のベクトルを作成します。各 theta には、p_i の独自の n 要素ベクトルがあります。

ここに画像の説明を入力

再現可能な例 (分子の)

p <- sample(seq(0.001,0.999,by=0.01), 1500, replace=T)
y <- sample(c(0,1), 1500, replace=T)
dbern(y, p) # 1500-element vector, each element < 1
prod(dbern(y, p)) # produces 0
exp(sum(log(dbern(y, p)))) # produces 0

編集 (コンテキスト): ベイジアン変化点分析を行っています (jstor.org/stable/25791783 - Western and Kleykamp 2004)。論文の連続 y とは異なり、私の y は 2 進数なので、Albert と Chib (1993) のデータ拡張法を使用しています。この方法では、y の尤度はベルヌーイ (p = cdf-normal(x'B)) です。

では、p はどのようにシータに依存するのでしょうか? シータが変化点だからです。x の 1 つはタイム ダミーです。たとえば、theta=10 の場合、10 日目以降のすべての観測ではタイム ダミー = 1、10 日より前のすべての観測では 0 になります。

したがって、p は x に依存し、x は theta に依存します。したがって、p は theta に依存します。

Gibbs サンプリングのシータの完全な条件であるため、上記の数量が必要です。

4

3 に答える 3

5

このような精度の問題に取り組む 1 つの方法は、対数空間で作業することです。ただし、それは分母から対数和積を導入するため、一般的には苦痛になる可能性があります。

最適化のために事後分布を計算している場合は、分母を完全に削除できる可能性があることに注意してください。正規化して を見つける必要はありませんargmax。

于 2013-04-21T02:44:37.677 に答える
3

さて、私はあなたの例を実行しまし0た。0

p <- sample(seq(0,1,by=0.01), 1500, replace=T)
y <- sample(c(0,1), replace=T)
x <- dbern(y, p)
any(x == 0)
## [1] TRUE
于 2013-04-21T06:37:07.903 に答える
1

Cross Validated でもこの質問をしたところ、glen_bから以下の (テスト済みの) 解決策が得られました。

これは、あらゆる種類のモデルの尤度の計算に共通の問題です。一般的に行われる種類のことは、ログで作業することと、値をより妥当な範囲にする共通の倍率を使用することです。

この場合、次のことをお勧めします。

ステップ 1: かなり「典型的な」θ、θ0 を選択します。一般項の分子と分母の両方の式を θ=θ0 の分子で割ると、アンダーフローの可能性がはるかに低くなります。

ステップ 2: 対数スケールで作業します。これは、分子が対数の差の合計の exp であり、分母が対数の差の合計の exp の合計であることを意味します。

注意: p のいずれかが 0 または 1 の場合は、それらを個別に取り出して、それらの用語のログを取らないでください。そのまま評価しやすい！

分子の通常の項は、より適度なサイズになる傾向があるため、多くの場合、分子と分母はどちらも比較的妥当です。

分母にサイズの範囲がある場合は、大きいものを追加する前に小さいものを合計します。

1 つまたはいくつかの用語が大きく支配している場合は、それらの用語を比較的正確に計算することに注意を向ける必要があります。

于 2013-04-22T06:15:23.057 に答える