Comment puis-je échantillonner à partir d'une distribution avec CDF incomputable?

Problème lié à la simulation semi-informatique ici.

J'ai une distribution où

P (x) =$\frac{(e^b-1) e^{b (n-x)}}{e^{b n+b}-1}$

pour certaines constantes b et n, et x est un entier tel que .$0\leq x \leq n$

Maintenant, je dois échantillonner à partir de cette distribution. Il a un CDF inversible, il est donc possible de le faire directement en théorie. Le problème est que les nombres impliqués sont GRANDS. En fait, ils sont tellement volumineux qu'ils débordent tous les deux de formats conventionnels et prennent au moins quelques minutes (à un moment donné, j'ai abandonné ...) pour calculer en utilisant des formats de précision arbitraires. Fondamentalement, le CDF inverse implique toujours un terme de , pour . Malgré cela, les numéros de sortie seront toujours compris entre , il semble donc qu'il devrait y avoir un moyen de le faire.$e^{b(n+1)}$$350 < n < 3500$$0-n$

Ce que je recherche, c'est un moyen approximatif d'échantillonnage de cette distribution qui soit calculable. Existe-t-il d'autres méthodes d'échantillonnage? Que sont-ils?

Le CDF est facilement inversible. Une formule pour l'inversion conduit à ce qui doit être l'une des solutions les plus simples et les plus rapides possibles.

Commençons par observer que la probabilité du résultat , , est proportionnelle à . Ainsi, si nous générons une valeur uniforme comprise entre et = , il suffit de trouver le plus grand pour lequel$k$$0 \le k \le n$$e^{-b k}$$q$$0$$q_{\max}=\sum_{k=0}^{n} e^{-b k}$$(1 - e^{-b(n+1)})/(1 - e^{-b})$$k$

L'algèbre simple donne la solution

Voici une Rimplémentation construite comme tous les autres générateurs de nombres aléatoires: son premier argument spécifie le nombre de valeurs iid à générer et le reste des arguments nomme les paramètres ( as et as ):$b$b$n$n.max

rgeom.truncated <- function(n=1, b, n.max) {
  a <- 1 - exp(-b)
  q.max <- (1 - exp(-b*(n.max+1))) / a
  q <- runif(n, 0, q.max)
  return(-ceiling(log(1 - q*a) / b))
}

À titre d'exemple de son utilisation, générons un million de variations aléatoires selon cette distribution:

b <- 0.001
n.max <- 3500
n.sim <- 10^6
set.seed(17)
system.time(sim <- rgeom.truncated(n.sim, b,n.max))

( seconde était nécessaire.)$0.10$

h <- hist(sim+1, probability=TRUE, breaks=50, xlab="Outcome+1")
pmf <- exp(-b * (0: n.max)); pmf <- pmf / sum(pmf)
lines(0:n.max, pmf, col="Red", lwd=2)

( été ajouté à chaque valeur afin de créer un meilleur histogramme: la procédure de a une idiosyncrasie (= bug) dans laquelle la première barre est trop haute lorsque le point d'extrémité gauche est réglé à zéro.) La courbe rouge est la distribution de référence que cette simulation tente de reproduire. Évaluons la qualité de l'ajustement avec un test du chi carré:$1$Rhist

observed <- table(sim)
expected <- n.sim * pmf
chi.square <- (observed-expected)^2 / expected
pchisq(sum(chi.square), n.max, lower.tail=FALSE)

La valeur de p est de : un bel ajustement.$0.84$

Vous avez affaire à une distribution géométrique tronquée avec . Il existe différentes manières d'aborder ce problème.$p = 1-e^{-b}$

Je conseillerais différentes options dans différentes situations; certaines options impliqueraient de simuler à partir d'une géométrie et de se régénérer lorsqu'elle est en dehors de la plage, en prenant la partie entière d'une exponentielle tronquée appropriée ( comme ici ), ou en utilisant l'une quelconque de plusieurs techniques rapides adaptées à des distributions discrètes sur une plage finie. Étant donné que est grand, la prise de parole d'une exponentielle tronquée est susceptible d'être relativement rapide, mais son choix dépend également de .$n$$b$

Voici une question connexe sur math.SE

Avant d'essayer des suggestions spécifiques, quelle est une plage typique de valeurs pour ?$b$

Notons tout d'abord que qui, si était continu, serait lié à une distribution exponentielle. Ensuite, ce que vous pouvez faire est de simuler à partir d'une distribution exponentielle tronquée et de prendre la (partie entière) des observations.$P(x)\propto e^{-bx}$$x$floor()

Le cdf d'une exponentielle tronquée est

Ensuite, si nous faisons , nous obtenons que . Si est grand, alors qui suggère d'approximer .$F(x;n,b)=u$$x=-\dfrac{1}{b}\log[1-u(1-e^{-bn})]$$bn$$e^{-bn}\approx 0$$x\approx -\dfrac{1}{b}\log[1-u]$

rweirdp <- function(ns,n,b){
u <- runif(ns)
samp <- - log(1-u*(1-exp(-n*b)))/b
return(floor(samp))
}

rweirdp(1000,10,1)

Un moyen d'échantillonner à partir de la distribution cible est de$p(k)\propto \exp\{-bk\}$

1. exécuter une expérience Metropolis-Hastings pour déterminer le support (intéressant) de la distribution, c'est-à-dire dans quel sous-ensemble de il se concentre;$\{0,1,\ldots,n\}$

   metro=function(N,b,n){
   x=sample(0:n,N,rep=TRUE)
   for (t in 2:N){
     x[t]=prop=x[t-1]+sample(c(-1,1),1)
   
     if ((prop<0)||(prop>n)||(log(runif(1))>b*(x[t]-prop)))
         x[t]=x[t-1]
     }
   return(x)
   }
   

2. Utilisez le support ainsi déterminé, disons, pour calculer les probabilités exactes comme pour éviter les débordements.$\{k_0,\ldots,k_1\}$$p(k)\propto \exp\{-bk+bk_0\}$

Mise à jour: En y réfléchissant davantage, puisque diminue en k, le support effectif de la distribution commencera toujours à . Si est assez grand, ce support se terminera très rapidement, auquel cas importe peu car les grandes valeurs de ne seront jamais visitées. Si est très petit, le pdf est presque plat, ce qui signifie que l'on peut utiliser une distribution uniforme sur comme proposition d'acceptation-rejet. Et utilisez les journaux dans l'étape d'acceptation pour éviter les débordements.$p(\cdot)$$k_0=0$$b$$n$$k$$b$$\{0,1,\ldots,n\}$