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

12

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?

ultrajohn
źródło
BTW jest niezwykle łatwe do wygenerowania listy posortowanych liczb losowych w R. W celu wytworzenia tablicy zestawów n liczb losowych nad równomiernie rozmieszczonych [ , b ] następujący kod działania: . kn[za,b]rand_array <- replicate(k, sort(runif(n, a, b))
RobertF

Odpowiedzi:

18

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!(za-b)/nb-za<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 ), n1 / 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 %nza=0b=1(1-1/n)n1/mi37%1-1/n1100%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(za,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.1000n=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 .

Rysunek: histogramy

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")
Whuber
źródło
Co sądzisz o algorytmie (opartym na statystykach kolejności rang) w mojej odpowiedzi? ;-)
Ma ZAKOŃCZENIE - Anony-Mousse
@Anony Jest to mniej wydajna wersja mojego algorytmu 3. (Wydaje się, że wiąże się to z wieloma niepotrzebnymi przeskalowaniami.) Generuje się wykładnicze zmienne, biorąc logi mundurów, co jest standardem.
whuber
6

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]Beta[1,n]n1-XBeta[n,1]-ln(1-X)Wykładniczy[n]-ln(U[0;1])n

-ln(1-x)=-ln(1-u)n1-x=u1nx=1-u1n

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(nlogn)

Ma ZAKOŃCZENIE - Anony-Mus
źródło
1
Mogą istnieć powody, aby unikać sortowania. Jednym z nich jest generowanie ogromnej liczby zmiennych losowych, tak wielu, że standardowa procedura sortowania nie jest w stanie ich obsłużyć.
whuber
Myślę, że problemy numeryczne z sumami wykorzystującymi matematykę zmiennoprzecinkową stały się problemem znacznie wcześniej. (I problemy z cyklicznymi wzorami w liczbach pseudolosowych!) Dość łatwo skalować podejście sortowania do terabajtów i eksabajtów w systemach rozproszonych.
Ma ZAKOŃCZENIE - Anony-Mousse
1012
Ok, brak konieczności ich przechowywania to argument. Ale wtedy będziesz potrzebować mojego podejścia, twój wariant 3 wykorzystujący łączną sumę nie zadziała.
Ma ZAKOŃCZENIE - Anony-Mousse
To jest doskonała uwaga. Teraz widzę zalety dodatkowych obliczeń! (+1)
whuber
5

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ć.

użytkownik67054
źródło
2
+1 Myślę, że jest to przydatny wkład do pytania, zwłaszcza poprzez scharakteryzowanie algorytmu 1 pod względem stratyfikacji.
whuber