momente şi schiţe de informatică şi matematică
To attain knowledge, write. To attain wisdom, rewrite.

Simularea aruncării zarului

Excel | Gnumeric | Octave | limbajul R
2016 jan

Dacă am repeta suficient de mult experienţa aruncării unui aceluiaşi zar, ar trebui să obţinem valori foarte apropiate una de alta, pentru numărul de apariţii ale fiecăreia dintre cele şase feţe ale zarului; altfel, dacă tot experimentând astfel constatăm discrepanţe clare privind frecvenţa apariţiei feţelor - atunci putem conchide că zarul sau generatorul de numere aleatoare implicat este defectuos.

Putem instrumenta astfel de experimente fie în manieră interactivă "standard" cu Excel (sau cu altă aplicaţie de calcul tabelar), fie printr-un program adecvat într-un limbaj de programare sau altul. Întâi, descriem o simulare a aruncării zarului folosind Gnumeric - cu următoarea imagine finală:

Am început prin a introduce în celula A1 formula =randbetween(1, 6), care va selecta aleatoriu un număr întreg din intervalul [1, 6]; valoarea 2 rezultată eventual în celula A1 poate fi interpretată în sensul apariţiei feţei cu 2 puncte, la prima aruncare a zarului.

Apoi, am selectat un număr de coloane (primele cinci, în imaginea noastră) şi am pastat formula din celula A1; în Gnumeric, o coloană are în mod implicit 216 celule - încât ne-a rezultat un număr de 5*216 celule conţinând fiecare punctajul apărut (în mod aleatoriu) la câte o aruncare a zarului.

Mai departe - pentru a obţine frecvenţa de apariţie pentru fiecare faţă a zarului - am selectat celulele G1:G6 şi am tastat formula =frequency(A:E, {1,2,3,4,5,6}) în care primul parametru A:E vizează concis întreaga zonă de celule din cele cinci coloane, iar al doilea parametru {1,2,...,6} este tabloul valorilor a căror frecvenţă ne interesează. Prin combinaţia de taste ctrl şi enter obţinem frecvenţele respective, în celulele G1:G6. Iar reselectând G1:G6 şi alegând din meniul Insert - obţinem histograma corespunzătoare acestor valori, redată (după unele retuşuri de format) în imaginea de mai sus.

Pentru eventuale confruntări de rezultate, pot fi utile unele valori standard: în celula H10 înregistrăm numărul total de repetări ale experienţei =count(A:E); în H11 avem =average(G:G), adică valoarea teoretică a numărului de apariţii pentru fiecare faţă; în H12 calculăm dispersia numărului de apariţii, dată de =stdevp(G:G). Putem adăuga desigur, media punctelor apărute =average(A:E), abaterea standard de la medie =stdevp(A:E), etc.

În mod standard, Gnumeric recalculează toate formulele dacă se apasă tasta F9 (corespunzătoare meniului Edit/Recalculate); deci putem repeta (printr-o singură apăsare de tastă) întreaga simulare descrisă mai sus, de oricâte ori dorim. Redăm o asemenea recalculare (instantanee), în care ne-a rezultat o dispersie ceva mai mare decât în simularea redată în imaginea precedentă:

Desigur că poate să nu fie cazul de a invoca direct şi pe rând, funcţiile menţionate mai sus; de obicei aplicaţiile de calcul tabelar oferă şi instrumente pentru analiza imediată a datelor - sintetizatoare care înlănţuie în mod automat formulele potrivite pentru obţinerea şi redarea informaţiilor statistice tipice, asupra datelor selectate. De exemplu, Google Sheets oferă instrumentul sintetic "Explore":

Google Sheets prevede 1000 de rânduri, pe foaia de calcul; am folosit însă formularul de adăugare a unui număr de rânduri, prezentat după ultimul rând - încât în coloana A am obţinut 52000 de valori aleatoare 1..6 (invocând desigur, funcţia between()); am implicat şi coloana B, dar "Explore" oferă statistici şi o histogramă pentru câte o singură coloană - în imaginea tocmai redată fiind vizate datele din coloana B (cea selectată). De altfel, lucrând în cadrul browser-ului (cum este cazul pentru aplicaţia Google Sheets) - operarea unui volum mare de date (mai ales dacă vrem şi reprezentări grafice) poate deveni mult prea costisitoare în privinţa timpului necesitat.

