Generator liczb losowych, który wytwarza dystrybucję energii?

Question

Piszę niektóre testy dla aplikacji Linuksowej linii poleceń C++. Chciałbym wygenerować paczkę liczb całkowitych z rozkładem mocy/długiego ogona. Znaczy, dostaję kilka liczb bardzo często, ale większość z nich stosunkowo rzadko.

Idealnie byłoby tylko kilka równań magicznych, których mogę użyć z rand() lub jedną z funkcji losowych stdlib. Jeśli nie, to łatwy w użyciu fragment C/C++ byłby świetny.

Dzięki!

Answer 1

Ten numer page at Wolfram MathWorld omawia, w jaki sposób uzyskać dystrybucję energii z rozkładu jednolitego (co zapewnia większość generatorów liczb losowych).

Krótki odpowiedź (pochodzenie z wyżej odnośnik)

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

gdzie Y jest jednolity variate, n jest moc dystrybucyjnej x0 i x1 określenia zakresu dystrybucja, i jest twoją rozproszoną zmiennością prawa.

Answer 2

Jeśli znasz pożądany rozkład (zwany funkcją rozkładu prawdopodobieństwa (PDF)) i masz go poprawnie znormalizowany, możesz go zintegrować, aby uzyskać funkcję skumulowanej dystrybucji (CDF), a następnie zamienić CDF (jeśli to możliwe) na uzyskać transformację, której potrzebujesz od jednolitej dystrybucji [0,1] do pożądanego.

Zacznij od zdefiniowania pożądanego rozkładu.

P = F(x) 

(dla X w [0,1]), a następnie włączone do nadania

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

Jeśli ta może zostać odwrócona masz

y = F^{-1}(C) 

Więc zadzwoń rand() i podłącz wynik w jako C w ostatniej linii i użyj y.

Ten wynik nazywany jest podstawowym twierdzeniem próbkowania. Jest to problem związany z wymogiem normalizacji i potrzebą analitycznego odwrócenia funkcji.

Można również zastosować technikę odrzucania: wyrzuć liczbę równomiernie w żądanym zakresie, a następnie wyrzuć inną liczbę i porównaj z plikiem PDF w lokalizacji wskazanej przez pierwszy rzut. Odrzuć, jeśli drugi rzut przekroczy PDF. Wydaje się być nieefektywny w przypadku plików PDF z obszarem o niskim prawdopodobieństwie, takich jak te z długimi ogonami ...

Podejście pośrednie polega na odwróceniu CDF przez brutalną siłą: przechowuje się CDF jako tabelę odnośników i wykonuje odwrotność lookup, aby uzyskać wynik.

---

Prawdziwy gnojek jest to, że proste x^-n dystrybucje są dla normalizable od zakresu [0,1], więc nie można używać twierdzenie próbkowania. Spróbuj (x + 1)^- n zamiast ...

Answer 3

Nie mogę wypowiedzieć się na temat matematyki wymaganej do stworzenia rozdziału prawa energetycznego (inne posty mają sugestie), ale proponuję zapoznać się z obiektami liczb losowych Standard Library TR1 C++ w <random>. Zapewniają one większą funkcjonalność niż std::rand i std::srand. Nowy system określa modułowe API dla generatorów, silników i dystrybucji oraz dostarcza kilka presetów.

Zawarte ustawień dystrybucyjne:

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

Podczas definiowania swoją dystrybucję mocy prawa, powinna być w stanie podłączyć go do istniejących generatorów i silników. Książka Rozszerzenia biblioteki standardowej C++ autorstwa Pete Becker ma świetny rozdział na temat <random>.

Here is an article o tym, jak tworzyć inne dystrybucje (z przykładami dla Cauchy'ego, chi-kwadrat, t Studenta i Snedecora F)

Answer 4

Chciałem przeprowadzić rzeczywistą symulację jako uzupełnienie (słusznie) zaakceptowanej odpowiedzi . Chociaż w R, kod jest tak prosty, że jest (pseudo) -pseudokodem.

Jedna malutka różnica między Wolfram MathWorld formula w przyjętym odpowiedź i inne, być może bardziej powszechne, równania jest fakt, że prawo wykładnik mocn (która jest zazwyczaj oznaczona jako alfa) nie niesie wyraźny znak ujemny. Dlatego wybrana wartość alfa musi być ujemna, a zwykle pomiędzy 2 a 3.

x0 i x1 oznaczają dolną i górną granicę rozkładu.

Więc to jest tutaj:

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) 

lub wykreślono w skali logarytmicznej:

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

Oto podsumowanie danych:

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