Wie wird von

19

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.)

f(ein)ceindein-1Γ(ein)1(1,)(ein)
cd

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 .)einf(ein,u)[0,10ein]×[0,f(ein)]u>f(ein)f(ein)cdcdein=cd

Vielen Dank im Voraus für jede Hilfe!

NF
quelle
+1 gute Frage. Ich bin nicht sicher, ob es einen Standardansatz gibt.
Suncoolsu
Haben Sie (nach Ideen) an den "offensichtlichen" Stellen nachgesehen , z. B. bei Devroye ?
Kardinal
Ja, ich habe bereits eine Reihe von Ideen aus Devroye's Text ausprobiert. Das hat es für mich schwierig gemacht, mit den meisten von ihnen irgendwohin zu kommen ... die meisten Ansätze scheinen entweder Integration (um die cdf zu finden), Zerlegung in einfachere Funktionen oder Begrenzung durch einfachere Funktionen zu erfordern ... aber die Γ- Funktion macht all dies schwierig. Wenn jemand Ideen hat, wo er nach Ansätzen für diese Teilprobleme suchen kann - z. B. wo sonst wird die Γ- Funktion in Statistiken auf eine "wesentliche" Weise angezeigt (nicht nur als Normalisierungskonstante) -, könnte dies für mich sehr hilfreich sein ! Γ(ein)ΓΓ
NF
Es gibt einen großen Unterschied zwischen dem Fall und c d 2 . Müssen Sie beide Fälle abdecken? cd<2cd2
Whuber
1
Das stimmt - danke. Wir können annehmen, dass . cd2
NF

Antworten:

21

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

Um die Mathematik ein wenig zu vereinfachen, lassen Sie , schreiben Sie x = a und beachten Sie, dassk=cdx=ein

f(x)kxΓ(x)dx

für . Einstellen x = u 3 / 2 gibtx1x=u3/2

f(u)ku3/2Γ(u3/2)u1/2du

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 Sieu1kexp(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.Logf(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<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:kM>1

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

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

  3. Set .x=u3/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.)fGf1k2

Plot von f und g für k = 5

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)

Diagramm des log-Verhältnisses

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 gLog(exp(0,004)G(u)/f(u))M=exp(0,004)Mg(u)f(u)MMgdominiert 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.fMkM1

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>exp(2)Mexp(2)<k<exp(5)f(u)k=exp(2)gM=1

Auftragung für k = 2

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 ~ 30log(exp(1)g(u))log(f(u))fexp(1)gu>10x>103/230) 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)gfexp(20)109

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=cd

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.fk


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=u3/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 ergibtx=uααα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

f(u)=kuαΓ(uα)uα1.

Nimm seinen Logarithmus und verwende Stirlings asymptotische Expansion von :log(Γ)

log(f(u))log(k)uα+(α1)log(u)αuαlog(u)+uαlog(2πuα)/2+cuα

(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α

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!uc

12u(3+α)α(2α(2α3)u2α+(α25α+6)uα+12cα).

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:kuα2α

2α3=0.

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 ( fα=3/2u3exp(2k)kkk an seiner Spitze) zeigen, dass die Standardabweichung dieser Normalen Näherung etwas kleiner als 2 istlog(f(u)), wobei der Fehler proportional zuexp(-k/2) ist.23exp(k/6)exp(k/2)

whuber
quelle
(+1) Gute Antwort. Vielleicht könnten Sie kurz auf die Motivation für Ihre Wahl der Transformationsvariablen eingehen.
Kardinal
Schöne Ergänzung. Dies ist eine sehr, sehr vollständige Antwort!
Kardinal
11

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/(10N)

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:

Graph

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).

Erik P.
quelle
Übrigens - wenn Sie das PDF nur sehr selten auswerten müssen, weil Sie eine gute Untergrenze dafür haben, können Sie es sich leisten, länger dafür zu brauchen, also können Sie einfach eine Bignum-Bibliothek (vielleicht sogar MPFR?) Verwenden und auswerten die Gamma-Funktion in diesem ohne zu viel Angst vor Überlauf.
Erik P.
(+1) Dies ist ein schöner Ansatz. Danke, dass du es geteilt hast.
Whuber
1Γ(exp(cd))/Γ(x)xexp(k)Γ12
whuber
@whuber re: Gammas: Ah ja - ich sehe, dass du das oben auch vorgeschlagen hast. Vielen Dank!
Erik P.
3

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)>100000c,d

cd

Makro
quelle
1
Die Methode ist richtig, aber furchtbar schmerzhaft! Wie viele Funktionsbewertungen werden Ihrer Meinung nach für eine einzelne Zufallsvariable benötigt? Tausende? Zehntausende?
whuber
cd(cd)xx
1
Eine Minute für 1.000 Variationen ist nicht sehr gut: Sie werden Stunden auf eine gute Monte-Carlo-Simulation warten. Sie können mit der Zurückweisungsabtastung vier Größenordnungen schneller arbeiten. Der Trick besteht darin, mit einer engen Annäherung von f abzulehnenfeinLog(cd)-Log(Γ(ein))
Das ist, was ich für die Berechnung mache - es vermeidet immer noch keinen Überlauf. Sie können auf einem Computer keine größere Zahl als etwa 500 potenzieren. Diese Menge wird viel größer. Ich meine "ziemlich gut", wenn ich es mit der Ablehnungsprobe vergleiche, die das OP erwähnte.
Makro
1
cd