Mai departe ilustrăm unele experimente de simulare a aruncării zarului folosind GNU Octave:

Am lansat octave cu opţiunea "-q" ("--quiet" suprimă mesajul iniţial). Cu randi(6, 1, 60000) obţinem un vector cu 60000 de numere întregi selectate în mod aleatoriu dintre valorile 1..6; apoi, hist(..., 6) produce o histogramă (cu 6 bare) a valorilor din vectorul respectiv.

La prompt-ul interpretorului ("octave:2>") putem actualiza comanda precedentă folosind "tastele săgeată"; astfel, este foarte comod să repetăm experimentul - observând direct pe histogramă ce variaţii apar. Aceste variaţii pot fi sesizate exact afişând valorile returnate de funcţia hist() - vezi în imaginea de mai sus valorile variabilei seen (instituite prin seen = hist(...)): faţa 1 s-a obţinut de 10102 ori, faţa 2 de 10036 ori, etc. (teoretic, fiecare faţă trebuia să apară de 10000 de ori). Dacă mărim de 10 ori (de exemplu) numărul de aruncări (înlocuind 60000 cu 600000 şi repetând comenzile de mai sus), atunci vom constata că valorile din vectorul seen se nivelează până la ordinul unităţilor.

Am putea organiza mai bine lucrurile, într-un script; prevedem să zicem (doar) N=6000 de aruncări consecutive ale zarului, dar repetăm această serie de aruncări de (să zicem) R=1000 de ori:

clear  % elimină variabilele create anterior
N = 6000;  % de câte ori aruncăm succesiv, zarul
R = 1000;  % de câte ori repetăm aruncarea de câte N ori a zarului

for n = 1:R
    runs = randi(6, 1, N); % aruncă zarul de N de ori consecutiv
    seen(n, :) = hist(runs, 6); % seen(n, F) = numărul de apariţii ale feţei F=1..6
end

for F = 1:6
    avg(F) = mean(seen(:, F)); % media apariţiilor feţei F în cele R repetări
end

disp(avg)

Salvăm comenzile redate mai sus în fişierul "simulare.m" şi lansăm octave din directorul în care am salvat acest script; "comanda" simulare execută secvenţa de comenzi din scriptul nostru:

vb@Home:~/aaaa$ octave -q
octave:1> simulare  % cu R = 1000 repetări, ale aruncării zarului de câte 6000 de ori
   1000.55    999.42    998.98   1000.98   1001.06    999.01
octave:2> simulare
    999.21    999.00    999.32   1001.64   1001.05    999.79
octave:3> 

Am repetat "simulare" (fiind curios, dar fiind şi aşa de simplu: săgeata-sus şi "enter") de vreo 20 de ori; o singură dată am avut sub 998 (de genul 997.51) - marea majoritate a valorilor rezultate fiind între 999 şi 1001 (foarte rar întâlnind 998 şi 1002). Cu alte cuvinte, repetând de R=1000 de ori experienţa aruncării de N=6000 de ori a zarului - aproape că am obţinut rezultatul teoretic: fiecare faţă să apară de câte 1000 de ori; însă cu R=100 de exemplu, avem o variaţie sensibil mai mare (între 995 şi 1005) în jurul mediei teoretice:

octave:40> simulare  % numai R = 100 repetări, ale aruncării zarului de câte 6000 de ori
   1002.45    995.97    997.45   1000.51    997.64   1005.98
octave:41> simulare
    996.22   1001.96   1006.11   1001.65    991.83   1002.23
octave:42>

Mai departe complicăm puţin lucrurile, vizând suma punctelor apărute prin aruncarea a două zaruri; dar acum, pentru simulare vom folosi limbajul R.

Sub R, aruncarea unui zar poate fi simulată folosind sample(x, size, replace):

