Statistica e Calcolo delle probabilità - Sintesi
0. Premessa
1. Approssimazioni, istogrammi ed altre rappresentazioni di dati distribuiti in modalità di tipo
non numerico
2. Rappresentazioni statistiche di dati numerici non classificati. Valori medi
3. Il caso dei dati interi e di quelli già classificati
4. Il calcolo delle probabilità
5. Leggi di distribuzione
6. Il teorema limite centrale. Altre leggi di distribuzione
7. Dipendenza e indipendenza stocastica
8. Sistemi di variabili casuali
9. Esercizi
0. Premessa
In questa scheda riassumiamo brevemente gli argomenti di statistica e probabilità affrontati negli anni precedenti.
1. Approssimazioni, istogrammi ed altre rappresentazioni di dati distribuiti in modalità di tipo non numerico
Consideriamo la tabella (1.1), in cui è riportato quanto si è speso in beni di consumo (alimenti, vestiti, automobili, …) e in servizi (taglio dei capelli, viaggi in treno, …) in Italia in due anni particolari.
| (1.1) |
| ||||||||||||||||||||||||||||||||||||||||
Nella tabella (1.2) abbiamo riportato gli stessi dati ma arrotondati ai miliardi. Ad esempio
77 749 milioni è più vicino a 78 000 milioni, ossia a
78 miliardi, che a 77 000 milioni, ossia a 77 miliardi, quindi viene arrotondato al primo
valore. In generale, per arrotondare un numero a n cifre si guarda la cifra n+1-esima: se questa è minore di 5 si
si sostituiscono con 0 tutte le cifre a destra del posto n, altrimenti si aumenta di uno la cifra di posto n e
si sostituiscono con 0 tutte le cifre alla sua destra. Si dice, anche, che 78 000 è
l'arrotondamento di 77 749 a 2 cifre significative.
Analogamente 951 000 è
l'arrotondamento di 950 502 a 3 cifre significative.
| (1.2) |
| ||||||||||||||||||||||||||||||||||||||||
Invece il troncamento ai miliardi è
77 000 milioni, ossia a 77 miliardi: per troncare un numero a n cifre si
sostituiscono, in ogni caso, con 0 tutte le cifre a destra del posto n. Si dice, anche, che 77 000
è il troncamento di 77 749 a 2 cifre significative.
Analogamente 950 000 è il
troncamento di 950 502 a 3 cifre significative.
Il numero 77 749 milioni, ossia 77 749 000 000, viene scritto più brevemente, e in modo più comprensibile, in notazione scientifica, ossia come 7.7749·1010.
Anche il software scrive i numeri molto grandi o molto piccoli in notazione scientifica. Ecco, ad esempio, come R (in cui, ad es., 5·10³ può essere scritto 5e3) rappresenta 950 502 milioni e 4 milionesimi:
950502e6; 0.000004 #9.50502e+11 4e-06
Vediamo come con R , dato un numero, posso calcolarne il più grande intero minore o eguale ad esso, il suo arrotondamento a tre cifre significative, il suo arrotondamento alla 3ª cifra a destra di quella delle unità e quello alla prima cifra a sinistra di quella delle unità:
175/6; floor(175/6); signif(175/6, 3); round(175/6, 3); round(175/6, -1) #29.16667 29 29.2 29.167 30
Ecco come calcolare la distribuzione percentuale. Metto i singoli dati in una collezione, c(...), che assegno (<-) ad una variabile che chiamo, ad esempio, dati. Faccio il rapporto tra i singoli dati e la loro somma, e lo moltiplico per 100. Nell'esempio seguente arrotondo le percentuali ai decimi.
dati <- c(144291, 18461, 71352, 210285, 119857, 386256) dati/sum(dati) # 0.15180505 0.01942237 0.07506770 0.22123573 0.12609863 0.40637053 round( dati/sum(dati)*100, 1) # 15.2 1.9 7.5 22.1 12.6 40.6
| (1.3) |
| ||||||||||||||||||||||||
Se voglio estrarre un dato da una collezione devo metterne l'indice tra parentesi quadre. Un esempio, riferito ai dati precedenti:
dati[3]; dati[3]/sum(dati)*100 # 71352 7.50677
Ecco come rappresentare la distribuzione precedente con un istogramma a barre. A sinistra la rappresentazione di come sono distribuiti i dati, a destra quella della loro distribuzione percentuale, che abbiamo appena calcolato, con qualche "abbellimento".
Il diagramma a sinistra (relativo ai dati) è stato eseguito battendo semplicemente:
barplot(dati)
Quello a destra (relativo alla loro distribuzione percentuale) con i seguenti comandi, dove:
– in una variabile (classi) ho messo i nomi che voglio dare alle varie colonne poi, con names.arg, nel
comando barplot;
– con abline traccio delle rette, che con h=c(...) specifico essere orizzontali, alle quote
5, …, 40, di colore rosso e (lty=3) tratteggiate.
classi <- c("al","tb","v","ab","tr","a")
barplot(dati/sum(dati)*100,las=1,names.arg=classi)
abline(h = c(5,10,15,20,25,30,35,40), col="red", lty=3)
Ecco come rappresentare gli stessi dati con un areogramma circolare e con un diagramma a striscia, che hanno il vantaggio di facilitare il confronto tra i singoli dati e il totale, ma lo svantaggio di rendere meno facile il confronto tra i singoli dati.
pie(dati,classi) # L'areogramma barplot(height=cbind(dati/sum(dati)*100),horiz=TRUE) # Il diagramma a striscia rug(seq(0,100,5),ticksize=-1/4) # Per aggiungere delle tacche alla striscia
2. Rappresentazioni statistiche di dati numerici non classificati. Valori medi
Consideriamo i tempi, in secondi, tra una telefonata e l'altra in arrivo presso una organizzazione di vendite televisive, raccolti nel file difarr.htm allegato; se copio tale file in R vengono caricati nella variabile difarr i tempi. Ecco i primi dati raccolti nel file:
7.0, 6.0, 1.3, 2.1, 15.8, 5.3, 3.0, 28.0, 2.8, 10.4, 3.0, 2.8, 1.5, 10.3, ...
Per rappresentare graficamente il file mediante un istogramma basta che azioni il comando hist(difarr):
ottengo la rappresentazione sotto a sinistra:
In realtà l'istogramma raffigurato è colorato e ha una griglia che ne segna le ordinate. Ecco i comandi con cui è stato tracciato:
hist(difarr, col="yellow") abline(h = seq(5,50, 5),lty=3)
 |
