Wie kann man schnell X abtasten, wenn exp (X) ~ Gamma?

12

Ich habe ein einfaches Stichprobenproblem, bei dem meine innere Schleife wie folgt aussieht:

v = sample_gamma(k, a)

wobei sample_gammaProben aus der Gamma-Verteilung eine Dirichlet-Probe bilden.

Es funktioniert gut, aber für einige Werte von k / a läuft ein Teil der nachgeschalteten Berechnung unter.

Ich habe es angepasst, um Log Space-Variablen zu verwenden:

v = log(sample_gamma(k, a))

Nachdem Sie den Rest des Programms angepasst haben, funktioniert es ordnungsgemäß (zumindest in Testfällen erhalte ich die gleichen genauen Ergebnisse). Es ist jedoch langsamer als zuvor.

Gibt es eine Möglichkeit, direkt abzutasten , ohne langsame Funktionen wie log ( ) zu verwenden ? Ich habe versucht, dafür zu googeln, aber ich weiß nicht einmal, ob diese Distribution einen gemeinsamen Namen hat (log-gamma?).X,exp(X)Gammalog()

Luispedro
quelle
Sie müssen lediglich jede Gammavariable durch ihre Summe dividieren. Wie kommt es dann zum Unterlauf? Und wie löst die Verwendung des Logarithmus dieses Problem (Sie können die Summe nicht berechnen, ohne sie erneut zu potenzieren)?
whuber
@whuber Im Protokollbereich berechnen Sie die Summe und subtrahieren sie dann von jedem Element. So wird der erste Unterlaufpunkt vermieden. Wenn diese Dirichlets als Mischungsbestandteile dienen und wieder mit kleinen Zahlen multipliziert werden, ist eine gewisse Weiterverarbeitung erforderlich.
Luispedro
Das Hinzufügen der Protokolle ist mathematisch falsch: Es entspricht dem Multiplizieren der Gammas und nicht dem Hinzufügen. Ja, Sie erhalten möglicherweise Arbeitsergebnisse, aber sie werden definitiv keine Dirichlet-Verteilung haben! Wie genau ist der ursprüngliche Unterlauf beschaffen und welche Berechnungen führen Sie aus, wenn er auftritt? Mit welchen Werten arbeiten Sie aktuell?
whuber
@whuber Ich hätte in meiner Beschreibung vielleicht etwas zu viel vereinfacht. Ich tue für alle i {t = gamma (a, b); Summe + = t; d [i] = log (t)}; logsum = log (summe); forall i {d [i] - = logsum; }. Bisher ist dies unterlaufen, wenn a sehr klein war.
Luispedro
Verstanden: für nahe 0 wirst du in Schwierigkeiten sein, egal was passiert. Interessantes Problem! α
whuber

Antworten:

9

Betrachten wir einen kleinen Formparameter in der Nähe von 0 ist , wie beispielsweise α = 1 / 100 . Im Bereich zwischen 0 und α ist e - α ungefähr 1 , so dass das Gamma pdf ungefähr x α - 1 d x / Γ ( α ) ist . Dies kann zu einer annähernden CDF integriert werden, F α ( x ) = x ααα=1/100αeα1xα1dx/Γ(α) . Wenn wir es invertieren, sehen wir eine1/α-Potenz: einen riesigen Exponenten. FürαFα(x)=xααΓ(α)1/α Das bewirktgewisse Wahrscheinlichkeit des Unterschreitung (einen Wertdoppelter Genauigkeitweniger als 10 - 300 , mehr oder weniger). Hier ist eine Darstellung der Wahrscheinlichkeit eines Unterlaufens als Funktion des Zehn-Basis-Logarithmus von α :α=1/10010300α

Bildbeschreibung hier eingeben

Eine Lösung besteht darin, diese Näherung für die Erzeugung von Gamma-Variablen zu nutzen: Versuchen Sie im Endeffekt, eine Gamma-Variable zu erzeugen, und generieren Sie, wenn sie zu klein ist, ihren Logarithmus anhand dieser ungefähren Leistungsverteilung (siehe unten). (Wiederholen Sie diesen Vorgang, bis sich das Protokoll innerhalb des Unterlaufbereichs befindet, sodass es einen gültigen Ersatz für die ursprüngliche Unterlaufvariable darstellt.) Subtrahieren Sie für die Dirichlet-Berechnung das Maximum aller Logarithmen von jedem der Protokollwerte: Dadurch wird implizit alle neu skaliert Das Gamma variiert, sodass es die Dirichlet-Werte nicht beeinflusst. Behandeln Sie jedes resultierende Protokoll, das zu klein ist (z. B. weniger als -100), als Protokoll einer echten Null. Potenzieren Sie die anderen Protokolle. Jetzt können Sie ohne Unterlauf fortfahren.

Dies wird noch länger dauern als zuvor, aber zumindest wird es funktionieren!

