Schritte, um eine posteriore Verteilung herauszufinden, wenn es einfach genug sein könnte, eine analytische Form zu haben?

12

Dies wurde auch bei Computational Science gefragt .

Ich versuche, eine Bayes'sche Schätzung einiger Koeffizienten für eine Autoregression mit 11 Datenproben zu berechnen: wobei ist Gauß mit Mittelwert 0 und Varianz Die vorherige Verteilung auf dem Vektor ist Gauß mit Mittelwert und einer diagonalen Kovarianzmatrix mit diagonale Einträge gleich .

Yi=μ+αYi1+ϵi
ϵiσe2(μ,α)t(0,0)σp2

Basierend auf der Autoregressionsformel bedeutet dies, dass die Verteilung der Datenpunkte ( ) normal mit dem Mittelwert und der Varianz . Somit wäre die Dichte für alle Datenpunkte gemeinsam (unter der Annahme, dass die Unabhängigkeit für das von mir geschriebene Programm in Ordnung ist):Yiμ+αYi1σe2(Y)

p(Y|(μ,α)t)=i=21112πσe2exp(YiμαYi1)22σe2.

Nach dem Bayes'schen Theorem können wir das Produkt der obigen Dichte mit der vorherigen Dichte nehmen, und dann brauchen wir nur die Normalisierungskonstante. Meine Vermutung ist, dass dies eine Gaußsche Verteilung sein sollte, so dass wir uns am Ende um die Normalisierungskonstante kümmern können, anstatt sie explizit mit Integralen über und berechnen .μα

Dies ist der Teil, mit dem ich Probleme habe. Wie berechne ich die Multiplikation der vorherigen Dichte (die multivariat ist) und dieses Produkts der univariaten Datendichten? Der hintere Teil muss nur eine Dichte von und , aber ich kann nicht sehen, wie Sie das aus einem solchen Produkt herausholen können.μα

Alle Hinweise sind wirklich hilfreich, auch wenn Sie mich nur in die richtige Richtung weisen und ich dann die unordentliche Algebra ausführen muss (was ich bereits mehrmals versucht habe).

Als Ausgangspunkt ist hier die Form des Zählers aus Bayes 'Regel:

1(2πσe2)52πσp2exp[12σe2i=211(YiμαYi1)2μ22σp2α22σp2].

Das Problem ist, wie man sieht, dass dies auf eine Gaußsche Dichte von .(μ,α)t

Hinzugefügt

Letztendlich läuft dies auf das folgende allgemeine Problem hinaus. Wenn Sie einen quadratischen Ausdruck wie wie bringen Sie das in eine quadratische Form für eine 2x2-Matrix ? Es ist in einfachen Fällen einfach genug, aber welchen Prozess verwenden Sie, um die mittleren Schätzungen und ?

Aμ2+Bμα+Cα2+Jμ+Kα+L
(μμ^,αα^)Q(μμ^,αα^)tQμ^α^

Beachten Sie, dass ich die einfache Möglichkeit ausprobiert habe, die Matrixformel zu erweitern und dann zu versuchen, die Koeffizienten wie oben darzustellen. In meinem Fall ist das Problem, dass die Konstante Null ist und ich dann drei Gleichungen in zwei Unbekannten erhalte, so dass es unterbestimmt ist, nur Koeffizienten abzugleichen (selbst wenn ich eine symmetrische quadratische Formmatrix annehme).L

ely
quelle
Meine Antwort auf [diese Frage] ( stats.stackexchange.com/questions/22852/… ) kann hilfreich sein. Beachten Sie, dass Sie für Ihre erste Beobachtung ein Prior benötigen - die Iterationen hören dort auf.
Wahrscheinlichkeitslogik
Ich verstehe nicht, warum ich es in diesem Fall brauche. Ich soll die Zeitintervalle so behandeln, als wären sie bedingt unabhängig von der Beobachtung. Beachten Sie, dass das Produkt der Fugendichte gerade . Ich glaube nicht, dass ich hier eine sequentiell aktualisierte Formel bekommen soll, nur eine einzige Formel für das hintere . i=2..11p((μ,α)t|Y)
ely
Das "multivariate" im vorherigen steht nicht im Widerspruch zu dem "univariaten" in den Datendichten, da es sich um Dichten in den . p(α,μ)yi
Xi'an,

Antworten:

7

In meiner Antwort auf die vorherige Antwort ging es darum, zu untersuchen, wie ich die Parameter integriert habe - denn Sie werden hier genau die gleichen Integrale verwenden. Bei Ihrer Frage wird davon ausgegangen, dass die Varianzparameter bekannt sind, sodass es sich um Konstanten handelt. Sie müssen nur die Abhängigkeit von vom Zähler betrachten. Um dies zu sehen, beachten Sie, dass wir schreiben können:α,μ

= 1

