Comment échantillonner rapidement X si exp (X) ~ Gamma?

J'ai un problème d'échantillonnage simple, où ma boucle intérieure ressemble à:

v = sample_gamma(k, a)

où des sample_gammaéchantillons de la distribution Gamma pour former un échantillon de Dirichlet.

Cela fonctionne bien, mais pour certaines valeurs de k / a, certains des calculs en aval sont sous-jacents.

Je l'ai adapté pour utiliser des variables d'espace de log:

v = log(sample_gamma(k, a))

Après avoir adapté tout le reste du programme, il fonctionne correctement (au moins, il me donne les mêmes résultats exacts sur les cas de test). Cependant, il est plus lent qu'auparavant.

Existe-t-il un moyen d'échantillonner directement sans utiliser de fonctions lentes comme ? J'ai essayé de googler pour cela, mais je ne sais même pas si cette distribution a un nom commun (log-gamma?). $X, \exp(X) \sim \text{Gamma}$ $\log()$

sampling gamma-distribution luispedro
la source

Il vous suffit de diviser chaque variable gamma par leur somme. Comment se produit alors le sous-dépassement? Et comment le fait de prendre le logarithme résout-il ce problème (vous ne pouvez pas calculer la somme sans exposer à nouveau de toute façon)?

whuber

@whuber Dans l'espace journal, vous calculez la somme, puis la soustrayez de chaque élément. Ainsi, cela évite le premier point de sous-dépassement. Il y a un peu de traitement supplémentaire lorsque ces dirichlets servent de composants de mélange et sont à nouveau multipliés par de petits nombres.

luispedro

L'ajout des journaux est mathématiquement incorrect: cela correspond à multiplier les gammas plutôt qu'à les ajouter. Oui, vous pourriez obtenir des résultats de travail, mais ils n'auront certainement pas de distribution Dirichlet! Encore une fois, quelle est exactement la nature du débordement d'origine et quels calculs faites-vous quand cela se produit? Quelles sont les valeurs réelles avec lesquelles vous travaillez?

whuber

@whuber j'aurais peut-être simplifié un peu trop dans ma description. Je fais pour tout i {t = gamma (a, b); somme + = t; d [i] = log (t)}; logsum = log (somme); forall i {d [i] - = logsum; }. Auparavant, cela débordait si un était très petit.

luispedro

J'ai compris: pour près de 0, vous allez avoir des ennuis quoi qu'il arrive. Problème intéressant!

α

$\alpha$

whuber

Réponses:

Considérons un petit paramètre de forme près de 0, tel que . Dans la plage comprise entre 0 et , est d'environ , donc le pdf Gamma est d'environ . Cela peut être intégré à un CDF approximatif, . En l'inversant, nous voyons une puissance : un énorme exposant. Pour cela entraîne un risque de sous-dépassement (une valeur de double précision inférieure à , plus ou moins). Voici un graphique de la chance d'obtenir un dépassement de capacité en fonction du logarithme en base dix de $\alpha$ $\alpha = 1/100$ $\alpha$ $e^{-\alpha}$ $1$ $x^{\alpha-1}dx / \Gamma(\alpha)$ $F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$ $1/\alpha$ $\alpha = 1/100$ $10^{-300}$ $\alpha$ :

entrez la description de l'image ici

