Generatore di numeri casuali che produce una distribuzione di legge di potenza?

Question

Sto scrivendo alcuni test per la linea di comando di un C ++ Linux app. Mi piacerebbe per generare un po 'di numeri interi con una distribuzione a coda lunga legge di potenza /. Significato, ho un qualche numero molto frequentemente, ma la maggior parte di loro relativamente di rado.

Idealmente ci sarebbe solo alcune equazioni magia ho potuto usare con rand () o una delle funzioni stdlib casuali. In caso contrario, un facile da usare pezzo di C / C ++ sarebbe grande.

Grazie!

Answer 1

Questa pagina rel="noreferrer"> href="http://mathworld.wolfram.com/RandomNumber.html" illustra come ottenere una distribuzione di legge di potenza da una distribuzione uniforme (che è ciò che la maggior parte dei generatori di numeri casuali forniscono).

La risposta breve (di derivazione al link qui sopra):

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

dove y è una variate uniforme, n è il potere di distribuzione, x0 e x1 definiscono l'intervallo della distribuzione, e x è variata la tua legge di potenza distribuita.

Answer 2

Se si conosce la distribuzione che si desidera (chiamato la funzione di probabilità di distribuzione (PDF)) e lo hanno correttamente normalizzata, è possibile integrare per ottenere la funzione di distribuzione cumulativa (CDF), quindi invertire la CDF (se possibile) per ottenere la trasformazione è necessario dalla distribuzione uniforme [0,1] al tuo desiderato.

Così si inizia definendo la distribuzione che si desidera.

P = F(x)

(per x in [0,1]) poi integrato per dare

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

Se questo può essere invertita si ottiene

y = F^{-1}(C)

Così chiamata rand() e collegare il risultato come C nell'ultima riga e utilizzare y.

Questo risultato è chiamato il Teorema Fondamentale del campionamento. Questa è una seccatura a causa del requisito di normalizzazione e la necessità di invertire analiticamente la funzione.

In alternativa è possibile utilizzare una tecnica di rifiuto: gettare un numero uniformemente nell'intervallo desiderato, poi gettare un altro numero e confrontarlo al PDF nella posizione indeicated dal primo tiro. Rifiuta se il secondo tiro supera il PDF. Tende ad essere inefficiente per i PDF con un sacco di regione a bassa probabilità, come quelli con lunghe code ...

Un approccio intermedio comporta invertendo la CDF con la forza bruta: si memorizza il CDF come una tabella di ricerca, e fare una ricerca inversa per ottenere il risultato

.

---

La vera stinker qui è che semplici distribuzioni x^-n non sono normalizzabili sulla gamma [0,1], quindi non è possibile utilizzare il teorema del campionamento. Prova (x + 1) ^ - n invece ...

Answer 3

Non posso commentare la matematica necessaria per produrre una distribuzione di legge di potenza (gli altri posti hanno delle proposte), ma vorrei suggerire di familiarizzare con le TR1 C ++ Standard Library strutture di numeri casuali in <random>. Questi forniscono più funzionalità rispetto std::rand e std::srand. Il nuovo sistema specifica un'API modulare per generatori, motori e distribuzioni e fornisce un gruppo di preset.

I preset di distribuzione inclusi sono:

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

Quando si definisce la vostra distribuzione legge di potenza, si dovrebbe essere in grado di collegarlo con generatori e motori esistenti. Il libro Il C ++ standard libreria di estensioni da Pete Becker ha un grande capitolo sulla <random>.

Ecco un articolo su come creare altre distribuzioni (con esempi di Cauchy, Chi -squared, t di Student e Snedecor F)

Answer 4

Volevo solo di effettuare una simulazione reale come complemento alla risposta (giustamente) ha accettato. Anche se in R, il codice è così semplice da essere (pseudo) Pseudo-codice.

Una piccola differenza tra il Wolfram MathWorld formula nella risposta accettata e altri, forse più comuni, equazioni è il fatto che la legge esponente potenza n (che è tipicamente indicata come alfa) non porta un segno negativo esplicito. Quindi il valore alfa prescelto deve essere negativo, e tipicamente tra 2 e 3.

x0 e x1 indicano i limiti inferiore e superiore della distribuzione.

Così qui è:

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)

o tracciati in scala logaritmica:

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")

Ecco il riepilogo dei dati:

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