p(μ,α|Y)=p(μ,α)p(Y|μ,α)p(μ,α)p(Y|μ,α)dμdα
=1(2πσe2)52πσp2exp[12σe2i=211(YiμαYi1)2μ22σp2α22σp2]1(2πσe2)52πσp2exp[12σe2i=211(YiμαYi1)2μ22σp2α22σp2]dμdα

Beachten Sie, wie wir den ersten Faktor aus dem Doppelintegral auf dem Nenner, und es bricht mit dem Zähler ab. Wir können auch die Summe der Quadrateexp [ -1 ziehen1(2πσe2)52πσp2und es wird auch abgebrochen. Das Integral, mit dem wir übrig bleiben, ist jetzt (nach Erweiterung des quadratischen Terms):exp[12σe2i=211Yi2]

=exp[10μ2+α2i=110Yi22μi=211Yi2αi=211YiYi1+2μαi=110Yi2σe2μ22σp2α22σp2]exp[10μ2+α2i=110Yi22μi=211Yi2αi=211YiYi1+2μαi=110Yi2σe2μ22σp2α22σp2]dμdα

Now we can use a general result from the normal pdf.

exp(az2+bzc)dz=πaexp(b24ac)
This follows from completing the square on az2+bz and noting that c does not depend on z. Note that the inner integral over μ is of this form with a=102σe2+12σp2 and b=i=211Yiαi=110Yiσe2 and c=α2i=110Yi22αi=211YiYi12σe2+α22σp2. After doing this integral, you will find that the remaining integral over α is also of this form, so you can use this formula again, with a different a,b,c. Then you should be able to write your posterior in the form 12π|V|12exp[12(μμ^,αα^)V1(μμ^,αα^)T] where V is a 2×2 matrix

Let me know if you need more clues.

update

(note: correct formula, should be 10μ2 instead of μ2)

if we look at the quadratic form you've written in the update, we notice there is 5 coefficients (L is irrelevant for posterior as we can always add any constant which will cancel in the denominator). We also have 5 unknowns μ^,α^,Q11,Q12=Q21,Q22. Hence this is a "well posed" problem so long as the equations are linearly independent. If we expand the quadratic (μμ^,αα^)Q(μμ^,αα^)t we get:

Q11(μμ^)2+Q22(αα^)2+2Q12(μμ^)(αα^)
=Q11μ2+2Q21μα+Q22α2(2Q11μ^+2Q12α^)μ(2Q22α^+2Q12μ^)α+
+Q11μ^2+Q22α^2+2Q12μ^α^

Comparing second order coefficient we get A=Q11,B=2Q12,C=Q22 which tells us what the (inverse) covariance matrix looks like. Also we have two slightly more complicated equations for α^,μ^ after substituting for Q. These can be written in matrix form as:

(2ABB2C)(μ^α^)=(JK)

Thus the estimates are given by:

(μ^α^)=(2ABB2C)1(JK)=14ACB2(BK2JCBJ2KA)

Showing that we do not have unique estimates unless 4ACB2. Now we have:

A=102σe2+12σp2B=i=110Yiσe2C=i=110Yi22σe2+12σp2J=i=211Yiσe2K=i=211YiYi1σe2

Note that if we define Xi=Yi1 for i=2,,11 and take the limit σp2 then the estimates for μ,α are given by the usual least squares estimate α^=i=211(YiY¯)(XiX¯)i=211(XiX¯)2 and μ^=Y¯α^X¯ where Y¯=110i=211Yi and X¯=110i=211Xi=110i=110Yi. So the posterior estimates are a weighted average between the OLS estimates and the prior estimate (0,0).

probabilityislogic
quelle
This isn't particularly helpful because I mentioned specifically that it's not the denominator that matters here. The denominator is just a normalizing constant, which will be obvious once you reduce the numerator to a Gaussian form. So tricks for evaluating the integrals in the denominator are mathematically really cool, but just not needed for my application. The only issue I need resolution with is manipulating the numerator.
ely
This answer gives you both numerator and denominator. The numerator exhibits the proper second degree polynomial in (α,μ) that leads to the normal quadratic form, as stressed by probabilityislogic.
Xi'an
@ems - by calculating the normalising constant you will construct the quadratic form required. it will contain the terms needed to compllete the square
probabilityislogic
I don't understand how this gives you the quadratic form. I've worked out the two integrals in the denominator using the Gaussian integral identity that you posted. In the end, I just get a huge, messy constant. There doesn't seem to be any clear way to take that constant and turn it into something times a determinant to the 1/2 power, etc. Not to mention I don't see how any of this explains how to calculate the new 'mean vector' (μ^,α^)t.. This is what I was asking for help for in the original question.
ely
Thanks tremendously for the detailed addition. I was making some silly errors when trying to do the algebra to figure out the quadratic form. Your comments about the relation to the OLS estimator are highly interesting and appreciated as well. I think this will speed up my code because I'll be able to draw samples from an analytic form that has built-in, optimized methods. My original plan was to use Metropolis-Hastings to sample from this, but it was very slow. Thanks!
ely