Gerador de números aleatórios que produz uma distribuição de poder-lei?
Pergunta
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!
Solução
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
Outras dicas
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 ...
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)
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