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

Question

J'écris quelques tests pour une application Linux en ligne de commande C++. Je voudrais générer un tas d'entiers avec une loi de puissance/longue queue. C'est-à-dire, je reçois quelques chiffres très fréquemment mais la plupart d'entre eux relativement rarement.

Idéalement, il y aurait juste quelques équations magiques que je pourrais utiliser avec rand() ou l'une des fonctions aléatoires de stdlib. Sinon, un morceau facile à utiliser de C/C++ serait génial.

Merci!

Answer 1

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

La réponse courte (dérivation à l'adresse 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 variable distribuée de loi de puissance.

Answer 2

Si vous connaissez la distribution souhaitée (appelée fonction de distribution de probabilité (PDF)) et que vous la normalisez correctement, vous pouvez l'intégrer pour obtenir la fonction de distribution cumulative (CDF), puis inverser le CDF (si possible) Obtenez la transformation dont vous avez besoin de la distribution uniforme [0,1] à votre désiré.

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

P = F(x) 

(pour x dans [0,1]), puis intégré pour donner

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

Si cela peut être inversée vous obtenez

y = F^{-1}(C) 

Alors, appelez rand() et branchez le résultat dans comme C dans la dernière ligne et utilisez y.

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

Alternativement, vous pouvez utiliser une technique de rejet: lancer un nombre uniformément dans la gamme désirée, puis jeter un autre nombre et comparer au PDF à l'emplacement indeicated par votre premier lancer. Rejeter si le second lancer dépasse le PDF. Tend à être inefficace pour les PDF avec beaucoup de région de faible probabilité, comme ceux avec des queues longues ...

Une approche intermédiaire consiste à inverser le CDF par force brute: vous stockez le CDF comme une table de recherche, et faites un reverse recherche pour obtenir le résultat.

---

Le vrai salaud ici est aussi simple que cela x^-n distributions 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 à la place ...

Answer 3

Je ne peux pas commenter les maths requises pour produire une distribution de loi de puissance (les autres messages ont des suggestions) mais je vous suggérerais de vous familiariser avec les facilités de numérotation aléatoire de bibliothèque standard TR1 C++ 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 ensemble de préréglages.

Les présélections de distribution sont inclus:

• 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 de bibliothèque standard C++ par Pete Becker a un grand chapitre sur <random>.

Here is an article sur la façon de créer d'autres distributions (avec des exemples pour Cauchy, Chi-carré, é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 soit si simple qu'il soit pseudo-pseudo-code.

Une petite différence entre le Wolfram MathWorld formula dans la réponse acceptée et d'autres, peut-être plus fréquentes, des équations est le fait que l'exposant de loi de puissance n (qui est généralement désigné comme alpha) ne porte pas de signe négatif explicite. Donc la valeur alpha choisie doit être négative, et typiquement entre 2 et 3.

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

Alors 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