Numeri casuali in R: rnorm, runif, sample e set.seed

Prima di tutto set.seed()

Estrazioni casuali che cambiano a ogni esecuzione sono inutili per insegnare, correggere, fare debug o fare scienza: se la tua simulazione dice 0.146 oggi e 0.153 domani, quale numero va nel report? set.seed() fissa il punto di partenza del generatore e rende esattamente ripetibile tutta la sequenza di "casualità" che segue:

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

Stesso seed, estrazioni identiche, sempre e su ogni macchina. I numeri sono pseudocasuali: una sequenza deterministica costruita per superare ogni test statistico di casualità, in cui il seed sceglie da quale punto della sequenza parti. Il valore del seed in sé non ha alcun significato: 42, 7, 20260807, scegli quello che vuoi, basta annotarlo. L'abitudine da prendere: un set.seed() all'inizio di ogni script che usa la casualità. (Nota che ogni estrazione fa avanzare lo stato, quindi anche l'ordine delle chiamate conta per la riproducibilità.)

Il sistema d/p/q/r: una tabella per decifrare tutta la libreria

R chiama ogni funzione di distribuzione con prefisso + famiglia, e una volta vista la griglia riesci a leggere tutta la libreria statistica:

p e q sono l'una l'inversa dell'altra, e sono la coppia che hai incontrato negli intervalli di confidenza come 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

Cambia il nome della famiglia e tutto resta valido: runif/dunif/punif/qunif, rbinom/..., rpois/..., rt/..., rexp/.... Impari quattro prefissi e ottieni decine di distribuzioni.

rnorm(): estrazioni da una normale

rnorm(n, mean, sd) estrae valori da una curva a campana: è lo strumento principale per simulare dati simili a misurazioni:

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

La media e la deviazione standard del campione finiscono vicino a 100 e 15, non esattamente lì: quello scarto è il rumore di campionamento, che si riduce al crescere di n. L'ultima riga è un trucco da rubare: mean() di un vettore logico è la proporzione di TRUE, e la quota oltre 130 viene vicina al valore teorico 1 - pnorm(130, 100, 15) ≈ 2.3%. I valori di default sono mean = 0, sd = 1 (la normale standard). Un istogramma di iq mostra la classica campana.

runif(), rbinom(), rpois()

Altre tre famiglie coprono la maggior parte delle esigenze di simulazione:

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() distribuisce le estrazioni in modo uniforme su un intervallo ("uniform", non "run if": prima o poi tutti lo leggono male). rbinom(n, size, prob) simula n esperimenti da size prove ciascuno e restituisce il numero di successi di ognuno: ogni valore qui sopra è il numero di teste su 10 lanci. rpois(n, lambda) genera conteggi di eventi che avvengono in modo indipendente a un tasso medio noto: ticket di assistenza all'ora, refusi per pagina.

sample(): campionare e mescolare

Mentre le funzioni r* inventano valori a partire da una distribuzione, sample() estrae da valori che hai già:

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

Tutto ruota attorno all'argomento replace: FALSE (il default) distribuisce carte, ogni valore può comparire una sola volta e chiederne più di quanti ne hai è un errore; TRUE lancia dadi, ogni estrazione riparte da zero. Campionare con reinserimento dai tuoi dati è il motore del bootstrap. Chiamata con il solo vettore, sample(x) restituisce una permutazione casuale: è il modo idiomatico per mescolare.

Per campionare le righe di un data frame si usa sample() dentro l'indicizzazione delle righe:

set.seed(42)

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

sample(nrow(mtcars), 5) sceglie 5 numeri di riga a caso; l'indicizzazione estrae quelle righe. È la mossa standard per controllare a campione un grande dataset o per dividere i dati in training e test.

Una mini simulazione Monte Carlo

Ecco il frutto di tutti questi strumenti. La domanda: un processo produce misurazioni distribuite come N(100, 15); ne fai la media di 10. Qual è la probabilità che la media superi 105? Invece di ricavare la risposta con i calcoli, simulala: ripeti l'esperimento 10.000 volte e conta:

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

La simulazione arriva a pochi millesimi dal valore esatto (circa 0.146). Questo è il Monte Carlo in poche parole: scrivi una prova come funzione, ripetila migliaia di volte con replicate(), calcola mean() dei successi. Qui la risposta esatta esisteva perché il problema era semplice da manuale; appena la domanda si complica (distribuzioni strane, il massimo di estrazioni correlate, un gioco basato su regole), la via analitica si chiude e la ricetta della simulazione continua a funzionare senza modifiche. Nota il 15 / sqrt(10) nella verifica: le medie oscillano meno delle singole estrazioni, la stessa legge della radice quadrata che determina l'ampiezza degli intervalli di confidenza.

Adatto alle simulazioni, non ai segreti

Un limite da rispettare: il generatore predefinito di R (Mersenne Twister) è pensato per la qualità statistica e la velocità, non per la segretezza. Il suo output è deterministico dato il seed e il suo stato interno si può ricostruire osservando l'output: difetti fatali per password, token o qualsiasi cosa legata alla sicurezza. Per simulazioni, bootstrap e didattica è ottimo; per la crittografia usa una libreria apposita (per esempio il pacchetto openssl), mai sample() o runif().

Cosa ti porti a casa

• Un set.seed() all'inizio rende riproducibile ogni risultato "casuale": stesso seed, stesse estrazioni.
• I prefissi d/p/q/r decifrano tutta la libreria delle distribuzioni: estrazione casuale, densità, probabilità cumulata, quantile.
• rnorm(n, mean, sd), runif(n, min, max), rbinom(n, size, prob), rpois(n, lambda) coprono la maggior parte delle esigenze di simulazione.
• sample() estrae dai tuoi valori: replace = TRUE per i dadi, il default per le carte, senza size per mescolare; df[sample(nrow(df), k), ] campiona le righe.
• Monte Carlo = funzione per una prova + replicate() + mean(): la ricetta che risponde a domande di probabilità che la matematica non raggiunge comodamente.
• Il generatore di R serve per le simulazioni, non per la crittografia.

Prossimo passo: il debug, cioè cosa significano davvero i messaggi di errore di R e come leggere un traceback.