Comment échantillonner à partir de ?

Je veux échantillonner selon une densité où et sont strictement positifs. (Motivation: cela pourrait être utile pour l'échantillonnage de Gibbs lorsque le paramètre de forme d'une densité gamma a une priorité uniforme.)

F (une) \propto \frac{c^{une} {ré}^{une - 1}}{Γ (une)} 1_{(1, \infty)} (une)

$f(a) \propto \frac{c^a d^{a-1}}{\Gamma(a)} 1_{(1,\infty)}(a)$

c

$c$

d

$d$

Quelqu'un sait-il facilement échantillonner à partir de cette densité? Peut-être que c'est standard et juste quelque chose que je ne sais pas?

Je peux penser à un stupide algorithme d'échantillonnage de rejet qui fonctionnera plus ou moins (trouver le mode de , échantillon d'uniforme dans une grande boîte et rejeter si ), mais (i) ce n'est pas du tout efficace et (ii) sera trop gros pour qu'un ordinateur puisse le gérer facilement même modérément grand et . (Notez que le mode pour les grands et est approximativement à .) $a^*$ $f$ $(a,u)$ $[0,10a^*]\times [0,f(a^*)]$ $u>f(a)$ $f(a^*)$ $c$ $d$ $c$ $d$ $a=cd$

Merci d'avance pour votre aide!

+1 bonne question. Je ne sais pas s'il existe une approche standard.

suncoolsu

Avez-vous déjà vérifié (des idées) dans les endroits "évidents", comme, par exemple, le texte de Devroye ?

cardinal

Oui, j'ai déjà essayé un certain nombre d'idées du texte de Devroye. Le

a rendu difficile l'accès à la plupart d'entre eux, bien que ... la plupart des approches semblent nécessiter soit l'intégration (pour trouver le cdf), la décomposition en fonctions plus simples, soit la délimitation par des fonctions plus simples ... mais la

fonction rend toutes ces difficultés. Si quelqu'un a des idées sur où chercher des approches pour ces sous-problèmes - par exemple, où la fonction

apparaît-elle de manière "essentielle" comme ici (pas seulement comme constante de normalisation) dans les statistiques - cela pourrait être très utile pour moi !

Γ (a)

$\Gamma(a)$

Γ

$\Gamma$

Γ

$\Gamma$

Il y a une énorme différence entre le cas

. Devez-vous couvrir ces deux cas?

c d < 2

$c d \lt 2$

c d \geq 2

$c d \ge 2$

whuber

C'est vrai - merci. On peut supposer que

c d \geq 2

$cd\geq 2$

Réponses:

L'échantillonnage de rejet fonctionnera exceptionnellement bien lorsque et est raisonnable pour . $c d \ge \exp(5)$ $c d \ge \exp(2)$

Pour simplifier un peu les mathématiques, soit , écrivez et notez que $k = c d$ $x = a$

F (X) \propto \frac{k^{X}}{Γ (X)} ré X

$f(x) \propto \frac{k^x}{\Gamma(x)} dx$

pour . Réglage donne $x \ge 1$ $x = u^{3/2}$

F (u) \propto \frac{k^{u^{3 / 2}}}{Γ (u^{3 / 2})} u^{1 / 2} ré u

$f(u) \propto \frac{k^{u^{3/2}}}{\Gamma(u^{3/2})} u^{1/2} du$

pour . Lorsque , cette distribution est extrêmement proche de la normale (et se rapproche à mesure que grandit). Plus précisément, vous pouvez $u \ge 1$ $k \ge \exp(5)$ $k$

Trouvez le mode de numériquement (en utilisant, par exemple, Newton-Raphson). $f(u)$
Développez au deuxième ordre sur son mode. $\log{f(u)}$

Cela donne les paramètres d'une distribution normale très approximative. Avec une grande précision, cette normale approximative domine sauf dans les queues extrêmes. (Lorsque , vous devrez peut-être augmenter légèrement le pdf normal pour assurer la domination.) $f(u)$ $k \lt \exp(5)$

Ayant effectué ce travail préliminaire pour une valeur donnée de et ayant estimé une constante (comme décrit ci-dessous), l'obtention d'une variable aléatoire est une question de: $k$ $M \gt 1$

Tirez une valeur de la distribution normale dominante . $u$ $g(u)$
Si ou si une nouvelle variable uniforme dépasse , retournez à l'étape 1. $u \lt 1$ $X$ $f(u)/(M g(u))$
Set . $x = u^{3/2}$