Le ordinate della griglia sono state introdotte con seq(5,50,5) che indica
la sequenza che parte da 5 e arriva a 50 con passo 5 (con seq(5,50,1) l'avrei ottenuta più fitta, con passo 1). Altrimenti
avrei potuto usare, come nell'esempio precedente, c(5,10,15,20,25,30,35,40,45,50). Il numero della classi è scelto automaticamente. Potrei variarlo; ecco come ottenere l'istogramma a sinistra: hist(difarr,col="yellow",nclass=20). La scala verticale dipende dalla quantità di dati esaminati. Se voglio una scala che non dipenda da essa ma sia tale che l'istogramma abbia area pari al 100%, ossia ad 1, in modo che l'estensione di ogni colonnina rappresenti la corrispondente percentuale di dati, devo usare i comandi seguenti. Ottengo l'istogramma sopra a destra. Le ordinate rappresentate, ossia le altezze delle colonnine, vengono chiamate densità di frequenza. |
dev.new() hist(difarr,freq=FALSE, col="yellow") abline(h = seq(0.01,0.1, 0.01), lty=3) | Abbiamo aggiunto il comando dev.new() che apre una nuova finestra grafica, in modo da mantenere la visualizzazione anche dell'istogramma precedente. |
Vediamo alcuni modi per sintetizzare e mettere meglio in luce alcuni aspetti dei dati precedenti. Partiamo da dati più semplici, usando un modo comodo per rappresentare sequenze di numeri interi: m:n rappresenta i numeri successivi che partono da m ed arrivano ad n.
-5:10; length(-5:10) # -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 # 16 a <- c(-5:10,19,-8); a; sort(a) # -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 19 -8 # -8 -5 -4 -3 -2 -1 0 1 2 3 4 5 6 7 8 9 10 19 n <- 1:9; summary(n) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 1 3 5 5 7 9 n <- 1:8; summary(n) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 1.00 2.75 4.50 4.50 6.25 8.00 min(n); max(n); mean(n); median(n) # 1 8 4.5 4.5 quantile(n,c(0.25,0.5,0.75), name=FALSE) # 2.75 4.50 6.25
length dà la lunghezza di una sequenza di dati, sort li mette in ordine,
summary ne sintetizza alcune informazioni: il minimo e il massimo, la mediana (ossia il valore che sta al centro dei dati, messi in
ordine di grandezza), il 1º e il 3º quartile (ossia i valori che delimitano il primo quarto e il terzo quarto dei dati, messi
in ordine - il 2º quartile, ovviamente, è la mediana), e la media, che è pari a sum(n)/length(n).
Nel caso dei numeri da 1 ad 8, visto che al centro vi sono due dati, 4 e 5, R stampa il valore a metà di
essi, 4.5; si comporta in modo analoghi con i quartili.
Veniamo ai nostri dati, che, come vediamo, sono 135.
In questo caso la mediana è più piccola della media.
Infatti la mediana corrisponde all'ascissa che divide l'istogramma in due parti di area uguale mentre la media risente
della presenza della coda a destra: se allontano un dato a destra ancora di più la mediana, ossia il dato centrale, non cambia mentre aumenta
il valore della media.
E il primo quartile è molto vicino al minimo, mentre il terzo quartile, risentendo della coda a destra, è molto
lontano dal massimo.
length(difarr); summary(difarr) #135 #Min. 1st Qu. Median Mean 3rd Qu. Max. #0.100 2.400 6.100 9.021 12.700 45.100 |
 |
 |
boxplot(difarr,horizontal=TRUE,range=0,col="yellow") points(quantile(difarr,0.05),1,pch=20) points(quantile(difarr,0.95),1,pch=20) |
Il diagramma sopra a sinistra, chiamato boxplot, rappresenta graficamente le informazioni precedenti. Il "box" centrale va dal primo quartile, 2.4, al terzo quartile, 12.7, e la barra che lo divide corrisponde alla mediana, 6.1. I "baffi" partono dal minimo, 0.1, e arrivano al massimo, 45.1. Ci sono, poi, due pallini, tracciati col comando points. L'opzione pch=20 è stata aggiunta per fare pieno il pallino; come ascissa dei pallini sono stati presi i valori che delimitano il primo 5% dei dati e il primo 95% dei dati; l'ordinata dei pallini è 1 (se si fossero rappresentati più boxplot, al successivo sarebbe corrisposta l'ordinata 2).
I dati precedenti erano arrotondati ai decimi di secondo e gli istogrammi realizzati avevano classi
ampie 5 secondi. Quindi possiamo interpretare i dati come se fossero "esatti". In altri situazioni occorre tener conto di come i dati
sono approssimati, se per arrotondamento o troncamento. Ad esempio se le età dei giocatori di una squadra di calcio fossero le
seguenti:
eta <- c(31,28,23,29,25,33,24,21,27,33,20,31,24,20,25,23,20,26,24,28,28,26,27,32)
siccome si tratta di dati troncati (quando una persona dice di avere 23 anni intende dire che potrebbe anche avere 24 anni meno 1 giorno),
nel calcolare la media, arrotondata ai decimi di anno, devo aggiungere 1/2 e concludere che è 26.7:
mean(eta); mean(eta+1/2)
#26.16667 26.66667
Anche l'istogramma dovrei tracciarlo di eta+1/2.
|
3. Il caso dei dati interi e di quelli già classificati
|  | ||
| Questa volta in abline ho messo h=axTicks(2) senza specificare le ordinate: axTicks(2) sceglie automaticamente un certo numero di valori nell'intervallo dell'asse y rappresentato. | |||
Se dispongo di dati già classificati posso procedere, per semplicità, usando un file per R già predisposto. Vediamo come, riferendoci ai dati relativi alle età dei morti di sesso femminile nel 2006 classificate negli intervalli di anni [0,5), [5,10), …, [110,115):
intervalli <- seq(0,115,5) #0, 5, 10, ..., 115: gli estremi degli intervalli
morti <- c(39,4,5,9,10,12,16,25,40,64,102,155,236,355,593,1004,1666,2223,2077,1040,297,27,1)
# i morti nelle varie classi di età (sono un dato in meno degli estremi degli intervalli)
source("http://macosa.dima.unige.it/R/daticlas.txt") # carico il file
aiuto # e un apposito help
interv <- intervalli; freq <- morti # metto i dati nelle variabili richieste dall'help
source("http://macosa.dima.unige.it/R/daticlas.txt")
# ottengo le uscite numeriche e grafiche seguenti:
# Min. 1st Qu. Median Mean 3rd Qu. Max.
# 0.00 79.15 86.49 83.98 92.27 115.00
In questo caso la media è più bassa della mediana in quanto vi è una coda a sinistra. Se i dati fossero stati troncati avrei dovuto procedere aggiungendo mezza unità dell'ordine di grandezza dell'ultima cifra, come visto alla fine del paragrafo precedente. |  |
4. Il calcolo delle probabilità
Consideriamo i dati sulle età delle femmine morte in Italia nel 2006 riportati alla fine del paragrafo precedente. Qual è la probabilità che, estratto del tutto a caso il nome di una di esse, la sua età di morte sia maggiore o eguale a 90 anni? Basta che calcoli la relativa percentuale; ottengo 34.42%:
sum(c(2077,1040,297,27,1))/sum(morti)*100 # 34.42
Se indico con Età l'età di morte di una generica donna morta nel 2006, posso esprimere con la formula
Consideriamo un'altra situazione. Sta per disputarsi la partita Roma-Lazio. Gigi ritiene che la Roma 25 su 100 vincerà e 40 su 100 pareggerà.
Qual è la probabilità per Gigi che vinca la Lazio? La situazione,
indicato con R il risultato della partita ("1", "2" o "X"), può essere sintetizzata
così:
poiché Pr(R = "1") + Pr(R = "2") + Pr(R = "X") = 100%,
Pr(R = "2") = 100% − Pr(R = "1") − Pr(R = "X") = 100% − 25% − 40% = 35%.
In entrambi gli esempi ho associato ad alcuni eventi A un numero compreso tra 0 e 1 (=100%) come Pr(A) (probabilità di A). Ho poi dedotto le probabilità relative ad altri eventi applicando a Pr alcune delle proprietà che si erano già usate per le frequenze percentuali.
 |
• Pr(not A) = 100% – Pr(A) • Pr(A or not A) = 100% =1 • Pr(A and not A) = 0 |
 |
• Pr(A1 or A2 or A3 or …) =
Pr(A1) + Pr(A2) + Pr(A3) + … se A1, A2, A3, … sono tra loro incompatibili, cioè se due qualunque eventi Ai e Aj non possono essere veri contemporaneamente (proprietà additiva) |
 |
Di fronte a valutazioni del tipo Pr(A OR B) con A e B eventi non incompatibili, si usa la proprietà: • Pr(A or B) = Pr(A) + Pr(B) – Pr(A and B) |
Naturalmente, a seconda di come si scelgono le valutazioni iniziali, per la stessa situazione si possono ottenere diverse misure di probabilità. Le valutazioni iniziali possono essere dedotte dall'esperienza o da considerazioni di tipo fisico o da propri convincimenti o …. Devono comunque essere tali da non condurre a contraddizioni: a partire da esse, applicando ripetutamente le proprietà sopra elencate, non posso ottenere valutazioni diverse per uno stesso evento, non posso ottenere probabilità negative o superiori al 100%, … (ad es. non posso valutare 60% la probabilità che nella prossima partita Roma-Lazio vinca la Roma e 50% che pareggino; verrebbe contraddetta la prima proprietà). Si osservi che il ruolo delle valutazioni iniziali mostra come anche in questo caso, come in altri discussi in altre voci, le conoscenze matematiche non sono di per sé sufficienti per modellizzare o risolvere "razionalmente" un problema.
Facciamo un esempio in cui è facile fare valutazioni probabilistiche. Il lancio di un dado equo, ossia un dado che,
diversamente da quello costruito col cartoncino considerato nel paragrafo precedente, abbia tutte le facce "equiprobabili", ossia, indicata con U l'uscita, tale che
Pr(U=1) = Pr(U=2) = Pr(U=3) = Pr(U=4) = Pr(U=5) = Pr(U=6). Queste sono tutte le 6 possibili uscite. Sia P la probabilità di ciascuna di esse; per la
proprietà additiva P+P+P+P+P+P = 1, ossia P = 1/6.
Se lancio due dadi equi, qual è la probabilità di avere un'uscita pari?
Posso procedere con un grafo ad albero, avendo indicato con U1 ed U2 le possibili uscite dei due dadi:
|
− rappresento con successive diramazioni i diversi esiti possibili per U1 e per U2 (eventualmente raggruppando gli esiti
"sfavorevoli" in un'unica diramazione); − associo agli archi che corrispondono a esiti "favorevoli" la relativa probabilità; − calcolo, per ogni percorso (dal nodo iniziale a un nodo finale) costituito solo da archi "favorevoli", il prodotto delle probabilità associate ai vari archi e lo scrivo a fianco del nodo finale; − sommo i valori così calcolati. I percorsi favorevoli sono U1 = 1, U2 = 3; U1 = 2, U2 = 2; U1 = 3, U2 = 1. Nella figura a fianco sono indicati A, B e C, e corrispondono ciascuno alla probabilità Complessivamente, 1/36 + 1/36 + 1/36 = 1/12. |  |
Potrei simulare il fenomeno anche con il generatore di numeri casuali, ossia usando il termine runif(n) che assume n valori reali che cadono in [0,1) ciascuno indipendentemente da dove cade il precedente e con distribuzione uniforme: presi comunque due sottointervalli di [0,1) ugualmente ampi, la probabilità che runif(1) cada nell'uno o nell'altro sono eguali. Vediamo come.
# L'istogramma a sinistra, con distribuzione uniforme hist(runif(1e6),freq=FALSE,col="yellow") # Quello a destra, relativo al lancio di 2 dadi equi: # [ floor(runif(1)*6)+1) è un numero intero a caso tra 1 e 6 ] n <- 10000; U1 <- floor(runif(n)*6)+1; U2 <- floor(runif(n)*6)+1 hist(U1+U2, seq(1.5, 12.5, 1), probability=TRUE, col="yellow") summary(U1+U2) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 2.000 5.000 7.000 6.973 9.000 12.000
5. Leggi di distribuzione
Considero ancora due leggi di distribuzione, una continua ed una discreta. Considero ancora il caso delle telefonate all'organizzazione di vendite televisive considerato nel pargrafo 2. Esaminiamo, ora, la durata delle telefonate, raccolte nel file durtel.htm allegato; copio tale file in R caricando i dati nella variabile durtel.
hist(durtel,freq=FALSE,col="yellow") summary(durtel) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 3.60 37.20 49.70 50.07 62.73 107.00
Ottengo l'istogramma sotto a sinistra, che ha una forma quasi a campana, con i tre quartili centrali pari a circa 37, 50, 63, ossia col primo e terzo quartile quasi equidistanti dalla mediana, e con la mediana quasi coincidente con la media. Vedremo poi che tale istogramma è approssimabile con una opportuna curva, di cui con dei pallini abbiamo segnato i punti di flesso.
Il secondo caso è discreto. Rappresenta come è distribuito il numero N dei lanci di una moneta equa da effettuare fino a che esce testa.
Lo studiamo teoricamente: al 50% testa esce al primo lancio (N=1), nei casi rimanenti al 50% (ossia, complessivamente, al 25%) esce al secondo lancio (N=2),
al 12.5% esce per N=3, al 6.25% esce per N=4, ….
[il grafico è stato tracciato per gli input 1, 2,..., 9: ho messo n=9 in modo siano scelti 9 punti;
type="h" - da histogram - serve per tracciare solo le altezze,
lwd=3 per fare le linee più spesse, xlim=c(0,9) per scegliere l'intervallo da rappresentare]
f <- function(x) 1/2^x plot(f,1,9,n=9,col="brown",type="h", lwd=3, xlim=c(0,9)) abline(h=axTicks(2), v=axTicks(1), lty=3); abline(h=0, v=0, col="blue")
La somma delle altezze di tutte le colonnine deve essere 100%, ossia 1. E infatti: 1/2 + 1/4 + 1/8 + … fa 1. Non deve stupire che la somma di infiniti numeri sia un numero: sin dai primi anni scolastici avete imparato che, ad esempio, 0.3 + 0.03 + 0.003 + … = 0.333… = 1/3.
Qual è il numero medio di lanci da effettuare affinché esca testa? Devo fare la media pesata. Facciamo un esempio. Se il valore 3 esce 5 volte e il 9 esce 10 volte, la media è (3·5+9·10)/(5+10) = 7. Ma posso anche fare 3·5/(5+10) + 9·10/(5+10), ossia fare la somma dei valori moltiplicati per le frequenze relative. Nel nostro caso:
# una prova per 4 valori 1*1/2 + 2*1/2^2 + 3*1/2^3 + 4*1/2^4 # 1.625 # generalizziamo: n <- 1:4; sum(n*1/2^n); n <- 1:10; sum(n*1/2^n); n <- 1:50; sum(n*1/2^n) # 1.625 1.988281 2
La media, 1/2 + 2/2^2 + 3/2^3 + 4/2^4 + 5/2^5 + …, è 2.
A questo punto dobbiamo mettere a punto degli strumenti per verificare quando una curva rappresenta
una funzione di densità, ossia una funzione sul cui grafico tende a stabilizzarsi l'istogramma della
frequenza relativa di una variabile casuale continua all'aumentare del numero delle prove.
Nel caso di un fenomeno discreto occorre che la somma delle frequenze relative sia 1, nel caso continuo occorre che l'area
tra tale grafico e l'asse delle x sia 1. Lo strumento per rappresentare quest'area è l'integrale.
Nel caso del generatore di numeri casuali considerato alla fine del paragrafo precedente la cosa è facile: l'area di
un quadrato di base 1 e altezza 1 è 1. Vediamo i casi che corrispondono ai tempi di arrivo e alle durate delle telefonate, considerati nel paragrafo 2 e all'inizio
di questo paragrafo.
La distribuzione esponenziale negativa
|
A destra è riprodotto l'istogramma di distribuzione dei tempi tra una telefonata e l'altra e, sovrapposto, il grafico di una funzione che sembra approssimarlo. Il tempo medio di attesa era di circa 9 secondi. Fenomeni di questo tipo (come ad es. anche la distanza temporale tra la venuta al semaforo di un'auto e quella della successiva, nel caso di un semaforo preceduto da un lungo tratto di strada senza impedimenti) hanno una distribuzione, chiamata esponenziale, che ha come funzione di densità la seguente, dove a è il reciproco della media (nel nostro caso a è circa 1/9): f : x → a·e− a x (x > 0) m <- mean(difarr); a <- 1/m; f <- function(x) a*exp(-a*x) hist(difarr,freq=FALSE, col="yellow", ylim=c(0,a)) abline(h=axTicks(2),v=axTicks(1),lty=3); abline(h=0,v=0,col="blue") curve(f, add=TRUE,lwd=2); points(m,0,pch=19) |
 |
Verifichiamo che l'area sottesa è 1. La funzione esponenziale ha la caratteristica di avere Dx exp(x) = exp(x), e, quindi,
Dunque: ∫ [0,∞) f =
Abbiamo già osservato che la media è il reciproco di a (m = 1/a). Verifichiamolo precisando il significato di "media" nel caso continuo.
Nel caso discreto essa è la somma dei valori moltiplicati per le frequenze relative, ovvero
moltiplicati per le probabilità.
Nel caso continuo diventa:
∫ [0,∞) x·f(x) dx =
∫ [0,∞) x·a·exp(-a·x) dx = 1/a
[qui vediamo solo come ottenere questo valore con WolframAlpha:
integrate x*a*exp(-a*x) dx from 0 to oo
o con R:
a <- 1/9; f <- function(x) x*a*exp(-a*x); integrate(f,1,Inf) - ottengo 9]
La distribuzione gaussiana (o normale) Torniamo alla durata delle telefonate, di cui all'inizio del paragrafo abbiamo visto l'istogramma. La funzione col cui
grafico è approssimabile è una particolare funzione di densità f,
detta normale o gaussiana, così definita,
dove:
Chiamata varianza la media dei "quadrati" degli scarti dal valor medio, si chiama scarto quadratico medio (e si indica spesso con la lettera greca σ) la sua radice quadrata:
Si può dimostrare che, indipendentemente dai valori che posso avere nei vari casi, il grafico è simmetrico rispetto a y=m, σ è la distanza dei punti di flesso da m e l'integrale di f tra m-kσ m+kσ dipende solo dal valore di k (a destra ne sono riportati alcuni valori approssimati). |  |
Vediamo come è stato fatto il grafico della gaussiana sovrapposto all'istogramma riportato all'inizio del paragrafo. In R la gaussiana (o distribuzione normale) può essere introdotta usando semplicemente il nome dnorm.
hist(durtel,freq=FALSE,col="yellow") # La definizione della varianza e dello s.q.m. V <- function(dati) sum((dati-mean(dati))^2)/length(dati) sqm <- function(dati) sqrt(V(dati)) M <- mean(durtel); S <- sqm(durtel) # la gaussiana nel nostro caso (sd sta per "deviazione standard" e indica lo sqm) f <- function(x) dnorm(x, mean=M, sd=S ); curve(f, add=TRUE) # l'integrale tra m e infinito (deve essere la metà di 1): integrate(f, M, Inf) #0.5 with absolute error < 7.5e-06 # l'integrale tra m-σ e m+σ (nel grafico arrotondato a 68.3%) integrate(f, M-S, M+S ) #0.6826895 with absolute error < 7.6e-15 # Il 1º e il 3º quartile corrispondono a m-kσ e m+kσ con k=0.67448975 integrate(f, M-0.67448975*S, M+0.67448975*S) #0.5 with absolute error < 5.6e-15 # i "pallini": points( M+S,f(M+S), pch=19 ); points( M-S,f(M-S), pch=19 )
6. Il teorema limite centrale. Altre leggi di distribuzione
Consideriamo un ulteriore esempio di legge di distribuzione. Supponiamo che una fabbrica di biscotti disponga di un forno che bruciacchi i biscotti con la frequenza p (ossia questa è la probabilità che un biscotto sia bruciacchiato) e che venda i biscotti in confezioni da n pezzi. Qual è la probabilità che in una confezione il numero N dei biscotti bruciacchiati sia k?
Se n = 6 e p = 1/8 la probabilità che esattamente i primi 4 biscotti siano bruciacchiati è (1/8)·(1/8)·(1/8)·(1/8)·(7/8)·(7/8) = (1/8)4·(7/8)2. Questo valore dobbiamo moltiplicarlo per i possibili sottoinsiemi di 4 elementi che possono essere formati da un insieme di 6 elementi. Questo numero viene in genere indicato C(6,4) e chiamato numero delle combinazioni di 6 elementi 4 a 4, ed è pari al numero dei quartetti ordinati (6·5·4·3: 6 modi di prendere il primo elemento, 5 di prenderne il secondo; …) diviso per i modi in cui posso ordinare 4 elementi (4·3·2·1).
C(6,4) = (6·5·4·3)/(4·3·2·1) = 6/4·5/3·4/2·3/1 = 6·5/2/1 = 3·5 = 15
Nel software in genere C(n,k) è indicato choose(n,k). Ecco i sottoinsiemi di 0, 1, …, 6 elementi di un insieme di 6 elementi:
k <- 0:6; choose(6, k) #1 6 15 20 15 6 1
Dunque, nel nostro caso particolare, la probabilità che vi siano 4 biscotti bruciacchiati è C(6,4)·(1/8)4·(7/8)2 = 15·(1/8)4·(7/8)2 = 0.0028038 = 0.28% (arrotondando).
In generale:
Pr(N = k) = C(n, k) · pk · (1 – p)n–k
Questa legge di distribuzione viene chiamata legge di distribuzione binomiale (o di Bernoulli). Si applica a tutte le situazioni in cui si ripete n volte la prova su una variabile casuale che può assumere solo due valori, in cui p è la probabilità di uno di questi due valori e N è il numero delle volte in cui questo valore esce.
Ecco qualche grafico. Il primo, tracciato con plot(0:5, dbinom(0:5, 5,0.2), type="h",lwd=5), rappresenta la binomiale con p = 0.2 e n = 5. Gli altri, tracciati analogamente, rappresentano quelle per p = 0.2 e n = 10, n = 30.
Si noti come all'aumentare di n l'istogramma di distribuzione tende ad assumere una forma simile alla
gaussiana. Come mai? La binomiale di ordine n è ottenibile come somma di n termini uguali ad una variabile casuale ad uscite in 0 ed 1
(ad esempio, nel caso dei 6 biscotti, è la somma di 6 variabili ad uscite in 0 od 1). Consideriamo un'altra variabile casuale che rappresenta
la ripetizione di esperimenti, ad esempio la somma di n termini pari a 9·√RND+RND², dove con RND abbiamo indicato
le uscite del generatore di numeri casuali. Ecco i grafici, per n pari ad 1, 2 e 20, ottenuti con:
x <- NULL; for(i in 1:2000) x[i] <- sum(9*sqrt(runif(n))+runif(n)^2); hist(x,freq=FALSE)
Si può dimostrare che se Ui (i intero positivo) sono variabili casuali (numeriche)
indipendenti con la stessa legge di distribuzione, allora
al crescere di n la variabile casuale
Tale proprietà, nota come teorema limite centrale, oltre ad essere utile per approssimare la binomiale nel caso in cui n sia molto grande, è fondamentale nelle applicazioni. Vediamo un esempio.
Voglio determinare il valor medio M(P) dove P è il "peso di un abitante adulto maschio" (di un certo paese). Indico con σ lo sqm di P. Rilevo i pesi P1, P2, …, Pn di n persone.
Dividendo per n ho
Quanto qui detto per P vale per ogni variabile casuale.
Il valore di σ devo già conoscerlo in base a considerazioni di qualche tipo oppure posso
man mano approssimarlo con la radice quadrata
della varianza sperimentale:
si può dimostrare che, fissato n, la varianza di
Ovvero come σ devo prendere
il valore sperimentale moltiplicato per
√(n/(n−1)). Ovvero devo prendere il secondo dei valori che, in modo non
molto corretto ma ormai diffuso, vengono in genere indicati nel modo seguente (dove
xi e μ sono dati e media):
|
|
Questi due termini sono spesso chiamati rispettivamente deviazione standard teorica e deviazione standard corretta o non distorta o statistica. Spesso sono entrambi chiamati semplicemente deviazione standard: sta al lettore capire quale uso si sta facendo. Il software R, ad es., indica con sd la deviazione standard corretta. Comunque quando n è abbastanza grande i due numeri hanno una piccola differenza relativa.
Un'altra legge di distribuzione che ha andamento abbastanza simile a quello della binomiale e che trova applicazione soprattutto in fisica e in biologia, in situazioni in cui gli eventi accadono abbastanza "raramente", è la legge di Poisson. Rimandiamo agli Oggetti Matematici chi voglia prenderne conoscenza.
7. Dipendenza e indipendenza stocastica
(a) Qual è la probabilità che alzando 2 volte un mazzo (nuovo) di carte da scopa ottenga sempre
una carta di denari?
(b) Qual è la probabilità che estraendo 2 carte dal mazzo queste siano entrambe di denari?
|
• Nel caso della alzata, avendo supposto il mazzo nuovo (e non truccato e mescolato bene) posso ritenere che, tagliandolo, le carte, e quindi (essendoci 10 carte per ogni seme) anche i semi, a valori in {♥, ♦, ♣, ♠}, escano con distribuzione uniforme: l'uscita di una carta di denari ha la stessa probabilità di quella di una di fiori o … Posso rappresentare queste due alzate col grafo ad albero a fianco, a due diramazioni. Ho 1/4 di probabilità di estrarre denari alla prima alzata ed 1/4 di estrarlo alla seconda. La probabilità cercata è dunque 1/4·1/4 = 1/16. |
|
|
• Anche nel caso della estrazione posso ritenere equiprobabili le carte, e i semi, del mazzo. Ma mentre alla prima estrazione ho 1/4 di probabilità di estrarre una carta di denari, alla seconda estrazione la probabilità cambia. Le carte da cui effettuare l'estrazione sono, ora, una in meno, e, se ho estratto una carta di denari alla prima estrazione, le carte di denari rimaste sono 9. Il grafo a destra illustra la situazione. La probabilità cercata in questo caso è 1/4·9/39 = 3/4/13 = 0.0576923… = 5.8%. |
 |
Indichiamo con le variabili casuali S1 e S2 il seme della prima uscita e quello della seconda. Nel caso della alzata S1 e S2 sono indipendenti: qualunque seme abbia la 1ª carta, la probabilità che la 2ª abbia un certo seme è sempre la stessa. Ciò corrisponde al fatto che il grafo relativo all'alzata si riproduce allo stesso modo passando da una diramazione alla successiva. Per calcolare Pr(S1=♦ and S2=♦) posso fare direttamente Pr(S1=♦)·Pr(S2=♦) = 1/4·1/4 = 1/16.
Nel caso della estrazione S1 e S2 non sono indipendenti: ad es. Pr(S2=♦) (la probabilità che la 2ª carta sia di ♦) dipende dal valore assunto da S1 (cioè dal seme della 1ª carta). Ciò corrisponde al fatto che il grafo relativo alla estrazione non si riproduce allo stesso modo passando da una diramazione alla successiva: al primo arco "♦" è associata la probabilità 1/4, al secondo arco "♦" è associata la probabilità 9/39.
Due variabili casuali X e Y sono probabilisticamente indipendenti se sono indipendenti gli eventi A e B comunque prenda A evento relativo a X (condizione in cui compare solo la variabile X) e B evento relativo a Y (condizione in cui compare solo variabile Y): conoscere qualcosa su come si manifesta X non modifica le mie aspettative sui modi in cui può manifestarsi Y, e viceversa. Altrimenti sono probabilisticamente dipendenti. Esempio:
− sapere qualcosa a proposito del seme della 1ª carta estratta cambia le mie valutazioni sul seme che potrebbe avere la 2ª carta estratta: il seme della 1ª estrazione e quello della 2ª sono variabili casuali dipendenti.
Ricordiamo che il concetto di dipendenza ora introdotto è diverso da quello impiegato per esprimere il legame tra due grandezze quando una varia in funzione dell'altra. L'avverbio "probabilisticamente" (o l'equivalente avverbio "stocasticamente") evidenzia questa differenza. Se non ci sono ambiguità, questo avverbio viene omesso.
|
8. Sistemi di variabili casuali Consideriamo tre diverse situazioni in cui abbiamo una coppia U = (X,Y) di variabili casuali e la rappresentazione grafica di come esse si distribuiscono. Nel caso di una uscita avevamo delle curve; in questo caso abbiamo delle superfici. Nel primo caso avevamo che l'area tra curva e asse x valeva 1; ora abbiamo che il volume tra superficie e piano xy vale 1. Il grafico (A) è riferito alla caduta di proiettili in un bersaglio circolare (il cerchio centrato nell'origine e di raggio 1), nell'ipotesi che la distribuzione sia uniforme, ossia che i proiettili arrivino senza privilegiare alcuna parte del bersaglio: vedi la distribuzione di un po' di uscite nella figura (A'). X ed Y sono le coordinate dei punti in cui cadono i proiettili. In parti del cerchio di eguale superficie i proettili cadono con eguale probabilità; a ciò corrisponde il fatto che la superficie rappresentata ha altezza costante. Il solido che sta tra il cerchio e il piano xy ha volume 1 (la sua altezza h vale 1/π); lo spicchio con le x e le y positive ha volume 1/4. Il grafico (B) rappresenta la distribuzione di (X,Y) con X e Y altezze di un uomo e una donna sorteggiati a caso. I valori sono stati traslati in modo che le altezze medie valgano 0. Il grafico (C) rappresenta in modo analogo la distribuzione di (X,Y) con X e Y altezze di marito e moglie di una coppia sorteggiata a caso (in (C') sono rappresentate un po' di coppie): l'altezza di uomini sposati con donne di una certa altezza ha andamento più o meno gaussiano, ma la loro altezza media è maggiore di quella degli uomini sposati con donne più basse (uomini più alti tendenzialmente sposano donne più alte: non è affatto vero che l'amore è cieco!). |
 |
Nel caso (A) i valori che può assumere una delle due variabili
è condizionato da quello che assume l'altra: se X è vicino ad 1 Y per forza deve essere vicino a 0.
Nel caso (B) X e Y sono indipendenti: comunque sezioni la superficie con piani paralleli ai piani xz e yz
ottengo grafici con andamenti simili: hanno massimo e punto di flesso collocati nella stessa posizione.
Nel caso (C), come abbiamo già osservato, X ed Y sono dipendenti, ma la dipendenza
è in un qualche senso "più forte" di quella del caso (A): al crescere di
X anche Y tende a crescere, ossia X ed Y sono "correlate".
Vediamo come si può quantificare questa idea di correlazione. Si introduce il:
| coefficiente di correlazione: r X,Y = |
|
Si può dimostrare che se X
e Y sono dipendenti deterministicamente e legate da una relazione lineare Y = aX + b il coefficiente di
correlazione assuma il valore assoluto massimo. Vale 1 se l'andamento è crescente e −1
se è decrescente. Quindi, in generale,
In R il coefficiente di correlazione tra X ed Y è indicato cor(X,Y).
Ecco, ad esempio, la parte iniziale del file battito.txt:
# Indagine sugli studenti di un corso universitario, dal manuale di MiniTab # battiti prima di eventuale corsa di 1 min # battiti dopo # fatta corsa (1 si`;0 no; a seconda di esito di lancio moneta) # fumatore (1 si`;0 no) # sesso (1 M; 2 F) # altezza # peso # attivita` fisica (0 nulla;1 poca;2 media; 3 molta) Num,BatPrima,BatDopo,Corsa,Fumo,Sesso,Alt,Peso,Fis 01,64,88,1,0,1,168,64,2 02,58,70,1,0,1,183,66,2 ... |
che potrei leggere col comando:
readLines("http://macosa.dima.unige.it/R/battito.txt",n=11)
Se raccolti su una usuale tabella i dati assumerebbero questo aspetto:
| bat | bat.dopo | corsa | fumo | sesso | alt | peso | attiv |
| 64 | 88 | 1 | 0 | 1 | 168 | 64 | 2 |
| 58 | 70 | 1 | 0 | 1 | 183 | 66 | 2 |
| ... | ... | ... | ... | ... | ... | ... | ... |
Vedo che si sono righe di commento inizianti con # e 1 riga di "intestazioni", e che
dati sono separati da ",". Il software R salta automaticamente le righe inizianti
con #. Metto, dunque, i dati in una "table" con:
dati <- read.table("http://macosa.dima.unige.it/R/battito.txt",sep=",",header=TRUE)
e ne visualizzo la struttura con:
str(dati)
con cui ottengo
'data.frame': 92 obs. of 9 variables: $ Num : int 1 2 3 4 5 6 7 8 9 10 ... $ BatPrima: int 64 58 62 66 64 74 84 68 62 76 ... $ BatDopo : int 88 70 76 78 80 84 84 72 75 118 ... $ Corsa : int 1 1 1 1 1 1 1 1 1 1 ... $ Fumo : int 0 0 1 1 0 0 0 0 0 0 ... $ Sesso : int 1 1 1 1 1 1 1 1 1 1 ... $ Alt : int 168 183 186 184 176 184 184 188 184 181 ... $ Peso : int 64 66 73 86 70 75 68 86 88 63 ... $ Fis : int 2 2 3 1 2 1 3 2 2 2 ...
Col comando summary posso avere una informazione sintetica di tutte le variabili; ma la prima colonna contiene solo il numero d'ordine, quindi invece di usare summary(dati) esamino solo le colonne dalla 2 alla 9 con:
summary(dati[2:9])
BatPrima BatDopo Corsa Fumo
Min. : 48.00 Min. : 50 Min. :0.0000 Min. :0.0000
1st Qu.: 64.00 1st Qu.: 68 1st Qu.:0.0000 1st Qu.:0.0000
Median : 71.00 Median : 76 Median :0.0000 Median :0.0000
Mean : 72.87 Mean : 80 Mean :0.3804 Mean :0.3043
3rd Qu.: 80.00 3rd Qu.: 85 3rd Qu.:1.0000 3rd Qu.:1.0000
Max. :100.00 Max. :140 Max. :1.0000 Max. :1.0000
Sesso Alt Peso Fis Min. :1.00 Min. :154.0 Min. :43.00 Min. :0.000 1st Qu.:1.00 1st Qu.:167.8 1st Qu.:57.00 1st Qu.:2.000 Median :1.00 Median :175.0 Median :66.00 Median :2.000 Mean :1.38 Mean :174.4 Mean :65.84 Mean :2.109 3rd Qu.:2.00 3rd Qu.:183.0 3rd Qu.:70.25 3rd Qu.:2.000 Max. :2.00 Max. :190.0 Max. :97.00 Max. :3.000Come analizzare un singolo campo (le altezze) e suoi sottocampi (quelle femminili e quelle maschili):
Ecco la matrice di correlazione, che sintetizza le correlazioni tra tutte le diverse variabili di "battito" [avrei potuto battere solo cor(dati[2:9]) ottenendo valori di 8 cifre; con print(…,2) posso "accorciare" i valori; avrei potuto usare, con esiti diversi, round(…,3)]:
print(cor(dati[2:9]),2)
BatPrima BatDopo Corsa Fumo Sesso Alt Peso Fis
BatPrima 1.000 0.616 0.0523 0.129 0.29 -0.211 -0.203 -0.0626
BatDopo 0.616 1.000 0.5768 0.046 0.31 -0.153 -0.166 -0.1411
Corsa 0.052 0.577 1.0000 0.066 -0.11 0.224 0.224 0.0073
Fumo 0.129 0.046 0.0656 1.000 -0.13 0.043 0.201 -0.1202
Sesso 0.285 0.309 -0.1068 -0.129 1.00 -0.709 -0.710 -0.1050
Alt -0.211 -0.153 0.2236 0.043 -0.71 1.000 0.783 0.0893
Peso -0.203 -0.166 0.2240 0.201 -0.71 0.783 1.000 -0.0040
Fis -0.063 -0.141 0.0073 -0.120 -0.10 0.089 -0.004 1.0000
Tra altezza e peso vi è un alto coefficiente di correlazione: 0.78.
Se mi restringo a una sottopopolazione più omogenea (quella femminile o quella maschile,
che hanno pesi e altezze con medie abbastanza diverse) potrei aspettarmi un
coefficiente maggiore. Ma se estraggo la popolazione femminile ottengo 0.52. Perché?
cor(dati[7],dati[8]) # Peso # Alt 0.7826331 cor(subset(dati[7],dati$Sesso==1),subset(dati[8],dati$Sesso==1)) # Peso # Alt 0.5904648 cor(subset(dati[7],dati$Sesso==2),subset(dati[8],dati$Sesso==2)) # Peso # Alt 0.5191614Se traccio in rettangoli cartesiani uguali, i grafici di dispersione
par(mfrow=c(1,3), mar=c(3,3,2,1))
plot(c(150,200),c(40,100),type="n",xlab="", ylab=""); griglia("blue")
points(dati$Alt,dati$Peso)
plot(c(150,200),c(40,100),type="n",xlab="", ylab=""); griglia("blue")
points(subset(dati$Alt,dati$Sesso==1),subset(dati$Peso,dati$Sesso==1),col="blue")
plot(c(150,200),c(40,100),type="n",xlab="", ylab=""); griglia("blue")
points(subset(dati$Alt,dati$Sesso==2),subset(dati$Peso,dati$Sesso==2),col="red")
Questo esempio mette in luce come le statistiche che si ottengono sono spesso ingannevoli. In casi come questo, abbastanza frequenti, il problema è dovuto alla presenza di due sottopopolazioni con caratteristiche differenti.
Poi occorre tener conto che quelle individuate sono solo relazioni statistiche, non di
causa-effetto. Ad esempio nel caso della correlazione
tra le colonne "battito dopo" e "corsa" di "battito"
c'è effettivamente una relazione causale (l'aver fatto la corsa influenza il battito cardiaco).
Ma quando nel caso di uno studio statistico sulle condizioni delle famiglie è emersa una forte correlazione
negativa fra il loro consumo di patate e la superficie dell'abitazione in cui vivono, essa non è
da interpretare come conseguenza di una relazione di causa-effetto: è semplicemente dovuta al fatto che le
famiglie benestanti abitano in genere in case di maggiori dimensioni e, nello stesso tempo, consumano meno patate
delle altre famiglie privilegiando cibi più costosi, come la carne e il pesce.
Purtroppo, specie nei campi medico e socio-psicologico, spesso si fanno collegamenti di questo genere.
Osserviamo, infine, che il coefficiente di correlazione è rilevante se i dati sono molti;
basti pensare che avere tre punti più o meno allineati ha senz'altro un significato diverso dall'averne molti.
Di fronte a dati sperimentali relativi a un sistema (X,Y) per cui si ritiene che Y vari in funzione di X, si può cercare di trovare una funzione F tale che il suo grafico approssimi i punti sperimentali. Vediamo come procedere nel caso in cui X ed Y siano casuali. Si cerca di individuare il tipo di funzione (lineare, polinomiale, esponenziale, …) che si vuole utilizzare. Se si ipotizza che ci sia una relazione lineare che esprima Y in funzione di X, e non si hanno altre informazioni, la tecnica in genere usata è quella dei minimi quadrati, che consiste nel trovare la retta, generica o passante per un punto fissato, a seconda dei casi, che rende minima la somma dei quadrati degli scarti tra i valori sperimentali di Y e quelli che sarebbero stati associati ai valori di X dalla equazione della retta.
Il caso illustrato a fianco è relativo alla ricerca della retta passante per (0,0) X <- c(10,25,35); Y <- c(5,8,11)
plot(c(0,40),c(0,13),type="n"); points(X,Y,pch=19)
r1 <- lm(Y ~ X-1); r1$coefficients # il "-1" specifica che si cerca
# una retta vincolata a passare (0,0)
# X
#0.325641
abline(r1,col="blue",lwd=2)
Nel caso in cui avessi cercato una retta generica, senza imporre il passaggio per (0,0), avrei battuto: r2 <- lm(Y ~ X); r2$coefficients #(Intercept) X # 2.4736842 0.2368421 abline(r2,col="brown",lwd=2) |   |
9. Esercizi Vai qui.