Gerador de números aleatórios que produz uma distribuição de poder-lei?

Question

Eu estou escrevendo alguns testes para aplicação Linux linha de comando a C ++. Eu gostaria de gerar um monte de números inteiros com uma distribuição de lei de potência / long-tail. Significado, eu recebo um alguns números muito frequentemente, mas a maioria deles com pouca frequência.

O ideal seria apenas algumas equações mágicas que eu poderia usar com rand () ou uma das funções aleatórias stdlib. Se não, um fácil de usar pedaço de C / C ++ seria ótimo.

Obrigado!

Answer 1

Esta página no discute Wolfram MathWorld como obter uma distribuição de lei de potência a partir de uma distribuição uniforme (que é o que a maioria dos geradores de números aleatórios fornecer).

A resposta curta (derivação no link acima):

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

onde y é uma variável uniforme, n é o poder de distribuição, x0 e x1 definir o intervalo da distribuição, e x variate é o seu poder-lei distribuído.

Answer 2

Se você sabe a distribuição que você quer (chamado a função de probabilidade de distribuição (PDF)) e tê-lo devidamente normalizado, você pode integrá-lo para obter a função de distribuição cumulativa (CDF), então invertido a CDF (se possível) para obter a transformação que você precisa de distribuição [0,1] uniforme para o seu desejado.

Assim que você começar por definir a distribuição que você deseja.

P = F(x)

(para x em [0,1]), em seguida, integrado para dar

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

Se isto pode ser invertido você começa

y = F^{-1}(C)

Assim chamada rand() e ligue o resultado em como C na última linha e uso y.

Este resultado é chamado de teorema fundamental da amostragem. Este é um aborrecimento devido à exigência de normalização ea necessidade de analiticamente invertido a função.

Como alternativa, você pode usar uma técnica de rejeição: atirar um número uniformemente no intervalo desejado, em seguida, jogar outro número e comparar com o PDF no local indeicated pelo seu primeiro lance. Rejeitar se o segundo lance excede o PDF. Tende a ser ineficiente para PDFs com um monte de região de baixa probabilidade, como aqueles com caudas longas ...

Uma abordagem intermediária envolve invertendo a CDF por força bruta:. Você armazenar o CDF como uma tabela de pesquisa, e fazer uma pesquisa inversa para obter o resultado

---

O verdadeiro stinker aqui é que as distribuições x^-n simples são não-normalizable na [0,1] gama, então você não pode usar o teorema da amostragem. Tente (x + 1) ^ - n vez ...

Answer 3

Eu não posso comentar sobre a matemática necessária para produzir uma distribuição de lei de potência (as outras mensagens tiver sugestões), mas eu sugiro que você se familiarizar com os TR1 biblioteca padrão C ++ instalações de números aleatórios em <random>. Estes fornecem mais funcionalidade que std::rand e std::srand. O novo sistema especifica um API modular para os geradores, motores e distribuições e suprimentos um grupo de presets.

Os presets de distribuição incluídos são:

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

Quando você define sua distribuição lei de potência, você deve ser capaz de ligá-lo com geradores e motores existentes. O livro As extensões de biblioteca C ++ padrão por Pete Becker tem um grande capítulo sobre <random>.

Aqui está um artigo sobre como criar outras distribuições (com exemplos de Cauchy, Chi -squared, t de Student e Snedecor F)

Answer 4

Eu só queria realizar uma simulação real como um complemento para o (legitimamente) resposta aceita. Embora em R, o código é tão simples como a ser (pseudo) -pseudo-código.

Uma diferença pequena entre a Wolfram MathWorld fórmula na resposta aceita e outro, talvez mais comum, equações é o fato de que o lei de potência expoente n (que normalmente é indicado como alfa) não carrega um sinal negativo explícito. Assim, o valor de alfa escolhido tem de ser negativo, e, tipicamente, entre 2 e 3.

x0 e estande x1 para os limites inferior e superior da distribuição.

Então aqui está:

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 plotados em escala logarítmica:

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

Aqui está o resumo dos dados:

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