Liczby losowe w R: rnorm, runif, sample i set.seed

Najpierw set.seed()

Losowe wartości, które zmieniają się przy każdym uruchomieniu, są bezużyteczne w nauczaniu, ocenianiu, debugowaniu i nauce: jeśli twoja symulacja dziś daje 0.146, a jutro 0.153, która liczba trafi do raportu? set.seed() ustala punkt startowy generatora, dzięki czemu cały następujący po nim strumień „losowości” da się dokładnie powtórzyć:

set.seed(42)
print(rnorm(3))

set.seed(42)          # same seed...
print(rnorm(3))       # ...same "random" numbers

set.seed(7)           # different seed, different stream
print(rnorm(3))

To samo ziarno daje identyczne wartości, zawsze i na każdym komputerze. Liczby są pseudolosowe: to deterministyczny ciąg zaprojektowany tak, żeby przechodził każdy statystyczny test losowości, a ziarno wybiera, od którego miejsca ciągu zaczynasz. Sama wartość ziarna nic nie znaczy: 42, 7, 20260807, wybierz cokolwiek, byle to zapisać. Nawyk do wyrobienia: jedno set.seed() na początku każdego skryptu, który korzysta z losowości. (Pamiętaj, że każde losowanie przesuwa stan generatora, więc dla powtarzalności liczy się też kolejność wywołań).

System d/p/q/r: jedna tabela rozszyfrowuje całą bibliotekę

R nazywa każdą funkcję rozkładu według schematu przedrostek + rodzina. Gdy zobaczysz tę siatkę, przeczytasz całą bibliotekę statystyczną:

p i q są wzajemnie odwrotne i to właśnie tę parę spotkasz w artykule o przedziałach ufności jako qt(0.975, df):

print(pnorm(1.96))      # P(Z <= 1.96), about 0.975
print(qnorm(0.975))     # which z has 97.5% below it? About 1.96
print(dnorm(0))         # height of the bell curve at its peak

Zmień nazwę rodziny, a wszystko działa tak samo: runif/dunif/punif/qunif, rbinom/..., rpois/..., rt/..., rexp/.... Nauczysz się czterech przedrostków i masz dziesiątki rozkładów.

rnorm(): losowanie z rozkładu normalnego

rnorm(n, mean, sd) losuje z krzywej dzwonowej. To podstawowe narzędzie do symulowania danych przypominających pomiary:

set.seed(42)
iq <- rnorm(1000, mean = 100, sd = 15)

print(mean(iq))
print(sd(iq))
print(mean(iq > 130))   # share more than 2 sd above the mean

Średnia i odchylenie standardowe z próby wychodzą blisko 100 i 15, ale nie dokładnie: ta różnica to właśnie szum próbkowania, który maleje wraz ze wzrostem n. Ostatnia linia to sztuczka, którą warto przejąć: mean() z wektora logicznego to odsetek wartości TRUE, a udział powyżej 130 wychodzi blisko teoretycznego 1 - pnorm(130, 100, 15) ≈ 2.3%. Wartości domyślne to mean = 0, sd = 1 (standardowy rozkład normalny). Histogram dla iq pokazuje znajomy dzwon.

runif(), rbinom(), rpois()

Trzy kolejne rodziny pokrywają większość potrzeb symulacji:

set.seed(42)

print(runif(5))                    # uniform on [0, 1]
print(runif(3, min = 10, max = 20))

print(rbinom(8, size = 10, prob = 0.5))   # successes in 10 coin flips, 8 times

print(rpois(8, lambda = 3))        # event counts averaging 3 per interval

runif() rozkłada wartości równomiernie w przedziale („uniform”, a nie „run if”, każdy raz to źle przeczyta). rbinom(n, size, prob) symuluje n eksperymentów po size prób każdy i zwraca liczbę sukcesów w każdym, więc każda wartość powyżej to liczba orłów w 10 rzutach. rpois(n, lambda) generuje liczby zdarzeń, które występują niezależnie ze znaną średnią częstością: zgłoszenia do supportu na godzinę, literówki na stronę.

sample(): próbkowanie i tasowanie

Funkcje r* wymyślają wartości z rozkładu, a sample() losuje z wartości, które już masz:

set.seed(42)