Berechnen Sie C = log ( Γ ( α ) ) +, um eine ungefähre log Gamma-Variation mit dem Formparameter erstellenα. Dies ist einfach, da es Algorithmen gibt, mit denen dieWerte von log Gamma direkt berechnet werden können. Erzeugen Sie einen gleichmäßigen Zufalls-Float zwischen 0 und 1, nehmen Sie seinen Logarithmus, dividieren Sie durch α und addieren Sie C dazu.C=log(Γ(α))+log(α)αC

Da der Parameter scale die Variable lediglich neu skaliert, ist es kein Problem, sie in diesen Prozeduren unterzubringen. Sie brauchen es nicht einmal, wenn alle Skalenparameter gleich sind.

Bearbeiten

In einer weiteren Antwort beschreibt das OP ein Verfahren, bei dem die Potenz einer gleichförmigen Variation (a B ( α ) -Variate) mit einer Γ ( α + 1 ) -Variate multipliziert wird . Dies funktioniert, weil das PDF der gemeinsamen Verteilung dieser beiden Variablen gleich ( α x α - 1 ) ist ( y α e - y d y / Γ ( α + 1 ) ) . Um das pdf von z = x y zu finden1/αB(α)Γ(α+1)(αxα1)(yαeydy/Γ(α+1))z=xywir setzen , dividieren durch das jakobinische x bis , da 0 y 1 istyz/xx, und integriere . Das Integral muss im Bereich von zxz0y1

pdf(z)=αΓ(α+1)z(xα/x)ex(z/x)α1dxdz=1Γ(α)zα1ezdz,

Das ist das PDF einer -Distribution.Γ(α)

Der springende Punkt ist, dass, wenn , es unwahrscheinlich ist, dass ein aus Γ ( α + 1 ) gezogener Wert unterschritten wird, und dass sein log und 1 /0<α<1Γ(α+1) mal das Protokolls eines unabhängigen einheitliches variate werden wir das Protokoll eines haben Γ ( α ) variieren. Das Protokoll ist wahrscheinlich sehr negativ, aber wir haben die Erstellung des Antilog umgangen, der in einer Gleitkommadarstellung unterlaufen wird.1/αΓ(α)

whuber
quelle
1
Nur ein Argument, um Ihre Bearbeitung ein bisschen eleganter zu gestalten. Sie müssen hier nicht wirklich an die Integration appellieren. Verwenden Sie einfach die Tatsache, dass , und daßΓ(α)+Γ(1)~Γ(α+1). Dies sind sowohl Standardeigenschaften der Beta-Verteilung als auch der Gamma-Verteilung. Auch wenn& agr;ap0haben wir etway~expo(1)Γ(α)Γ(α)+Γ(1)Beta(α,1)Γ(α)+Γ(1)Γ(α+1)α0yexpo(1), Die schneller sein kann , simulieren ( ) als eine allgemeine Γ ( α + 1 ) Zufallsvariable. log(u)Γ(α+1)
Wahrscheinlichkeitsrechnung
7

Ich beantworte meine eigene Frage, habe aber eine ziemlich gute Lösung gefunden, auch wenn ich sie nicht vollständig verstehe. Wenn Sie sich den Code aus der GNU Scientific Library ansehen, sehen Sie hier, wie er Gammavariablen abtastet ( rist der Zufallszahlengenerator, aist undαb ist ):β

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
   }

gsl_ran_gammaist die Funktion, die eine Gamma-Zufallsstichprobe zurückgibt (obiger Aufruf ist also rekursiv), während gsl_rng_uniform_poseine gleichmäßig verteilte Zahl in (diesist streng positiv, da garantiert wird, dass sie nicht 0.0 zurückgibt).(0,1)_pos

Daher kann ich das Protokoll des letzten Ausdrucks nehmen und verwenden

return log(gsl_ran_gamma(r, 1.0 + a, b)) + log(u)/a;

Um zu bekommen, was ich wollte. Ich habe jetzt zwei log()Anrufe (aber einen weniger pow()), aber das Ergebnis ist wahrscheinlich besser. Zuvor hatte ich, wie Whuber betonte, etwas mit der Potenz von erhöht , möglicherweise eine riesige Zahl. Jetzt multipliziere ich im Logspace mit 1 / a . Es ist also weniger wahrscheinlich, dass es unterläuft.1/a1/a

Luispedro
quelle
α
Ich habe meine Antwort so bearbeitet, dass sie jetzt detaillierter ist.
Luispedro
Danke: aber was ist "r"? (Beachten Sie, dass die Rekursion begrenzt ist: Es wird höchstens ein rekursiver Aufruf ausgeführt, da a> 0 1,0 + a> 1 impliziert.)
whuber
r ist der Zufallszahlengenerator (von dem Sie die Zufallszahlen erhalten).
Luispedro
Ah, das ist klug: das Produkt von a Γ(α+1) und eine unabhängige B(α,1) variate entpuppt sich als a Γ(α) variate. I edited my reply so it points to your solution and explains why it works.
whuber