Le nombre prévu d'évaluations de raison des écarts entre et n'est que légèrement supérieur à 1. (Certaines évaluations supplémentaires se produiront en raison des rejets de variances inférieures à , mais même lorsque est aussi faible que la fréquence de telles les occurrences sont petites.) $f$ $g$ $f$ $1$ $k$ $2$

Tracé de f et g pour k = 5

Ce graphique montre les logarithmes de g et f en fonction de u pour . Parce que les graphiques sont si proches, nous devons inspecter leur ratio pour voir ce qui se passe: $k=\exp(5)$

tracé du rapport de log

Ceci affiche le rapport logarithmique ; le facteur été inclus pour garantir que le logarithme est positif dans toute la partie principale de la distribution; c'est-à-dire pour assurer sauf éventuellement dans des régions de probabilité négligeable. En rendant suffisamment grand, vous pouvez garantir que $\log(\exp(0.004)g(u)/f(u))$ $M = \exp(0.004)$ $Mg(u) \ge f(u)$ $M$ $M \cdot g$ domine dans toutes les queues sauf les plus extrêmes (qui n'ont pratiquement aucune chance d'être choisies dans une simulation de toute façon). Cependant, plus est grand, plus les rejets se produisent fréquemment. Lorsque devient grand, peut être choisi très près de , ce qui n'entraîne pratiquement aucune pénalité. $f$ $M$ $k$ $M$ $1$

Une approche similaire fonctionne même pour , mais des valeurs assez grandes de peuvent être nécessaires lorsque , car est sensiblement asymétrique. Par exemple, avec , pour obtenir un raisonnablement précis, nous devons définir : $k \gt \exp(2)$ $M$ $\exp(2) \lt k \lt \exp(5)$ $f(u)$ $k = \exp(2)$ $g$ $M=1$

Tracer pour k = 2

La courbe rouge supérieure est le graphique de tandis que la courbe bleue inférieure est le graphique de . L'échantillonnage de rejet de rapport à entraînera le rejet d' environ 2/3 de tous les tirages d'essai, triplant l'effort: toujours pas mal. La queue droite ( ou $\log(\exp(1)g(u))$ $\log(f(u))$ $f$ $\exp(1)g$ $u \gt 10$ $x \gt 10^{3/2} \sim 30$ ) sera sous-représentée dans l'échantillonnage de rejet (car n'y domine plus ), mais cette queue comprend moins que de la probabilité totale. $\exp(1)g$ $f$ $\exp(-20) \sim 10^{-9}$

Pour résumer, après un effort initial pour calculer le mode et évaluer le terme quadratique de la série de puissance de autour du mode - un effort qui nécessite au plus quelques dizaines d'évaluations de fonction - vous pouvez utiliser l'échantillonnage de rejet à un coût prévu compris entre 1 et 3 (ou plus) évaluations par variable. Le multiplicateur de coût tombe rapidement à 1 lorsque augmente au-delà de 5. $f(u)$ $k = c d$

Même lorsqu'un seul tirage de est nécessaire, cette méthode est raisonnable. Il prend tout son sens lorsque de nombreux tirages indépendants sont nécessaires pour la même valeur de , car le surcoût des calculs initiaux est amorti sur de nombreux tirages. $f$ $k$

Addenda

@Cardinal a demandé, tout à fait raisonnablement, le soutien d'une partie de l'analyse de la main dans ce qui précède. En particulier, pourquoi la transformation faire la distribution à peu près normale? $x = u^{3/2}$