print(sample(1:10, 3))                              # 3 values, no repeats
print(sample(c("win", "lose"), 5, replace = TRUE))  # 5 draws, repeats allowed
print(sample(1:5))                                  # no size = shuffle everything

Cała sprawa sprowadza się do argumentu replace: FALSE (domyślnie) to rozdawanie kart, każda wartość może pojawić się raz, a prośba o więcej, niż masz, kończy się błędem; TRUE to rzucanie kostką, każde losowanie zaczyna od nowa. Losowanie ze zwracaniem z własnych danych napędza bootstrap. Wywołana z samym wektorem, sample(x) zwraca losową permutację: to idiom do tasowania.

Losowanie wierszy ramki danych wykorzystuje sample() w indeksowaniu wierszy:

set.seed(42)

random_cars <- mtcars[sample(nrow(mtcars), 5), ]
print(random_cars[, c("mpg", "wt", "hp")])

sample(nrow(mtcars), 5) wybiera 5 losowych numerów wierszy, a indeksowanie wyciąga te wiersze. To standardowy sposób na wyrywkowe sprawdzenie dużego zbioru danych albo podział na zbiór treningowy i testowy.

Mała symulacja Monte Carlo

Oto korzyść z całego zestawu narzędzi. Pytanie: proces daje pomiary o rozkładzie N(100, 15); uśredniasz 10 z nich. Jakie jest prawdopodobieństwo, że średnia przekroczy 105? Zamiast wyprowadzać odpowiedź, zasymuluj ją: przeprowadź eksperyment 10 000 razy i policz:

set.seed(42)

trial <- function() {
  mean(rnorm(10, mean = 100, sd = 15)) > 105
}

results <- replicate(10000, trial())
print(mean(results))

# The exact answer, for comparison:
print(1 - pnorm(105, mean = 100, sd = 15 / sqrt(10)))

Symulacja trafia w kilka tysięcznych od dokładnej wartości (około 0.146). To Monte Carlo w pigułce: zapisz jedną próbę jako funkcję, powtórz ją tysiące razy przez replicate(), weź mean() z sukcesów. Tutaj dokładna odpowiedź istniała, bo układ był podręcznikowo prosty. Gdy tylko pytanie się komplikuje (dziwne rozkłady, maksimum skorelowanych wartości, gra oparta na regułach), droga analityczna się zamyka, a przepis na symulację działa bez zmian. Zwróć uwagę na 15 / sqrt(10) w sprawdzeniu: średnie wahają się mniej niż pojedyncze wartości. To to samo prawo pierwiastka, które decyduje o szerokości przedziału ufności.

Do symulacji, nie do sekretów

Jedna granica, której trzeba przestrzegać: domyślny generator R (Mersenne Twister) jest zbudowany z myślą o jakości statystycznej i szybkości, a nie o tajności. Jego wynik jest deterministyczny przy danym ziarnie, a stan wewnętrzny da się odtworzyć z obserwowanych wartości. To dyskwalifikuje go przy hasłach, tokenach i wszystkim, co dotyka bezpieczeństwa. Do symulacji, bootstrapu i nauczania sprawdza się świetnie; do kryptografii użyj specjalnej biblioteki (np. pakietu openssl), nigdy sample() ani runif().

Najważniejsze informacje

• Jedno set.seed() na początku sprawia, że każdy „losowy” wynik jest powtarzalny: to samo ziarno, te same wartości.
• Przedrostki d/p/q/r rozszyfrowują całą bibliotekę rozkładów: losowanie, gęstość, prawdopodobieństwo skumulowane, kwantyl.
• rnorm(n, mean, sd), runif(n, min, max), rbinom(n, size, prob), rpois(n, lambda) pokrywają większość potrzeb symulacji.
• sample() losuje z twoich własnych wartości: replace = TRUE dla kostki, domyślnie dla kart, bez size do tasowania; df[sample(nrow(df), k), ] losuje wiersze.
• Monte Carlo = funkcja jednej próby + replicate() + mean(): przepis, który odpowiada na pytania o prawdopodobieństwo, do których matematyka nie dociera wygodnie.
• Generator liczb losowych R służy do symulacji, nie do kryptografii.

Dalej: debugowanie, czyli co naprawdę znaczą komunikaty o błędach w R i jak czytać traceback.