Une solution consiste à exploiter cette approximation pour générer des variables log (Gamma): en effet, essayez de générer une variable Gamma et si elle est trop petite, générez son logarithme à partir de cette distribution de puissance approximative (comme illustré ci-dessous). (Effectuez cette opération à plusieurs reprises jusqu'à ce que le journal se trouve dans la plage de sous-dépassement, afin qu'il soit un substitut valide pour la variable de sous-flux d'origine.) Pour le calcul de Dirichlet, soustrayez le maximum de tous les logarithmes de chacune des valeurs de journal: cela redimensionne implicitement tous le gamma varie donc il n'affectera pas les valeurs de Dirichlet. Traitez tout journal résultant qui est trop petit (par exemple, moins de -100) comme étant le journal d'un vrai zéro. Exponentiate les autres journaux. Vous pouvez maintenant continuer sans débordement.

Cela va prendre encore plus de temps qu'auparavant, mais au moins ça marchera!

Pour générer une variable Gamma logarithmique approximative avec le paramètre de forme , précalculez . C'est facile, car il existe des algorithmes pour calculer directement les valeurs du journal Gamma . Générer un flotteur aléatoire uniforme entre 0 et 1, prendre son logarithme, diviser par , et ajouter à elle. $\alpha$ $C = \log(\Gamma(\alpha)) + \log(\alpha)$ $\alpha$ $C$

Étant donné que le paramètre d'échelle ne fait que redimensionner la variable, il n'y a aucun problème à l'adapter dans ces procédures. Vous n'en avez même pas besoin si tous les paramètres d'échelle sont identiques.

Éditer

Dans une autre réponse, l'OP décrit une méthode dans laquelle la puissance d'une variable uniforme (une variable ) est multipliée par une variable . Cela fonctionne car le pdf de la distribution conjointe de ces deux variables est égal à . Pour trouver le pdf de nous substituons , divisons par le jacobéen et intégrons . L'intégrale doit aller de à car , d'où $1/\alpha$ $B(\alpha)$ $\Gamma(\alpha+1)$ $\left(\alpha x^{\alpha-1}\right) \left(y^{\alpha}e^{-y}dy/\Gamma(\alpha+1)\right)$ $z = xy$ $y \to z/x$ $x$ $x$ $z$ $\infty$ $0 \le y \le 1$

pdf (z) = \frac{α}{Γ (α + 1)} \int_{z}^{\infty} (x^{α} / x) e^{- x} (z / x)^{α - 1} d x d z = \frac{1}{Γ (α)} z^{α - 1} e^{- z} d z,

$\text{pdf}(z)=\frac{\alpha}{\Gamma(\alpha+1)}\int_z^{\infty}{\left(x^\alpha / x\right) e^{-x} (z/x)^{\alpha-1} dx} dz = \frac{1}{\Gamma(\alpha)}z^{\alpha-1}e^{-z}dz,$

qui est le pdf d'une distribution . $\Gamma(\alpha)$

Le fait est que lorsque , une valeur tirée de est peu susceptible de dépasser et en sommant son log et fois le log d'une variable uniforme indépendante, nous aura le journal d'une variable . Le journal est probablement très négatif, mais nous aurons contourné la construction de son antilog, qui débordera dans une représentation en virgule flottante. $0 \lt \alpha \lt 1$ $\Gamma(\alpha+1)$ $1/\alpha$ $\Gamma(\alpha)$

whuber
la source

Juste un argument pour rendre votre montage un peu plus élégant, vous n'avez pas vraiment besoin de faire appel à l'intégration ici. Utilisez simplement le fait que , plus que . Ce sont deux propriétés standard des distributions bêta et gamma. De plus, lorsque nous avons à peu près , qui peut être plus rapide à simuler ( ) qu'une variable aléatoire générale .

\frac{Γ (α)}{Γ (α) + Γ (1)} \sim B e t a (α, 1)

$\frac{\Gamma(\alpha)}{\Gamma(\alpha)+\Gamma(1)}\sim Beta(\alpha,1)$

Γ (α) + Γ (1) \sim Γ (α + 1)

$\Gamma(\alpha)+\Gamma(1)\sim \Gamma(\alpha+1)$

α \approx 0

$\alpha\approx 0$

y \sim e x p o (1)

$y\sim expo(1)$

- \log (u)

$-\log(u)$

Γ (α + 1)

$\Gamma(\alpha+1)$

Probabilogic

Je réponds à ma propre question, mais j'ai trouvé une assez bonne solution, même si je ne la comprends pas bien. En regardant le code de la bibliothèque scientifique GNU, voici comment il échantillonne les variables gamma ( rest le générateur de nombres aléatoires, aest et est ): $\alpha$ b $\beta$

  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_gammaest la fonction qui retourne un échantillon aléatoire gamma (donc ce qui précède est un appel récursif), tandis que gsl_rng_uniform_posretourne un nombre uniformément distribué en (le est pour strictement positif car il est garanti de ne pas retourner 0,0). $(0,1)$ _pos

Par conséquent, je peux prendre le journal de la dernière expression et utiliser

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

Pour obtenir ce que je voulais. J'ai maintenant deux log()appels (mais un de moins pow()), mais le résultat est probablement meilleur. Avant, comme l'a souligné Whuber, j'avais quelque chose élevé à la puissance de , potentiellement un nombre énorme. Maintenant, dans l'espace de journal, je multiplie par . Il est donc moins susceptible de déborder. $1/a$ $1/a$

luispedro
la source

Pourriez-vous expliquer ce que font gsl_rng_uniform_pos et gsl_ran_gamma? Je suppose que le premier renvoie une valeur aléatoire uniforme entre 0 et r et le second est lié à une valeur Gamma (1 + a, b) - c'est peut-être un Gamma incomplet? Dans l'ensemble, cela semble très proche de l'approximation que j'ai suggérée (sauf qu'en l'examinant, il est évident que j'ai oublié de spécifier la partie diviser par , ce qui est essentiel!)

α

$\alpha$

whuber

J'ai modifié ma réponse pour inclure plus de détails maintenant.

luispedro

Merci: mais qu'est-ce que "r"? (Notez que la récursivité est bornée: au plus un appel récursif sera effectué, car a> 0 implique 1.0 + a> 1.)

whuber

r est le générateur de nombres aléatoires (d'où vous obtenez les nombres aléatoires).

luispedro

Ah, c'est intelligent: le produit d'une et d'une variable indépendante se révèle être une variable . J'ai modifié ma réponse afin qu'elle pointe vers votre solution et explique pourquoi cela fonctionne.

Γ (α + 1)

$\Gamma(\alpha+1)$

B (α, 1)

$B(\alpha,1)$

Γ (α)

$\Gamma(\alpha)$

whuber