Jak efektywnie generować posortowane równomiernie rozłożone wartości w przedziale?

Powiedzmy, że chcę wygenerować zestaw liczb losowych z przedziału (a, b). Wygenerowana sekwencja powinna również mieć właściwość, która jest posortowana. Mogę wymyślić dwa sposoby na osiągnięcie tego.

Niech nbędzie długością sekwencji, która ma zostać wygenerowana.

Pierwszy algorytm:

Let `offset = floor((b - a) / n)`
for i = 1 up to n:
   generate a random number r_i from (a, a+offset)
   a = a + offset
   add r_i to the sequence r

Drugi algorytm:

for i = 1 up to n:
    generate a random number s_i from (a, b)
    add s_i to the sequence s
sort(r)

Moje pytanie brzmi: czy algorytm 1 wytwarza sekwencje równie dobre, jak te generowane przez algorytm 2?

Pierwszy algorytm zawodzi źle z dwóch powodów:

1. Zabranie głosu może drastycznie go zmniejszyć. Rzeczywiście, gdy b - a < n , wyniesie zero, dając ci zestaw, którego wartości są takie same!$(a-b)/n$$b-a \lt n$

2. Gdy nie zabierasz głosu, wynikowe wartości są zbyt równomiernie rozłożone. Na przykład, w dowolnym prostym losowej próbie IID jednolite zmiennych towarzyszących (na przykład od a = 0 i b = 1 ), to jest ( 1 - 1 / n ), n ≈ 1 / e ≈ 37 % szansę, że największy nie będzie w górnym przedziale od 1 - 1 / n do 1 . W przypadku algorytmu 1 istnieje 100 $n$$a=0$$b=1$$(1-1/n)^n\approx 1/e\approx 37\%$$1-1/n$$1$$100\%$szansa, że ​​maksimum będzie w tym przedziale. Dla niektórych celów ta superjednorodność jest dobra, ale ogólnie jest to straszny błąd, ponieważ (a) wiele statystyk zostanie zniszczonych, ale (b) ustalenie przyczyny może być bardzo trudne.

3. Jeśli chcesz uniknąć sortowania, zamiast tego generuj niezależne zmienne wykładniczo rozłożone. Normalizuj ich sumę skumulowaną do zakresu ( 0 , 1 ) , dzieląc przez sumę. Usuń największą wartość (która zawsze będzie wynosić 1 ). Przeskaluj do zakresu ( a , b ) .$n+1$$(0,1)$$1$$(a,b)$

Pokazane są histogramy wszystkich trzech algorytmów. (Każda przedstawia skumulowane wyniki niezależnych zestawów wartości n = 100 każda). Brak widocznej zmienności histogramu dla algorytmu 1 wskazuje na problem. Różnice w pozostałych dwóch algorytmach są dokładnie tym, czego można się spodziewać - i tym, czego potrzebujesz od generatora liczb losowych.$1000$$n=100$

Aby poznać wiele innych (zabawnych) sposobów symulowania niezależnych zmiennych jednolitych, zobacz Symulowanie losowań z rozkładu jednolitego przy użyciu losowań z rozkładu normalnego .

Oto Rkod, który utworzył figurę.

b <- 1
a <- 0
n <- 100
n.iter <- 1e3

offset <- (b-a)/n
as <- seq(a, by=offset, length.out=n)
sim.1 <- matrix(runif(n.iter*n, as, as+offset), nrow=n)
sim.2 <- apply(matrix(runif(n.iter*n, a, b), nrow=n), 2, sort)
sim.3 <- apply(matrix(rexp(n.iter*(n+1)), nrow=n+1), 2, function(x) {
  a + (b-a) * cumsum(x)[-(n+1)] / sum(x)
})

par(mfrow=c(1,3))
hist(sim.1, main="Algorithm 1")
hist(sim.2, main="Algorithm 2")
hist(sim.3, main="Exponential")

Pierwszy algorytm wytwarza zbyt równomiernie rozmieszczone liczby

Zobacz także serie o niskiej rozbieżności .

$[0;1]$

(Jak wskazano, może to być pożądana właściwość np. Do stratyfikacji. Szeregi o niskiej rozbieżności, takie jak Halton i Sobel , mają przypadki użycia.)

Właściwe, ale drogie podejście (dla prawdziwych wartości)

... ma używać losowych liczb dystrybuowanych w wersji beta. Statystyka kolejności szeregów rozkładu jednolitego jest rozkładem beta. Możesz użyć tego do losowego narysowania najmniejszego , a następnie drugiego najmniejszego, ... powtórz.

$[0;1]$$\text{Beta}[1,n]$$n$$1-X\sim \text{Beta}[n, 1]$$-\ln (1-X)\sim \text{Exponential}[n]$$\frac{-\ln(U[0;1])}{n}$

Co daje następujący algorytm:

x = a
for i in range(n, 0, -1):
    x += (b-x) * (1 - pow(rand(), 1. / i))
    result.append(x) 

Mogą występować niestabilności numeryczne, a obliczenia powi podział dla każdego obiektu mogą okazać się wolniejsze niż sortowanie.

W przypadku wartości całkowitych konieczne może być zastosowanie innego rozkładu.

Sortowanie jest niezwykle tanie, więc po prostu z niego korzystaj

$O(n \log n)$

Zależy to również od tego, co robisz z liczbami losowymi. W przypadku problemów z integracją numeryczną metoda pierwsza (skorygowana przez usunięcie operatora podłogi) dałaby lepszy zestaw punktów. To, co robisz, jest formą warstwowego próbkowania i ma tę zaletę, że pozwala uniknąć zlepiania. nie można na przykład uzyskać wszystkich wartości z zakresu 0– (ba) / n. To powiedziawszy dla innych aplikacji może to być bardzo złe, zależy to od tego, co chcesz z tym zrobić.