générateur de nombres aléatoires qui produit une distribution en loi de puissance?

Question

J'écris quelques tests pour une ligne de commande C de l'application Linux. Je voudrais créer un groupe d'entiers avec une loi de puissance / de distribution à longue queue. Ce qui signifie, je reçois un quelques chiffres très souvent, mais la plupart d'entre eux relativement rarement.

Idéalement, il serait juste quelques équations magiques que je pourrais utiliser avec rand () ou l'une des fonctions aléatoires stdlib. Dans le cas contraire, un outil facile à utiliser morceau de C / C ++ serait grand.

Merci!

Answer 1

à Wolfram MathWorld explique comment obtenir une distribution en loi de puissance d'une distribution uniforme (qui est ce que la plupart des générateurs de nombres aléatoires fournissent).

La réponse courte (dérivation sur le lien ci-dessus):

x = [(x1^(n+1) - x0^(n+1))*y + x0^(n+1)]^(1/(n+1))

où y est une variable aléatoire uniforme, n est la puissance de distribution, x0 et x1 définir la plage de la distribution, et x est votre loi de puissance distribuée variate.

Answer 2

Si vous connaissez la distribution que vous voulez (appelée la fonction Probabilité de distribution (PDF)) et ont correctement normalisé, vous pouvez l'intégrer pour obtenir la distribution cumulative Fonction (CDF), puis inverser la CDF (si possible) pour obtenir la transformation dont vous avez besoin de la distribution de [0,1] uniforme à votre choix.

Vous commencez par définir la distribution que vous voulez.

P = F(x)

(pour x dans [0,1]), puis intégrée à envoyer

C(y) = \int_0^y F(x) dx

Si cela peut être inversé vous obtenez

y = F^{-1}(C)

Alors, appelez rand() et branchez le résultat en tant que C dans la dernière ligne et y utiliser.

Ce résultat est appelé le théorème fondamental de l'échantillonnage. C'est embêtant en raison de l'exigence de normalisation et la nécessité d'inverser analytiquement la fonction.

Alternativement, vous pouvez utiliser une technique de rejet: jeter un nombre uniforme dans la plage souhaitée, puis jeter un autre numéro et comparer au format PDF à l'emplacement indeicated par votre premier jet. Rejeter si le second jet dépasse le PDF. Tend à être inefficace pour les fichiers PDF avec beaucoup de région à faible probabilité, comme ceux qui ont des longues queues ...

Une approche intermédiaire consiste à inverser la CDF par la force brute: vous stockez la CDF comme une table de consultation, et faire une recherche inversée pour obtenir le résultat

.

---

Le vrai salaud est ici que les distributions de x^-n simples sont non normalisable sur la plage [0,1], de sorte que vous ne pouvez pas utiliser le théorème d'échantillonnage. Essayez (x + 1) ^ - n au lieu ...

Answer 3

Je ne peux pas commenter sur le calcul nécessaire pour produire une distribution de loi de puissance (les autres postes ont des suggestions) mais je vous suggère de vous familiariser avec les TR1 C ++ bibliothèque standard installations de nombres aléatoires dans <random>. Ceux-ci fournissent plus de fonctionnalités que std::rand et std::srand. Le nouveau système spécifie une API modulaire pour les générateurs, les moteurs et les distributions et fournit un tas de paramètres prédéfinis.

Les paramètres prédéfinis de distribution inclus sont:

• uniform_int
• bernoulli_distribution
• geometric_distribution
• poisson_distribution
• binomial_distribution
• uniform_real
• exponential_distribution
• normal_distribution
• gamma_distribution

Lorsque vous définissez votre distribution de loi de puissance, vous devriez être en mesure de le brancher avec des générateurs et moteurs existants. Le livre Les extensions bibliothèque standard C ++ par Pete Becker a un grand chapitre sur <random>.

Voici un article sur la façon de créer d'autres distributions (avec des exemples pour Cauchy, Chi -squared, étudiant t et Snedecor F)

Answer 4

Je voulais juste réaliser une simulation réelle en complément à la (à juste titre) accepté réponse. Bien que dans R, le code est aussi simple que d'être (pseudo) -pseudo code.

Une petite différence entre le Wolfram MathWorld formule dans la réponse acceptée et d'autres, peut-être plus commun, les équations est le fait que la exposant de la loi de puissance n (qui est généralement désignée par alpha) ne porte pas un signe négatif explicite. Ainsi, la valeur alpha choisi doit être négative, et typiquement compris entre 2 et 3.

x0 et x1 représentent les limites inférieure et supérieure de la distribution.

Alors la voici:

x1 = 5           # Maximum value
x0 = 0.1         # It can't be zero; otherwise X^0^(neg) is 1/0.
alpha = -2.5     # It has to be negative.
y = runif(1e5)   # Number of samples
x = ((x1^(alpha+1) - x0^(alpha+1))*y + x0^(alpha+1))^(1/(alpha+1))
hist(x, prob = T, breaks=40, ylim=c(0,10), xlim=c(0,1.2), border=F, 
col="yellowgreen", main="Power law density")
lines(density(x), col="chocolate", lwd=1)
lines(density(x, adjust=2), lty="dotted", col="darkblue", lwd=2)

ou tracé en échelle logarithmique:

h = hist(x, prob=T, breaks=40, plot=F)
     plot(h$count, log="xy", type='l', lwd=1, lend=2, 
     xlab="", ylab="", main="Density in logarithmic scale")

Voici le résumé des données:

> summary(x)
   Min.   1st Qu.  Median    Mean   3rd Qu.    Max. 
  0.1000  0.1208  0.1584    0.2590  0.2511   4.9388