vb@Home:~/aaaa$ R -q
> sample(1:6, 1, replace=TRUE)
[1] 5  # a apărut faţa cu 5 puncte
> sample(1:6, 2, TRUE)  # selectează aleatoriu două valori 1..6 (cu repunere)
[1] 3 4

replicate(n, expr) repetă de n ori evaluarea expresiei indicate, returnând tabloul rezultatelor. table(x) calculează frecvenţa fiecărei valori din tabloul indicat, iar plot(...) produce o histogramă în fişierul şi cu mărimea indicate în apelul png(...). Constituim cu aceste elemente funcţia make() - într-un fişier "simulare.R" pe care îl marcăm ca executabil (prin comanda chmod +x simulare.R):

#!/usr/bin/Rscript
make <- function(N) {
    runs <- replicate(N, sum(sample(1:6, 2, TRUE)))  # N sume de 2 valori aleatorii 1..6
    seen <- table(runs) / length(runs)  # frecvenţa relativă a valorilor 2..12
    print(round(seen, 3))
    title <- sprintf("%d.png", N)
    png(filename=title, width=340, height=340, pointsize=10) 
    plot(seen, xlab='Suma punctelor', ylab='Frecvenţa relativă', main = title) 
    grid()
}

args <- commandArgs(TRUE)  # ./simulare.R 200 1000 5000 36000 360000
nar <- length(args)
if(nar == 0) {
    N3 = c(3600, 36000)
} else {
    N3 <- as.integer(args[1:nar])
}

for(N in N3) {
    make(N)
}

N-avem decât să inserăm (temporar) în scriptul de mai sus o linie print(commandArgs()), pentru a constata direct care sunt pentru R, argumentele liniei de comandă:

vb@Home:~/aaaa$ ./simulare.R  200  36000  360000
[1] "/usr/lib/R/bin/exec/R" "--slave"               "--no-restore"         
[4] "--file=./simulare.R"   "--args"                "200"                  
[7] "36000"                 "360000"               

În scriptul nostru am apelat commandArgs(TRUE) (cu un argument TRUE) - pentru a obţine numai valorile de după "--args" (tastate pe linia de comandă după numele scriptului). Pentru argumentele respective se afişează pe ecran frecvenţele relative corespunzătoare valorilor posibile ale sumei zarurilor şi se constituie câte un fişier conţinând histograma corespunzătoare:

vb@Home:~/aaaa$ ./simulare.R  200  36000  360000
runs pentru 200 de aruncări
    2     3     4     5     6     7     8     9    10    11    12 
0.020 0.070 0.120 0.120 0.115 0.200 0.130 0.120 0.040 0.060 0.005 
runs pentru 3600 de aruncări
    2     3     4     5     6     7     8     9    10    11    12 
0.027 0.055 0.084 0.113 0.139 0.164 0.140 0.110 0.084 0.056 0.028 
runs pentru 360000 de aruncări
    2     3     4     5     6     7     8     9    10    11    12 
0.028 0.055 0.083 0.111 0.139 0.167 0.138 0.112 0.084 0.055 0.028 

Ştim desigur că - notând cu V o valoare 2..7 - suma zarurilor este egală cu V cu aceeaşi probabilitate cu care este egală cu (12 - V) şi anume, cu probabilitatea (V-1)/36; în R "1:6" dă un vector conţinând valorile 1..6, iar operarea aritmetică a unui vector cu un scalar paralelizează operaţia respectivă pe componentele vectorului - încât putem obţine valorile (V-1)/36, rotunjite la 3 zecimale, astfel:

vb@Home:~/aaaa$ R  -q
> round(1:6 / 36, 3)
[1] 0.028 0.056 0.083 0.111 0.139 0.167

Vedem mai sus (inclusiv pe imaginea "200.png") că pentru un număr mic de aruncări ale zarurilor (N = 200, dar chiar N = 1000), rezultatele diferă sensibil faţă de cele teoretice; în schimb, pentru un număr suficient de mare de aruncări (mai sus, N=300000) diferenţele apar eventual după a treia zecimală sau în orice caz - sunt mai mici decât o sutime (fiind astfel neglijabile).

vezi Cărţile mele (de programare)

docerpro | Prev | Next