#
# - - - INCONTRO N. 2 - Statistica e probabilità (prima parte) - - -
#
# Vediamo prime statistiche elementari, su dati riferiti a CLASSI
# "QUALITATIVE". Supponiamo di sapere che gli studenti iscritti ad un
# corso di laurea sono di 4 differenti nazioni, 12 della nazione A, 72
# della B, 36 della C e 40 della D. Che cosa possiamo dire di questi
# dati?
# Possiamo mettere questi dati in una variabile, come collezione:
stud <- c(12, 72, 36, 40)
# Possiamo visualizzare le diverse provenienze con un diagramma a barre
# (un istogramma in cui le basi non sono intervalli numerici)
barplot(stud)
# Se vogliamo mettere una griglia (in questo caso ci interessa solo una
# griglia orizzontale) possiamo usare il comando che ci consente di
# tracciare rette (abline), e tracciare quelle di quota 10, 20,..., 70,
# in modo tratteggiato (per "seq" vedi)
abline(h=seq(10,70,10), lty=3)
# Se voglio dare dei nomi alle colonne posso usare il comando "names"
names(stud) <- c("A","B","C","D");barplot(stud);abline(h=seq(10,70,10),lty=3)
# Se voglio posso tracciare colonne orizzontali:
barplot(stud,horiz=TRUE);abline(v=seq(10,70,10),lty=3)
# e mettere orizzonali anche le scritte (con "las=1"):
# (las: labels axis style; altri valori: 2 e 3)
barplot(stud,horiz=TRUE,las=1);abline(v=seq(10,70,10),lty=3)
#
# Gli stessi dati posso rappresentarli con un diagramma a torta
pie(stud)
#
# I dati potrei averli anche in rete. Vediamo cinque possibili formati
# (vedremo più avanti il formato "foglio di calcolo"):
# 1: 'commento: dati A, B, C, D
# 12, 72, 36, 40
# 2: 'commento: dati A, B, C, D
# 12; 72; 36; 40
# 3: 'commento: dati A, B, C, D
# 12
# 72
# 36
# 40
# 4: 'commento: dati A, B, C, D
# A, B, C, D
# 12; 72; 36; 40
# 5: # commento: dati A, B, C, D
# dati <- c(12,72,36,40)
# Si tratta dei file abcd1.txt abcd2.txt abcd3.txt abcd4.txt abcd5.txt
# collocati nella stessa cartella del presente file. Vediamo come
# analizzarli. I modi sono vari; consideriamone uno.
#
# I file sono collocati in rete in:
# http://macosa.dima.unige.it/om/prg/R/agg
# [se i file fossero nel computer potrei cliccare su uno di essi,
# selezionare "Proprietà" e, dalla finestra che si apre, "Generale"
# e, quindi, copiare il "Percorso": "C:/ ... /agg/" oppure potrei
# (=> vedi) scegliere questa come cartella di riferimento; se si sta
# lavorando non in rete conviene procedere così (ossia con Change Dir
# scegliere .../R/agg come cartella di riferimento e al posto dei
# comandi in cui compare http:// ... considerarne una copia da cui
# togliere http://macosa.dima.unige.it/om/prg/R/agg/ - ad es. al posto
# del comando seguente considerare: abcd1 <- readLines("abcd1.txt") )]
# Posso vedere 3 righe del file con:
abcd1 <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/abcd1.txt")
abcd1[1]; abcd1[2]; abcd1[3]
abcd2 <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/abcd2.txt")
abcd2[1]; abcd2[2]; abcd2[3]
abcd3 <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/abcd3.txt")
abcd3[1]; abcd3[2]; abcd3[3]; abcd3[4]; abcd3[5]
abcd4 <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/abcd4.txt")
abcd4[1]; abcd4[2]; abcd4[3]
abcd5 <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/abcd5.txt")
abcd5[1]; abcd5[2]; abcd5[3]
#
# I file potrei copiarmeli da rete usando il browser o con vari comandi di
# R. Ad es. con
# download.file("http://macosa.dima.unige.it/om/prg/R/agg/abcd1.txt","...")
# potrei copiare abcd1.txt come file ... del mio computer
#
# Vediamo come analizzare i diversi file. Il primo ha 1 riga di commenti, che
# salto con skip=1, e poi una riga in cui i dati sono separati da ",": lo carico
# e lo metto nella variabile dati1 (per comodità metto il nome in file)
file <- "http://macosa.dima.unige.it/om/prg/R/agg/abcd1.txt"
dati1 <- scan(file,skip=1,sep=",")
dati1
# Poi posso analizzare il file come fatto sopra per stud
barplot(dati1)
# Analogamente procedo per abcd2 tenendo conto che il separatore è ";"
file <- "http://macosa.dima.unige.it/om/prg/R/agg/abcd2.txt"
dati2 <- scan(file,skip=1,sep=";")
dati2
# Per abcd3 in cui il separatore è "a capo" basta che batta
file <- "http://macosa.dima.unige.it/om/prg/R/agg/abcd3.txt"
dati3 <- scan(file,skip=1)
dati3
# Lo stesso otterrei se il separatore fosse uno "spazio"
#
# Per abcd4 la cosa è più articolata:
file <- "http://macosa.dima.unige.it/om/prg/R/agg/abcd4.txt"
dati <- scan(file,skip=2,sep=";")
# I dati erano nella terza riga: ne ho saltato 2. I nomi sono nella
# seconda; vedendo l'help per scan vedo che posso leggere anche caratteri:
help(scan)
#
nomi <- scan(file,skip=1,sep=",",what="character",n=4)
# Quindi (come fatto all'inizio con "stud"):
names(dati) <- nomi; barplot(dati)
#
# Vediamo, ora, come confrontare distribuzioni di dati analoghe. Consideriamo
# i seguenti: gli studenti iscritti a due diversi corsi di studi per provenienza
# dalle nazioni A, B, C e D.
stud <- c(12, 72, 36, 40); stud2 <- c(37, 105, 60, 28)
# Il primo modo, più semplice, e di confrontarne la distribuzione con
# degli areogrammi circolari, in due finestre diverse, in modo da poterli
# confrontare. Basta usare il comando dev.new() per aprire una nuova finestra
# (posso spostare agendo col mouse sulla barra del titolo una finestra per vedere
# le altre.
pie(stud); dev.new(); pie(stud2)
# Con barplot posso vedere meglio il confronto tra una classe e l'altra
# dello stesso diagramma
dev.new(); barplot(stud); dev.new(); barplot(stud2)
# Se voglio confrontare le distribuzioni percentuali posso fare:
stud_p <- stud/sum(stud); stud2_p <- stud2/sum(stud2)
stud_p*100; stud2_p*100
dev.new(); barplot(stud_p*100); dev.new(); barplot(stud2_p*100)
# Sono tutte cose particolarmente semplici.
#
# Vediamo come rappresentare i due istogrammi sovrapposti in una stessa
# finestra. Abbiamo visto che come quota arrivano a 40 e qualcosa,
# quindi possiamo tracciarli tra 0 e 50, usando il comando ylim per
# scegliere la scala verticale, invece che farla scegliere automaticamente
# al programma. Cancelliamo prima tutte le finestre grafiche precedenti.
barplot(stud_p*100,angle=45,density=7,ylim=c(0,50))
# Con angle e density ho specificato come tracciare dei segmenti sulle
# colonne dell'istogramma. Traccio, tratteggiate, delle tacche.
abline(h=c(10,20,30,40,50),lty=3)
# L'altro istogramma lo traccio sulla stessa scala (add=TRUE), evidenziandone
# diversamente le colonne, con bordo rosso
barplot(stud2_p*100,add=TRUE,angle=135,density=14,border="red")
#
# R è in grado di rappresentare dati memorizzati in quasi ogni formato.
# Se essi sono in un file occorre esaminarli prima di decidere come analizzarli.
# Vediamo un altro esempio di uso di un file recuperato da rete.
boh <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/alimger.txt")
boh[1:11]
# Ho visualizzato un po' di righe del file, e ho capito di che cosa si tratta
# Metto i dati in una tabella, saltando le prime tre righe:
alimger <- read.table("http://macosa.dima.unige.it/om/prg/R/agg/alimger.txt",skip=3)
alimger
# rappresento i dati, che sono nella colonna V1, con un istogramma a barre
barplot(alimger$V1,horiz=TRUE)
# se voglio le scritte, che sono nella colonna V2, posso aggiungere:
barplot(alimger$V1,horiz=TRUE,names.arg=alimger$V2)
# non le vedo tutte perché lunghe; le "accorcio":
barplot(alimger$V1,horiz=TRUE,names.arg=substr(alimger$V2,1,2))
help(substr)
# o le metto orizzontali
barplot(alimger$V1,horiz=TRUE,names.arg=alimger$V2,las=1)
# aggiungo una griglia verticale, tratteggiata:
abline(v=c(20,40,60,80), lty=3)
#
# Volendo, rappresento i dati con un areogramma circolare:
pie(alimger$V1, alimger$V2)
#
#
# Abbiamo visto come analizzare distribuzioni in classi qualitative.
# Vediamo ora dati di tipo "QUANTITATIVO", ossia dati che esprimono grandezze
# quantitative (misure, valori monetari, …). Per ora consideriamo dati non
# già classificati.
# Partiamo rivisitando i dati sugli arrivi ad uno sportello già considerati
# nella prima e nella seconda attività del primo incontro. Per cambiare,
# li vediamo in un altro formato: questo t_arrivi2.txt
file <- "http://macosa.dima.unige.it/om/prg/R/agg/t_arrivi2.txt"
# se non sono in rete il file t_arrivi2.txt lo posso caricare dal pc: vedi
readLines(file)
# visto come è il file, lo carichiamo:
DATI <- scan(file, skip=1, sep=",")
length(DATI); summary(DATI); IQR(DATI)
# IQR è la distanza interquartile, tra 1º e 3º (ossia tra 25º e 75º percentile):
# l'ampiezza dell'intervallo in cui cade il 50% centrale dei dati.
# Se voglio, posso avere altri percentili. Ad es.:
quantile(DATI,c(0.05,0.1,0.25,0.5,0.75,0.9,0.95))
#
boxplot(DATI,horizontal=TRUE, col="yellow")
# R traccia solo il 25º e il 75º percentile evidenziando poi con dei pallini
# gli eventuali dati che precedono il 25º o seguono il 75º per più di 1.5
# volte la distanza interquartile, dando l'idea della (eventuale) dispersione
# dei dati. Volendo, con l'opzione range=0 si elimina questa visualizzazione:
dev.new()
boxplot(DATI, horizontal=TRUE, col="yellow", range=0)
# Non mettere "range" equivale a mettere range=1.5, altrimenti posso specificare
# un valore diverso (con lo stesso significato descritto sopra per 1.5)
boxplot(DATI,horizontal=TRUE, col="yellow", range=2)
# Rivediamo le rappresentazioni stem-and-leaf:
stem(DATI)
# possiamo vedere solo le quantità delle varie colonne:
stem(DATI,width=0)
# se i dati sono molti (ad es. il triplo), non si ha una buona rappresentazione:
stem(c(DATI,DATI,DATI))
# a meno che non la si ottenga aumentando il numero delle classi:
stem(c(DATI,DATI,DATI),scale=2)
# Con hist traccio l'istogramma; specificando "right=FALSE" gli intervalli
# vengono presi [.,.) - aperti a destra - come è meglio fare:
hist(DATI); dev.new(); hist(DATI, right=FALSE)
#
# hist fa una scelta automatica degli intervalli; comunque possiamo sceglierli
# noi, specificando come prendere le classi:
dev.new(); hist(DATI,seq(0,180,10),right=FALSE,col="yellow",xlab="",main="")
#
# Aggiungendo probability=TRUE vegono rappresentate sulla scala verticale non
# le frequenze assolute ma quelle relative, ovvero la densità di frequenza
# (freq. diviso ampiezza delle classi): help di seq
hist(DATI,seq(0,180,10),right=FALSE,col="yellow",xlab="",main="",probability=TRUE)
# Potrei scegliere liberamente le classi. In questo caso (ovviamente - perché?)
# viene rappresentata la densità di frequenza.
hist(DATI,c(0,10,20,30,50,75,100,140,180),right=FALSE,col="grey")
#
# Può essere comodo analizzare l'evoluzione di una sequenze di dati o
# confrontare diverse distribuzioni; facciamolo nel caso di dati inventati:
DATI1 <- DATI/2; DATI2 <- DATI+50; DATI3 <- DATI*1.5
boxplot(DATI1,DATI,DATI2,DATI3, range=0, col="yellow",horizontal=TRUE)
# Volendo posso aggiungere una griglia, piccola:
rug(seq(10,260,10))
# o più estesa:
rug(seq(25,250,25), ticksize=1, col="grey")
# Ovviamente, usando rug, potrei visualizzare anche altri percentili. Ecco, ad
# es., come visulizzare accanto al 25°,50° e 75° percentile, il 5° e il 95°:
boxplot(DATI, horizontal=TRUE, col="yellow", range=0)
rug(quantile(DATI,c(0.05,0.95)), ticksize=1, col="grey"); rug(seq(10,170,10))
#
# Confrontiamo il grafico delle frequenze con quello delle frequenze cumulate
# (empirical cumulative distribution function)
hist(DATI,seq(0,180,10),right=FALSE,col="yellow",xlab="",main="",probability=TRUE)
dev.new()
plot(ecdf(DATI))
# possiamo rimpicciolire i "punti"
plot(ecdf(DATI),cex=0.5)
# o eliminarli
plot(ecdf(DATI),cex=0)
# e tracciare dei segmenti verticali
plot(ecdf(DATI), cex=0, verticals=TRUE)
abline(h=0,v=0,col="blue")
# Ovviamente, abbiamo una curva con asintoto orizzontale alla quota 1 (100%)
#
# Consideriamo ora dati non raccolti in classi ma dotati di FREQUENZA
# Prima richiamiamo il comando rep(x,n), che costruisce una sequenza di n oggetti
# uguali ad x (replica x n volte)
rep(1/8, 3)
# Analizziamo questo file:
altf <- readLines("http://macosa.dima.unige.it/om/prg/R/altezzef.txt")
# Ricordiamo, basta: altf <- readLines("altezzef.txt") se lavoriamo non in rete e
# se abbiamo scelto .../R/agg come cartella di riferimento
altf[1:5]
# Il commento inziale conferma che si tratta di dati dotati di frequenza
# Li carico come tabella (non indico come sono separati - dovrei farlo se lo
# fossero da ";" o da ",", non da " ")
dati <- read.table("http://macosa.dima.unige.it/om/prg/R/altezzef.txt")
# Non ho aggiunto "skip=1" per saltare la prima riga in quanto essa inizia
# con "#" e il programma la interpreta automaticamente come un commento.
#
# Ne visulizzo le dimensioni (39 righe di 2 dati ciascuna)
dim(dati)
# Controlliamo quanti sono i dati. Le colonne, come visto, sono V1 e V2
sum(dati$V2)
# In tutto sono 998 valori (essendo una distribuzione per 1000 è giusto
# che siano circa 1000; essendo 39 intervalli sarebbe non facile, a
# causa degli arrotondamenti, che la somma fosse proprio 1000).
# È una tabella di dati e relative frequenze. Vediamo uno dei modi
# in cui possiamo analizzarla.
# Metto in una variabile i dati e le loro frequenze usando il comando rep.
dati2 <- rep(dati$V1,dati$V2)
# Col comando hist ne costruisco l'istogramma
hist(dati2,right=FALSE)
# Se voglio avere le frequenze percentuali aggiungo probability=TRUE
hist(dati2,right=FALSE,probability=TRUE)
# Se voglio classificare in intervalli a mio piacimento faccio:
hist(dati2,c(145,150,155,160,161,162,163,164,165,170,175,180,185),right=FALSE)
# o:
hist(dati2,c(seq(145,160,5),seq(161,174,1),seq(175,185,5)),right=FALSE)
# In questi casi (in cui gli intervalli non sono di eguale ampiezza) vengono
# visualizzate solo le frequenze percentuali.
hist(dati2,c(seq(145,185,1)),right=FALSE,probability=TRUE)
hist(dati2,c(seq(145,185,1)),right=FALSE)
#
# Vedremo più avanti come calcolare lo scarto quadratico medio
# e la approssimazione del profilo superiore degli istogrammi con curve.
#
#
# Consideriamo, ora, l'analisi di DATI CLASSIFICATI.
# Prendiamo come dati la distribuzione dell'età dei morti in Italia nel 1990.
# I dati sono in centinaia di persone: ad es. sono morte 25 centinaia di persone
# nella fascia 1-14 anni (cioè in [1,15): avevano compiuto 1 anno e non
# ancora i 15).
# Metto in una variabile (freq) le frequenze, in un'altra (interv) gli estremi
# degli intervalli, che sono uno di più.
freq <- c(46,25,58,186,870,1071,3124)
interv <- c( 0, 1,15, 25, 45, 65, 75, 100)
# Poi prendiamo a "scatola nera" le righe "viola"
totale <- sum(freq); fr_perc <- freq/totale*100; fr_perc
rap <- 1e5/totale; fr <- freq*rap; num_int <- length(freq)
dati <- c(seq(interv[1],interv[2],by=(interv[2]-interv[1])/fr[1]))
for(i in 2:num_int){dati<- c(dati,seq(interv[i],interv[i+1],by=(interv[i+1]-interv[i])/fr[i]))}
hist(dati,interv,col="yellow",axes=FALSE,xlab="",main="")
axis(1,pos=0,col="blue",label=TRUE, at=interv)
axis(2,pos=interv[1],col="blue",label=TRUE)
# Posso, poi, decidere se tracciare una griglia per le "y":
tacchey <- c(0.005,0.01,0.015,0.02,0.025)
abline(h = tacchey, col="grey60", lty=3)
# Volendo, nel caso le tacche coprissero parte dell'istogramma,
# puoi ritracciare questo senza colore con:
hist(dati,interv,axes=FALSE,xlab="",main="",ylim=c(0,max(tacchey)),add=TRUE)
#
# A voce spieghiamo i comandi, che per un docente sono abbastanza facili da
# comprendere, ma che non vale la spesa di spiegare agli alunni: si possono
# copiare e incollare, nei pochi casi in cui può servire ricorrere ad essi
# (trasformo la popolazione di dati classificati in intervalli - [0,1) [1,15)
# ... - in 100 mila dati singoli a cui applicare poi i comandi standard di R).
# Non c'è modo di ottenere rappresentazioni come queste con un foglio di calcolo.
# Poi su questi "dati" possiamo fare tutte le elaborazioni viste sopra, tenendo
# conto, nell'usare le uscite numeriche, che si tratta di valori approssimativi,
# elaborati supponendo che i dati abbiano esattamente l'andamento descritto dall'
# istogramma. Vedremo la prossima volta come caricare facilmente questi comandi.
#
summary(dati)
#
# Le cose messe al fuoco sono molte. Svilupperemo, rapidamente, ancora alcuni
# aspetti legati all'uso del generatore di numeri casuali e al foglio di calcolo.
# Riprenderemo questi aspetti e approfondiremo le statistiche "multivariate",
# ossia riferite a più variabili casuali, in parte del prossimo incontro.
# La gestione di esso e del successivo incontro sarà più agevole, una volta
# che avremo preso più confidenza con R.
#
#
# runif è il GENRATORE di numeri pseudocasuali
runif(1); runif(1); runif(1)
# vengono generati numeri con valori "casuali" tra 0 ed 1; da esso possono
# essere ottenuti altri valori casuali, ad esempio relativi al lancio di un
# dado ("floor" è la "parte intera"):
help(floor)
floor(runif(1)*6)+1; floor(runif(1)*6)+1; floor(runif(1)*6)+1
floor(runif(30)*6)+1
# runif(n) genera n numeri casuali
# il comando set.seed(k) ("stabilire il seme") con k numero naturale fa partire la
# successione "casuale" da un dato valore (con lo stesso k ho la stessa sequenza).
set.seed(23); x <- floor(runif(1000)*6)+1; mean(x)
set.seed(23); x <- floor(runif(1000)*6)+1; mean(x)
set.seed(4); x <- floor(runif(1000)*6)+1; mean(x)
set.seed(4); x <- floor(runif(1000)*6)+1; mean(x)
# Per scegliere "a caso" il seme possiamo scegliere come seme il tempo in secondi
# dall'avvio proc.time()[3]:
help(proc.time)
set.seed(proc.time()[3]); runif(1)
#
# Ecco la media del lancio di 1000 dadi (che si avvicina a 3.5):
x <- floor(runif(1000)*6)+1; mean(x)
# Vedremo la prossima volta come usare runif per fare altre simulazioni; è evidente,
# comunque, la sua semplicità d'uso rispetto all'analogo "strumento" presente nei
# fogli di calcolo.
#
# Vediamo, infine, la gestione di una TABELLA, come questa. Essa è memorizzata
# in un foglio di calcolo così, ad esempio nel formato standard csv.
# Per vederne il contenuto direttamente da R possiamo usare uno dei modi già visti:
tabella <- readLines("http://macosa.dima.unige.it/om/prg/R/agg/tabella.csv")
tabella[1:10]
# (uso solo readLines("tabella.csv") se ho usato Change Dir)
# Capito il formato, gli facciamo leggere la tabella:
dati <- read.table("http://macosa.dima.unige.it/om/prg/R/agg/tabella.csv", sep=";")
dati
# Saltiamo la prima riga, e abbiamo:
dati <- read.table("http://macosa.dima.unige.it/om/prg/R/agg/tabella.csv", sep=";", skip=1)
dati
# Volendola vedere in formato "foglio di calcolo" basta che battiamo:
edit(dati)
# Facciamo la media dei dati nella seconda colonna:
mean(dati$V2)
#
# Nota: è bene impostare i fogli di calcolo con la notazione dei numeri standard,
# ossia "inglese", usando il "punto" e non la "virgola" come separatore tra parte
# intera e parte frazionaria di un numero. Tutto il software scientifico (e i
# fogli di calcolo impostati in lingua inglese) non riconoscono questa notazione.
# Ecco che cosa può succedere:
dati <- read.table("http://macosa.dima.unige.it/om/prg/R/agg/tabella2.csv", sep=";", skip=1)
dati
mean(dati$V2)
#
# Alla prossima ...
#
# Torniamo indietro (<- clicca).
#