À la lumière de la théorie des transformations de Box-Cox , il est naturel de rechercher une transformation de puissance de la forme (pour une constante , espérons-le pas trop différente de l'unité) qui rendra une distribution "plus" normale. Rappelons que toutes les distributions normales sont simplement caractérisées: les logarithmes de leurs pdfs sont purement quadratiques, avec zéro terme linéaire et aucun terme d'ordre supérieur. Par conséquent, nous pouvons prendre n'importe quel pdf et le comparer à une distribution normale en étendant son logarithme en tant que série de puissance autour de son pic (le plus élevé). Nous recherchons une valeur de qui fait (au moins) le troisième $x = u^\alpha$ $\alpha$ $\alpha$ la puissance s'évanouit, au moins approximativement: c'est le plus que l'on puisse raisonnablement espérer qu'un seul coefficient libre accomplira. Cela fonctionne souvent bien.

Mais comment maîtriser cette distribution particulière? En effectuant la transformation de puissance, son pdf est

F (u) = \frac{k^{u^{α}}}{Γ (u^{α})} u^{α - 1} .

$f(u) = \frac{k^{u^{\alpha}}}{\Gamma(u^{\alpha})} u^{\alpha-1}.$

Prenez son logarithme et utilisez l'expansion asymptotique de Stirling de : $\log(\Gamma)$

\log (f (u)) \approx \log (k) u^{α} + (α - 1) Journal (u) - α u^{α} Journal (u) + u^{α} - Journal (2 π u^{α}) / 2 + c u^{- α}

$\log(f(u)) \approx \log(k) u^\alpha + (\alpha - 1)\log(u) - \alpha u^\alpha \log(u) + u^\alpha - \log(2 \pi u^\alpha)/2 + c u^{-\alpha}$

(pour les petites valeurs de , ce qui n'est pas constant). Cela fonctionne à condition que soit positif, ce que nous supposerons être le cas (car sinon nous ne pouvons pas négliger le reste de l'expansion). $c$ $\alpha$

Calculez sa dérivée troisième (qui, lorsqu'elle est divisée par Sera le coefficient de la troisième puissance de dans la série de puissance) et exploitez le fait qu'au sommet, la dérivée première doit être nulle. Cela simplifie considérablement la troisième dérivée, donnant (approximativement, parce que nous ignorons la dérivée de ) $3!$ $u$ $c$

- \frac{1}{2} u^{- (3 + α)} α (2 α (2 α - 3) u^{2 α} + (α^{2} - 5 α + 6) u^{α} + 12 c α) .

$-\frac{1}{2} u^{-(3+\alpha)} \alpha \left(2 \alpha(2 \alpha-3) u^{2 \alpha} + (\alpha^2 - 5\alpha +6)u^\alpha + 12 c \alpha \right).$

Lorsque n'est pas trop petit, sera en effet grand au sommet. Parce que est positif, le terme dominant dans cette expression est la puissance , que nous pouvons mettre à zéro en faisant disparaître son coefficient: $k$ $u$ $\alpha$ $2\alpha$

2 α - 3 = 0.

$2 \alpha-3 = 0.$

Voilà pourquoi les fonctionne si bien: avec ce choix, le coefficient du terme cubique autour du pic se comporte comme , qui est proche de . Une fois que dépasse 10 environ, vous pouvez pratiquement l'oublier, et il est raisonnablement petit, même pour jusqu'à 2. Les puissances supérieures, à partir du quatrième, jouent de moins en moins un rôle lorsque devient grand, car leurs coefficients augmentent proportionnellement plus petit aussi. Par ailleurs, les mêmes calculs (basés sur la dérivée seconde de $\alpha = 3/2$ $u^{-3}$ $\exp(-2 k)$ $k$ $k$ $k$ à son apogée) montre que l'écart-type de cette approximation normale est légèrement inférieur à $log(f(u))$ , avec l'erreur proportionnelle à. $\frac{2}{3}\exp(k/6)$ $\exp(-k/2)$

whuber
la source

(+1) Excellente réponse. Peut-être pourriez-vous développer brièvement la motivation de votre choix de variable de transformation.

cardinal

Belle addition. Cela fait une réponse très, très complète!

cardinal

J'aime beaucoup la réponse de @ whuber; il est susceptible d'être très efficace et a une belle analyse. Mais cela nécessite une compréhension approfondie de cette distribution particulière. Pour les situations où vous n'avez pas cette idée (donc pour différentes distributions), j'apprécie également l'approche suivante qui fonctionne pour toutes les distributions où le PDF est deux fois différenciable et que la dérivée seconde a un nombre fini de racines. Cela nécessite un peu de travail à configurer, mais ensuite vous avez un moteur qui fonctionne pour la plupart des distributions que vous pouvez lui lancer.

Fondamentalement, l'idée est d'utiliser une limite supérieure linéaire par morceaux au PDF que vous adaptez lorsque vous effectuez un échantillonnage de rejet. Dans le même temps, vous avez une baisse linéaire par morceauxà destination du PDF, ce qui vous évite d'avoir à évaluer le PDF trop fréquemment. Les limites supérieures et inférieures sont données par des accords et des tangentes au graphique PDF. La division initiale en intervalles est telle qu'à chaque intervalle, le PDF est soit tout concave, soit tout convexe; chaque fois que vous devez rejeter un point (x, y), vous subdivisez cet intervalle en x. (Vous pouvez également effectuer une subdivision supplémentaire en x si vous deviez calculer le PDF car la limite inférieure est vraiment mauvaise.) Cela rend les subdivisions particulièrement fréquentes lorsque les limites supérieures (et inférieures) sont mauvaises, vous obtenez donc une très bonne approximation de votre PDF essentiellement gratuitement. Les détails sont un peu difficile à obtenir à droite, mais j'ai essayé d'expliquer la plupart d'entre eux dans cette série de blogs messages - en particulierle dernier .

Ces messages ne discutent pas de ce qu'il faut faire si le PDF n'est pas limité, ni dans le domaine ni dans les valeurs; Je recommanderais la solution quelque peu évidente soit de faire une transformation qui les rend finis (ce qui serait difficile à automatiser) soit d'utiliser une coupure. Je choisirais la coupure en fonction du nombre total de points que vous prévoyez de générer, disons N , et je choisirais la coupure de sorte que la partie retirée ait moins de probabilité. (C'est assez facile si vous avez un formulaire fermé pour le CDF; sinon cela pourrait aussi être délicat.) $1 / (10 N)$

Cette méthode est implémentée dans Maple comme méthode par défaut pour les distributions continues définies par l'utilisateur. (Divulgation complète - Je travaille pour Maplesoft.)

J'ai fait un exemple, générant 10 ^ 4 points pour c = 2, d = 3, en spécifiant [1, 100] comme plage initiale pour les valeurs:

graphique

Il y a eu 23 rejets (en rouge), 51 points "en probation" qui se situaient à l'époque entre la borne inférieure et le PDF réel, et 9949 points qui ont été acceptés après avoir vérifié uniquement les inégalités linéaires. C'est 74 évaluations du PDF au total, soit environ une évaluation PDF pour 135 points. Le rapport devrait s'améliorer à mesure que vous générez plus de points, car l'approximation s'améliore de plus en plus (et inversement, si vous ne générez que peu de points, le rapport est pire).

Erik P.
la source

Et en passant - si vous avez besoin d'évaluer le PDF uniquement très rarement parce que vous avez une bonne limite inférieure pour cela, vous pouvez vous permettre de prendre plus de temps pour cela, vous pouvez donc simplement utiliser une bibliothèque bignum (peut-être même MPFR?) Et évaluer la fonction Gamma en cela sans trop craindre de débordement.

Erik P.

(+1) C'est une belle approche. Merci de le partager.

whuber

Le problème de débordement est géré en exploitant les relations (simples) entre Gammas. L'idée est qu'après normalisation du pic autour de

, les seuls calculs qui importent sont de la forme

où

est assez proche de

1

$1$

Γ (\exp (c d)) / Γ (x)

$\Gamma(\exp(c d))/\Gamma(x)$

x

$x$

\exp (k)

$\exp(k)$

Γ

$\Gamma$

1

$1$

2

$2$

whuber

@whuber re: Gammas: Ah oui - je vois que vous l'avez également suggéré ci-dessus. Merci!

Erik P.

Vous pouvez le faire en exécutant numériquement la méthode d'inversion, qui dit que si vous branchez des variables aléatoires uniformes (0,1) dans le CDF inverse, vous obtenez un tirage de la distribution. J'ai inclus un code R ci-dessous qui fait cela, et d'après les quelques vérifications que j'ai faites, cela fonctionne bien, mais c'est un peu bâclé et je suis sûr que vous pouvez l'optimiser.

Si vous n'êtes pas familier avec R, lgamma () est le journal de la fonction gamma; integr () calcule une intégrale 1-D définie; uniroot () calcule la racine d'une fonction en utilisant la bissection 1-D.

# 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$

Macro
la source

La méthode est correcte, mais terriblement douloureuse! Combien d'évaluations de fonction pensez-vous être nécessaires pour une seule variable aléatoire? Milliers? Des dizaines de milliers?

whuber

c d

$cd$

(c d)^{x}

$(cd)^{x}$

x

$x$

f

$f$

a \log (c d) - \log (Γ (a))

$a \log(c d) - \log(\Gamma(a))$

C'est ce que je fais pour le calcul - cela n'évite toujours pas le débordement. Vous ne pouvez pas exposer un nombre supérieur à environ 500 sur un ordinateur. Cette quantité devient beaucoup plus grande que cela. Je veux dire "assez bien" en le comparant à l'échantillonnage de rejet mentionné par l'OP.

Macro

c d

$cd$