Wie wird von

Ich möchte nach einer Dichte wobei und sind streng positiv. (Motivation: Dies kann für die Gibbs-Abtastung nützlich sein, wenn der Formparameter einer Gammadichte eine einheitliche Priorität hat.)

Kann jemand von dieser Dichte leicht probieren? Vielleicht ist es Standard und nur etwas, von dem ich nichts weiß?

Ich kann mir einen dummen Ablehnungs-Abtastalgorithmus vorstellen, der mehr oder weniger funktioniert (finde den Modus von , probiere von Uniform in einer großen Kiste und ablehnen, wenn ), aber (i) es ist überhaupt nicht effizient und (ii) wird zu groß für einen Computer, um auch nur mäßig leicht zu handhaben großes und . (Beachten Sie, dass der Modus für großes und d ungefähr bei a = c d liegt .)$a^*$$f$$(a,u)$$[0,10a^*]\times [0,f(a^*)]$$u>f(a)$$f(a^*)$$c$$d$$c$$d$$a=cd$

Vielen Dank im Voraus für jede Hilfe!

Die Rückweisungsabtastung funktioniert außergewöhnlich gut, wenn und für c d ≥ exp ( 2 ) angemessen ist .$c d \ge \exp(5)$$c d \ge \exp(2)$

Um die Mathematik ein wenig zu vereinfachen, lassen Sie , schreiben Sie x = a und beachten Sie, dass$k = c d$$x = a$

für . Einstellen x = u 3 / 2 gibt$x \ge 1$$x = u^{3/2}$

für . Wenn k ≥ exp ( 5 ) ist , ist diese Verteilung extrem nahe an Normal (und kommt näher, wenn k größer wird). Insbesondere können Sie$u \ge 1$$k \ge \exp(5)$$k$

1. Ermitteln Sie den Modus von numerisch (z. B. mit Newton-Raphson).$f(u)$

2. Erweitern Sie in Bezug auf seinen Modus auf eine zweite Ordnung.$\log{f(u)}$

Dies ergibt die Parameter einer eng angenäherten Normalverteilung. Mit hoher Genauigkeit dominiert diese angenäherte Normale Ausnahme der extremen Schwänze. (Wenn k < exp ( 5 ) ist , müssen Sie möglicherweise das normale PDF-Dokument ein wenig vergrößern, um die Dominanz sicherzustellen.)$f(u)$$k \lt \exp(5)$

Nachdem Sie diese Vorarbeit für einen bestimmten Wert von und eine Konstante M > 1 (wie unten beschrieben) geschätzt haben, müssen Sie eine Zufallsvariable erhalten:$k$$M \gt 1$

1. Zeichnen Sie einen Wert aus der dominierenden Normalverteilung g ( u ) .$u$$g(u)$

2. Wenn oder wenn eine neue gleichförmige Variable X f ( u ) / ( M g ( u ) ) überschreitet , kehre zu Schritt 1 zurück.$u \lt 1$$X$$f(u)/(M g(u))$

3. Set .$x = u^{3/2}$

Die erwartete Anzahl von Bewertungen von aufgrund der Diskrepanzen zwischen g und f ist nur geringfügig größer als 1. (Einige zusätzliche Bewertungen werden aufgrund von Zurückweisungen von Variablen kleiner als 1 auftreten , aber selbst wenn k so niedrig wie 2 ist, ist die Häufigkeit von solchen Vorkommen ist klein.)$f$$g$$f$$1$$k$$2$

Dieses Diagramm zeigt die Logarithmen von g und f als Funktion von u für . Da die Grafiken so nahe beieinander liegen, müssen wir ihr Verhältnis überprüfen, um zu sehen, was los ist:$k=\exp(5)$

Dies zeigt das logarithmische Verhältnis ; Der Faktor M = exp ( 0,004 ) wurde einbezogen, um sicherzustellen, dass der Logarithmus im gesamten Hauptteil der Verteilung positiv ist. das heißt, es wird sichergestellt, dass M g ( u ) ≥ f ( u ) ist, außer möglicherweise in Bereichen mit vernachlässigbarer Wahrscheinlichkeit. Indem Sie M ausreichend groß machen, können Sie sicherstellen, dass M ⋅ $\log(\exp(0.004)g(u)/f(u))$$M = \exp(0.004)$$Mg(u) \ge f(u)$$M$$M \cdot g$dominiert in allen außer den extremsten Schwänzen (die ohnehin praktisch keine Chance haben, in einer Simulation ausgewählt zu werden). Je größer M ist, desto häufiger kommt es jedoch zu Ausschuss. Wenn k groß wird, kann M sehr nahe an 1 gewählt werden , was praktisch keine Nachteile mit sich bringt.$f$$M$$k$$M$$1$

Ein ähnlicher Ansatz funktioniert sogar für , aber es können ziemlich große Werte von M erforderlich sein, wenn exp ( 2 ) < k < exp ( 5 ) ist , weil f ( u ) merklich asymmetrisch ist. Zum Beispiel müssen wir mit k = exp ( 2 ) M = 1 setzen , um ein einigermaßen genaues g zu erhalten :$k \gt \exp(2)$$M$$\exp(2) \lt k \lt \exp(5)$$f(u)$$k = \exp(2)$$g$$M=1$

Die obere rote Kurve ist der Graph von während die untere blaue Kurve der Graph von log ( f ( u ) ) ist . Die Zurückweisungsabtastung von f relativ zu exp ( 1 ) g führt dazu, dass ungefähr 2/3 aller Versuchsziehungen zurückgewiesen werden, was den Aufwand verdreifacht: immer noch nicht schlecht. Die rechte tail ( u > 10 oder x > 10 3 / 2 ~ $\log(\exp(1)g(u))$$\log(f(u))$$f$$\exp(1)g$$u \gt 10$$x \gt 10^{3/2} \sim 30$) In dem Abstoßungs Probenahme (weil unterrepräsentiert seinen nicht länger vorherrscht f dort), aber das tail umfasst weniger als exp ( - 20 ) ~ 10 - 9 der Gesamtwahrscheinlichkeit.$\exp(1)g$$f$$\exp(-20) \sim 10^{-9}$

Zusammenfassend lässt sich sagen, dass Sie nach einem ersten Versuch, den Modus zu berechnen und den quadratischen Term der Potenzreihe von um den Modus herum zu bewerten - ein Versuch, der höchstens einige zehn Funktionsbewertungen erfordert -, die Ablehnungsabtastung bei verwenden können erwartete Kosten zwischen 1 und 3 (oder so) Bewertungen pro Variation. Der Kostenmultiplikator fällt schnell auf 1 ab, wenn k = c d über 5 hinaus ansteigt.$f(u)$$k = c d$

Auch wenn nur ein Draw von benötigt wird, ist diese Methode sinnvoll. Es kommt zum Tragen, wenn für den gleichen Wert von k viele unabhängige Ziehungen erforderlich sind , denn dann wird der Aufwand für die anfänglichen Berechnungen über viele Ziehungen amortisiert.$f$$k$

Nachtrag

@Cardinal hat vernünftigerweise darum gebeten, einen Teil der Hand-Waving-Analyse in der Vergangenheit zu unterstützen. Insbesondere warum soll die Transformation macht die Verteilung etwa normal?$x = u^{3/2}$

In Anbetracht der Theorie der Box-Cox-Transformationen ist es normal, eine Potenztransformation der Form (für eine Konstante α , die sich hoffentlich nicht zu stark von der Einheit unterscheidet) anzustreben, die eine Verteilung "normaler" macht. Denken Sie daran, dass alle Normalverteilungen einfach charakterisiert werden: Die Logarithmen ihrer pdfs sind rein quadratisch, mit einem linearen Term von Null und keinen Termen höherer Ordnung. Daher können wir jedes PDF mit einer Normalverteilung vergleichen, indem wir seinen Logarithmus als Potenzreihe um seinen (höchsten) Peak erweitern. Wir suchen einen Wert von α , der (mindestens) den dritten Wert ergibt$x = u^\alpha$$\alpha$$\alpha$Macht schwindet, zumindest annähernd: Das ist das Höchste, was wir zu Recht hoffen können, dass ein einziger freier Koeffizient erreicht wird. Oft funktioniert das gut.

Aber wie bekommt man diese bestimmte Distribution in den Griff? Nach der Leistungsumwandlung ist das PDF

Nimm seinen Logarithmus und verwende Stirlings asymptotische Expansion von :$\log(\Gamma)$

(für kleine Werte von , die nicht konstant sind). Dies funktioniert, vorausgesetzt α ist positiv, was wir als der Fall annehmen werden (ansonsten können wir den Rest der Erweiterung nicht vernachlässigen).$c$$\alpha$

Berechnen ihre dritte Ableitung (die, wenn sie durch unterteilt , Wird der Koeffizient der dritten Potenz sein , u in der Potenzreihe) und nutzt die Tatsache aus, dass an der Spitze, die erste Ableitung gleich Null sein muss. Dies vereinfacht die dritte Ableitung erheblich und gibt (ungefähr, weil wir die Ableitung von c ignorieren )$3!$$u$$c$

Wenn nicht zu klein ist, ist u in der Tat am Gipfel groß. Da α positiv ist, ist der dominante Term in diesem Ausdruck die 2 α- Potenz, die wir auf Null setzen können, indem wir ihren Koeffizienten verschwinden lassen:$k$$u$$\alpha$$2\alpha$

Deshalb funktioniert so gut: Mit dieser Wahl wird der Koeffizient des kubischen Begriffs um die Spitze verhält sich wie u - 3 , das in der Nähe ist exp ( - 2 k ) . Sobald k ungefähr 10 überschreitet, können Sie es praktisch vergessen, und es ist sogar für k bis zu 2 einigermaßen klein . Die höheren Potenzen spielen ab dem vierten eine immer geringere Rolle, wenn k groß wird, weil ihre Koeffizienten zunehmen auch proportional kleiner. Im Übrigen gelten die gleichen Berechnungen (basierend auf der zweiten Ableitung von l o g ( $\alpha = 3/2$$u^{-3}$$\exp(-2 k)$$k$$k$$k$ an seiner Spitze) zeigen, dass die Standardabweichung dieser Normalen Näherung etwas kleiner als$log(f(u))$, wobei der Fehler proportional zuexp(-k/2) ist.$\frac{2}{3}\exp(k/6)$$\exp(-k/2)$

Ich mag die Antwort von @ whuber sehr; Es ist wahrscheinlich sehr effizient und hat eine schöne Analyse. Es erfordert jedoch einige tiefe Einsichten in Bezug auf diese bestimmte Verteilung. In Situationen, in denen Sie diese Einsicht nicht haben (also für verschiedene Distributionen), gefällt mir auch der folgende Ansatz, der für alle Distributionen funktioniert, bei denen das PDF doppelt differenzierbar ist und diese zweite Ableitung endlich viele Wurzeln hat. Es erfordert eine Menge Arbeit, um es einzurichten, aber danach haben Sie eine Engine, die für die meisten Distributionen funktioniert, die Sie darauf werfen können.

Grundsätzlich besteht die Idee darin, eine stückweise lineare Obergrenze für das PDF zu verwenden, die Sie anpassen, wenn Sie das Ablehnungssampling durchführen. Gleichzeitig haben Sie eine stückweise lineare Absenkunggebunden für das PDF, so dass Sie das PDF nicht zu häufig auswerten müssen. Die oberen und unteren Grenzen werden durch Akkorde und Tangenten an das PDF-Diagramm angegeben. Die anfängliche Unterteilung in Intervalle erfolgt so, dass das PDF in jedem Intervall entweder konkav oder konvex ist. Wenn Sie einen Punkt (x, y) ablehnen müssen, unterteilen Sie dieses Intervall in x. (Sie können bei x auch eine zusätzliche Unterteilung vornehmen, wenn Sie die PDF-Datei berechnen mussten, da die Untergrenze wirklich schlecht ist.) Dadurch treten die Unterteilungen besonders häufig auf, wenn die Ober- (und Untergrenze) schlecht sind, sodass Sie eine wirklich gute erhalten Annäherung Ihres PDF im Wesentlichen kostenlos. Die Details sind etwas schwierig zu verstehen, aber ich habe versucht, die meisten davon in dieser Reihe von Blog- Beiträgen zu erklären - insbesondereder letzte .

In diesen Beiträgen wird nicht erläutert, was zu tun ist, wenn die PDF-Datei weder in Domänen noch in Werten beschränkt ist. Ich würde die etwas offensichtliche Lösung empfehlen, entweder eine Transformation durchzuführen, die sie endlich macht (was schwer zu automatisieren ist), oder einen Cutoff zu verwenden. Ich würde den Cutoff in Abhängigkeit von der Gesamtzahl der Punkte wählen, die Sie erwarten, sagen wir N , und den Cutoff so wählen, dass der entfernte Teil eine Wahrscheinlichkeit von weniger als . (Dies ist einfach genug, wenn Sie ein geschlossenes Formular für die CDF haben; andernfalls kann es auch schwierig sein.)$1 / (10 N)$

Diese Methode ist in Maple als Standardmethode für benutzerdefinierte kontinuierliche Verteilungen implementiert. (Vollständige Offenlegung - Ich arbeite für Maplesoft.)

Ich habe einen Beispiellauf durchgeführt und dabei 10 ^ 4 Punkte für c = 2, d = 3 generiert, wobei [1, 100] als Anfangsbereich für die Werte angegeben wurde:

Es gab 23 Ablehnungen (in rot), 51 Punkte "auf Bewährung", die zu der Zeit zwischen der Untergrenze und dem tatsächlichen PDF lagen, und 9949 Punkte, die akzeptiert wurden, nachdem nur lineare Ungleichungen überprüft wurden. Das sind 74 Bewertungen des PDFs insgesamt oder ungefähr eine PDF-Bewertung pro 135 Punkte. Das Verhältnis sollte besser werden, wenn Sie mehr Punkte generieren, da die Approximation immer besser wird (und umgekehrt, wenn Sie nur wenige Punkte generieren, ist das Verhältnis schlechter).

Sie können dies tun, indem Sie die Inversionsmethode numerisch ausführen. Wenn Sie einheitliche (0,1) Zufallsvariablen in die inverse CDF einfügen, erhalten Sie ein Unentschieden von der Verteilung. Ich habe unten einen R-Code eingefügt, der dies bewirkt, und nach den wenigen Überprüfungen, die ich durchgeführt habe, funktioniert er gut, aber er ist ein bisschen schlampig und ich bin sicher, dass Sie ihn optimieren können.

Wenn Sie nicht mit R vertraut sind, ist lgamma () das Protokoll der Gammafunktion. integriere () berechnet ein bestimmtes 1-D Integral; uniroot () berechnet eine Wurzel einer Funktion unter Verwendung einer 1-D-Halbierung.

# density. using the log-gamma gives a more numerically stable return for 
# the subsequent numerical integration (will not work without this trick)
f = function(x,c,d) exp( x*log(c) + (x-1)*log(d) - lgamma(x) )

# brute force calculation of the CDF, calculating the normalizing constant numerically
F = function(x,c,d) 
{
   g = function(x) f(x,c,d)
   return( integrate(g,1,x)$val/integrate(g,1,Inf)$val )
}

# Using bisection to find where the CDF equals p, to give the inverse CDF. This works 
# since the density given in the problem corresponds to a continuous CDF. 
F_1 = function(p,c,d) 
{
   Q = function(x) F(x,c,d)-p
   return( uniroot(Q, c(1+1e-10, 1e4))$root )
}

# plug uniform(0,1)'s into the inverse CDF. Testing for c=3, d=4. 
G = function(x) F_1(x,3,4)
z = sapply(runif(1000),G)

# simulated mean
mean(z)
[1] 13.10915

# exact mean
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*x/nc
integrate(h,1,Inf)$val
[1] 13.00002 

# simulated second moment
mean(z^2)
[1] 183.0266

# exact second moment
g = function(x) f(x,3,4)
nc = integrate(g,1,Inf)$val
h = function(x) f(x,3,4)*(x^2)/nc
integrate(h,1,Inf)$val
[1] 181.0003

# estimated density from the sample
plot(density(z))

# true density 
s = seq(1,25,length=1000)
plot(s, f(s,3,4), type="l", lwd=3)

$(1,10000)$$> 100000$$c, d$

$cd$