# # Vediamo come fare qualche altra analisi statistica di dati. # # Ecco un modo alternativo a quelli visti in g.htm e g_1.htm per costruire un # diagramma a torta ("pie" in inglese). Vediamo come rappresentare # le estensioni complessive delle zone di collina, montagna e pianura # del Piemonte espresse in decine di chilometri quadrati: # primo modo (semplicissimo): pie(c(770,1099,671)) # modo più sofisticato (aggiungo i nomi): piemonte <- c(770,1099,671) names(piemonte) <- c("col","mon","pia") pie(piemonte) # scelgo anche i colori: pie(piemonte,col=c("yellow","brown","green")) # # Se voglio mantenere il precedente diagramma posso fare il successivo # in una nuova finestra, che posso aprire col comando dev.new() # # La rappresentazione con un diagramma a barre ("barplot"), primo modo: barplot(piemonte) # Eccola riducendo lo spazio tra le colonne e mettendo una griglia barplot(piemonte,space=0) # Traccia linee orizzontali a certe quote, tratteggiate abline(h=c(200,400,600,800,1000),lty=3) # Eccola aggiungendo colori: barplot(piemonte,space=0,col=c("yellow","brown","green")) abline(h=c(200,400,600,800,1000),lty=3) # La distribuzione percentuale (sum fa la somma): piemonte/sum(piemonte)*100 col mon pia 30.31496 43.26772 26.41732 # Il suo arrotondamento agli interi round(piemonte/sum(piemonte)*100) col mon pia 30 43 26 # e ai decimi round(piemonte/sum(piemonte)*100, 1) col mon pia 30.3 43.3 26.4 # Ecco il diagramma con le percentuali: colori <- c("yellow","brown","green") barplot(piemonte/sum(piemonte)*100,space=0,col=colori) abline(h=c(10,20,30,40),lty=3) # # # Ecco le altezze di 3 ragazze di scuola media e l'esito di alcuni calcoli: x <- c(145,152,147) c( min(x), max(x), sum(x) ) [1] 145 152 444 c( length(x), sum(x)/length(x), mean(x) ) [1] 3 148 148 # si sono trovati il minimo, il massimo e la somma dei dati; poi si sono # scritti quanti sono i dati (la "lunghezza" di x), la loro media, espressa # come somma divisa per il numero dei dati o direttamente. # # Vediamo la analoga analisi fatta su più dati: alu <- c(146,158,152,140,157,147,161,147,149,154,155,153,155,156,150,153,152,145) c( min(alu), max(alu), length(alu) ) [1] 140 161 18 # Mettiamo i dati in ordine ("sort" in inglese significa "sorta", "specie" ma # anche "ordinamento") sort(alu) [1] 140 145 146 147 147 149 150 152 152 153 153 154 155 155 156 157 158 161 c( min(alu), max(alu), length(alu), mean(alu), median(alu) ) [1] 140.0000 161.0000 18.0000 151.6667 152.5000 # Si noti che quando si stampano più dati usando c(...) essi vengono tutti # scritti con cifre dopo il "." se ce n'è qualcuno che ne ha. # Alla fine abbiamo stampato anche il valore della mediana ossia il valore # che sta al centro dell'elenco. # # Ecco l'istogramma. Il programma sceglie automaticamente le classi in cui # classificare i dati. Il comando è hist. hist(alu, right=FALSE) # L'espressione "right=FALSE" serve per comunicare al programma che vogliamo che gli # intervalli siano del tipo [ , ), in modo che ad esempio 150 non sia classificato # tra 145 e 150 ma tra 150 e 155: 150 non sta in [145,150) ma sta in [150,155). # Vengono rappresentate le frequenze delle varie classi in cui sono stati # suddivisi i dati. Se aggiungiamo probability=TRUE otteniamo la stessa figura # (vedi immagine di sopra, a destra), ma cambia la scala verticale: sono # rappresentate non le frequenze ma le frequenze relative unitarie (o densità di # frequenza), ossia le frequenze relative divise per l'ampiezza degli intervalli. hist(alu, right=FALSE, probability=TRUE, col="yellow"); abline(h=0.06,lty=3) # Tra 150 e 155 c'erano 6 dati; i dati sono 18, quindi la frequenza relativa è # 6/18 = 0.333 (cioè 33.333 %), e quella unitaria è 6/18/5 = 0.06666 = 6.7%. # # Come faccio a trovare frequenze assolute e quelle percentuali (ossia le # frequenze relative in forma percentuale) facilmente? Ecco: hist(alu,right=FALSE)$counts [1] 1 5 6 5 1 hist(alu,right=FALSE)$counts/length(alu)*100 [1] 5.555556 27.777778 33.333333 27.777778 5.555556 # # # Volendo posso decidere una sequenza di classi in cui classificare i dati. # Basta che usi seq(A,B,h) per classificarli tra A e B in classi ampie h. # Utilizziamo i dati alu già impiegati in precedenza. alu <- c(146,158,152,140,157,147,161,147,149,154,155,153,155,156,150,153,152,145) hist(alu, right=FALSE, seq(140,161,3)) # Per scritte sugli assi più piccole potrei mettere hist( , cex.axis=0.88) # # Con summary ho una sintesi delle informazioni numeriche sui dati. # Vengono stampati anche il 1° e il 3° quartile, cioè i dati che stanno # a metà della prima metà e a metà della seconda metà dei dati; la # mediana sarebbe il 2° quartile. summary(alu) Min. 1st Qu. Median Mean 3rd Qu. Max. 140.0 147.5 152.5 151.7 155.0 161.0 # Volendo, se hai caricato il file seguente: # source("http://macosa.dima.unige.it/r.R") # hai gli stessi esiti e una loro rappresentazione grafica con: statistiche(alu) Min. 1st Qu. Median Mean 3rd Qu. Max. 140.0 147.5 152.5 151.7 155.0 161.0 I pallini sono il 5° e il 95° percentile # in cui il "box" azzuro rappresenta i dati dal 1° al 3° quartile e i # pallini rappresentano i valori che corrispondono a dove starebbero il # 5º e il 95º dato se i dati, distribuiti in modo simile, fossero 100. # # Se abbiamo più dati da analizzare li mettiamo su più righe: alu2 <- c(146,158,152,140,157,147,161,147,149,154,155,153, 155,156,150,153,152,145,158,140,147,149,155,156,153,145, 141,156,158,151,161,148,157,159,160,162,149,163,152,149, 157,155) # Ecco la loro analisi statistica: hist(alu2,right=FALSE) hist(alu2,right=FALSE,seq(140,164,3)) summary(alu2) Min. 1st Qu. Median Mean 3rd Qu. Max. 140.0 149.0 153.0 152.6 157.0 163.0 # Posso anche usare come sopra statistiche() fissando la stessa scala per alu e alu2, # col comando boxAB: BF=4;HF=1.2 boxAB=c(140,163); statistiche(alu) Min. 1st Qu. Median Mean 3rd Qu. Max. 140.0 147.5 152.5 151.7 155.0 161.0 I pallini sono il 5° e il 95° percentile boxAB=c(140,163); statistiche(alu2) Min. 1st Qu. Median Mean 3rd Qu. Max. 140.0 149.0 153.0 152.6 157.0 163.0 I pallini sono il 5° e il 95° percentile # # Ecco come possiamo confrontare graficamente i due insiemi di dati: hist(alu,seq(140,164,3),angle=45,density=7,probability=TRUE) hist(alu2,seq(140,164,3),probability=TRUE,add=TRUE,angle=135,density=10) # Ho aggiunto add=TRUE al comando per tracciare il 2° istogramma. # La rappresentazione a destra è stata ottenuta con questi comandi. # # Volendo posso introdurre una quantità variabile di input direttamente # da tastiera mediante il comando scan(file="",n=...) in cui al posto # di ... occorre mettere i dati da leggere. Due esempi. Il primo (faccio # la somma di 3 numeri): numeri <- scan(file="", n=3); sum(numeri) 1: 3 2: 7 3: 5 Read 3 items [1] 15 # Il secondo (introduco un numero N e, poi, una quantità N di dati, di # cui viene calcolata la media): nalu <- scan(file="",n=1); altalu <- scan(file="",n=nalu); mean(altalu) 1: 7 Read 1 item 1: 154 2: 156 3: 163 4: 170 5: 161 6: 156 7: 159 Read 7 items [1] 159.8571 # # # Un altro esempio di come i comandi possono essere spezzati su più # righe. Ecco le lunghezze di molte fave (ossia semi di fava) raccolte # da una classe di alunni di 12 anni. Occorre copiare e incollare tutte # le righe, da "fave <- c(" alla riga finale ")". fave <- c( 1.35,1.65,1.80,1.40,1.65,1.80,1.40,1.65,1.85,1.40,1.65,1.85,1.50,1.65,1.90, 1.50,1.65,1.90,1.50,1.65,1.90,1.50,1.70,1.90,1.50,1.70,1.90,1.50,1.70,2.25, 1.55,1.70,1.55,1.70,1.55,1.70,1.60,1.70,1.60,1.75,1.60,1.75,1.60,1.80,1.60, 1.80,1.60,1.80,1.60,1.80,1.00,1.55,1.70,1.75,1.30,1.55,1.70,1.75,1.40,1.60, 1.70,1.75,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40, 1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80, 1.45,1.60,1.70,1.80,1.50,1.60,1.70,1.80,1.50,1.60,1.70,1.85,1.50,1.60,1.70, 1.85,1.50,1.60,1.75,1.90,1.50,1.60,1.75,1.90,1.50,1.65,1.75,1.90,1.55,1.65, 1.75,1.95,1.55,1.65,1.75,2.00,1.55,1.65,1.75,2.30,1.35,1.65,1.80,1.40,1.65, 1.80,1.40,1.65,1.85,1.40,1.65,1.85,1.50,1.65,1.90,1.50,1.65,1.90,1.50,1.65, 1.90,1.50,1.70,1.90,1.50,1.70,1.90,1.50,1.70,2.25,1.55,1.70,1.55,1.70,1.55, 1.70,1.60,1.70,1.60,1.75,1.60,1.75,1.60,1.80,1.60,1.80,1.60,1.80,1.60,1.80, 1.00,1.55,1.70,1.75,1.30,1.55,1.70,1.75,1.40,1.60,1.70,1.75,1.40,1.60,1.70, 1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60, 1.70,1.80,1.40,1.60,1.70,1.80,1.40,1.60,1.70,1.80,1.45,1.60,1.70,1.80,1.50, 1.60,1.70,1.80,1.50,1.60,1.70,1.85,1.50,1.60,1.70,1.85,1.50,1.60,1.75,1.90, 1.50,1.60,1.75,1.90,1.50,1.65,1.75,1.90,1.55,1.65,1.75,1.95,1.55,1.65,1.75, 2.00,1.55,1.65,1.75,2.30 ) # e il loro istogramma hist(fave,right=FALSE) # e quello delle frequenze relative, con aggiunte delle righe # tratteggiate orizzontali: hist(fave,right=FALSE,probability=TRUE) abline(h=c(1/2, 1, 1.5, 2, 2.5), lty=3) # (il 2.5% dei dati sta tra 1.6 e 1.7 e un altro 2.5% tra 1.7 e 1.8) # La sintesi delle informazioni: summary(fave) # ovvero, caricato: source("http://macosa.dima.unige.it/r.R") statistiche(fave) Min. 1st Qu. Median Mean 3rd Qu. Max. 1.000 1.550 1.650 1.659 1.750 2.300 # # Volendo posso scegliere io gli estremi degli intervalli. Se non li scelgo # di uguale ampiezza viene tracciato l'istogramma con le frequenze relative. hist(fave,c(1,1.3,1.5,1.7,1.9,2.3)) # # I dati possono essere copiati da rete e incollati in R o possono essere # caricati da rete. Un esempio (il file"t-sec.txt" nella cartella "R"). # Leggo alcune righe (ad es. 4) di un file (ne leggo poche in quanto il file può # essere lunghissimo, ad es. contenere qualche migliaia di dati): readLines("http://macosa.dima.unige.it/R/t-sec.txt",n=4) #[1] "# misure (troncate ai centesimi di sec) della durata di 1 s cronometrata a mano" #[2] "111" #[3] "103" #[4] "109" # Vedo che sono singoli dati interi e che c'è 1 riga da saltare. Lo faccio col # comando skip. I dati sono separati da "aCapo"; li leggo col comando "scan" # (se fossero separati da, ad es., ";" aggiungerei sep=";" ) dati <- scan("http://macosa.dima.unige.it/R/t-sec.txt", skip=1) # Read 47 items # In alternativa, posso esaminare il file da un browser (programma per navigare # in Internet), copiarlo e poi usare scan("clipboard", skip=1) # Il comando str visualizza la struttura di oggetti. Lo uso per esaminare dati str(dati) # num [1:47] 111 103 109 97 99 110 99 103 109 106 ... # sono 47 dati. I dati sono troncati agli interi. Per fare valutazioni corrette # devo avere dati arrotondati: aggiungo 1/2 DATI <- dati+1/2 summary(DATI) # Min. 1st Qu. Median Mean 3rd Qu. Max. # 68.50 96.50 98.50 99.86 108.00 129.50 hist(DATI,right=FALSE,seq(65,135,10)) # visto il grafico, scelgo la scala verticale in modo da vedere la quota 25 hist(DATI,right=FALSE,seq(65,135,10),ylim=c(0,25)) abline(h=c(5,10,15,20,25), lty=3) # # Un'ultima osservazione. # Se dispongo di dati già classificati, in intervalli di diversa ampiezza, # posso riccorrere al comando istoclas presente nel file richiamabile con: source("http://macosa.dima.unige.it/r.R") # Vedi qui per un esempio: # # Possiamo generare dei numeri a caso. Ad es. runif(4, min=1, max=10) # [posso ottenere 9.139260 3.747637 3.847650 2.825768] # genera 4 numeri casuali distribuiti "uniformemente" nell'intervallo (1,10). # Invece: trunc(runif(10, min=1, max=7)) # [ottengo ad esempio 2 4 1 4 5 3 4 6 2 5; "trunc" tronca un numero agli interi] # genera 10 numeri casuali interi compresi tra 1 e 6 # Vediamo il lancio, 10 mila volte, di due dadi "equi": n <- 10000; dadi <- trunc(runif(n, min=1, max=7))+trunc(runif(n, min=1, max=7)) # I valori sono stati messi in "dadi". Eccone alcuni: dadi[1]; dadi[1:10] # [1] 7 # [1] 7 9 9 7 8 6 4 10 7 7 hist(dadi, seq(1.5,12.5,1),probability=TRUE,col="yellow") abline(h=c(0.05,0.1,0.15), lty=3) # Le uscite possibili sono 6*6 = 36; avendo supposto il dado "equo" sono tutte # equiprobabili, ossia con probabilità 1/36. Ad es. 3 può uscire in 2 casi: # 1 e 2, 2 ed 1; quindi la probabilità che esca è 2*1/36, ossia 0.0555 ; e con # 10 mila lanci ho ottenuto una frequenza relativa è vicina a 0.0555 # L'istogramma a destra è stato ottenuto col comando Istogramma: noClassi=1; Istogramma(dadi, 1.5,12.5, 1) # Frequenze e frequenze percentuali: #277, 554, 857, 1084, 1442, 1663, 1337, 1126, 862, 532, 266 #2.77,5.54,8.57,10.84,14.42,16.63,13.37,11.26,8.62,5.32,2.66 # La freq. relativa di 3 è 5.54% (ossia 0.0554). # # Posso trovare anche probabilità relative a fenomeni non facilmente studiabili in # modo teorico. Vediamo come studiare le uscite del lancio di una coppia di dadi # equi, assumendo come uscita la differenza dei valori che escono, che può # andare da 0 a 5. Ci aspettiamo che 5 sia l'uscita meno probabile (viene solo se # escono un 6 ed un 1), ma non è facile valutare gli altri casi. # Per fare la differenza posso prendere il valore assoluto della sottrazione della # seconda uscita dalla prima. n=1e7; dadi=abs(trunc(runif(n, min=1,max=7))-trunc(runif(n, min=1,max=7))) Istogramma(dadi,-0.5,5.5,1) Frequenze e frequenze percentuali: 1666379, 2778367, 2221573, 1665591, 1111816, 556274 16.66,27.78,22.22,16.66,11.12,5.56 # Per n = 10^7 (così come per n = 10^6) ottengo che l'uscita più frequente è 1 # (27.78% di probabilità). # A questo punto posso provare a studiare il problema anche teoricamente ... # # Questo era un esempio di simulazione. Per valutare la precisione della stima si # può ricorrere al seguente comando: Pr <- function(n) {f <- 0; for (i in 1:n) f <- f + ifelse(Evento(),1,0); fr <- f/n; S <- sqrt(fr*(1-fr)/(n-1)); cat(fr, "+/-", 3*S,'\n') } # Introdotte queste righe e messo in Evento() una rappresentazione di un fenomeno # di cui si vuole valutare la probabilità che valga 1 se esso si verifica e 0 # altrimenti, valutiamo, riferendoci all'esempio precedente, la probabilità che # lanciando una coppia di dadi equi la loro differenza sia 1. Evento <- function() abs(floor(runif(1)*6)-floor(runif(1)*6)) == 1 Pr(1e4) #0.2825 +/- 0.01350713 Pr(1e5) #0.27771 +/- 0.004248885 Pr(1e6) #0.277448 +/- 0.001343219 Pr(1e7) #0.2777537 +/- 0.000424907 # Possiamo dedurre che la probabilità è 0.2825±0.0135, anzi 0.2777537±0.000425, # ovvero 0.27775±0.00043. Andando avanti (se ho tempo ) posso trovare valutazioni # più precise (queste valutazioni non sono certe, ma sono probabili al 99.7%). # Una valutazione più complessa: # qual è la probabilità, pescando 10 carte da un mazzo da 40, di ottenere almeno un # tris (cioè 3 o 4 carte dello stesso valore)? Evento <- function() { ca <- array(rep(0,40),dim=c(4,10)); nu <- array(rep(0,10),dim=10) tris <- 0; for (i in 1:10) {ripeti <- 1; while(ripeti==1) {seme <- floor(runif(1)*4)+1; valore <- floor(runif(1)*10)+1; if (ca[seme,valore]==0) ripeti <- 0 } ca[seme,valore] <- 1; nu[valore] <- nu[valore]+1; for(valore in 1:10) if (nu[valore]==3 | nu[valore]==4) {tris <- 1; valore <- 10} }; ifelse(tris==1,1,0) } # Descritto l'evento, valutiamo la probabilità: Pr(1e3) #0.409 +/- 0.04666528 Pr(1e4) #0.3765 +/- 0.01453596 Pr(1e5) #0.38396 +/- 0.004613929 Pr(1e6) # Per questa uscita occorre aspettare qualche minuto #0.384978 +/- 0.001459771 # La probabilità cercata è 0.3850±0.0015, ovvero (38.50±0.15)%. # # Volendo generare P uscite casuali intere diverse tra loro che cadono tra M,M+1,...,N # (N = P-1) uso sample(M:N). Esempi sample(14:20) #14 16 20 15 17 19 18 # La cosa è utile in molte situazioni. Es.: un generico ordine in cui posso estrarre # le tredici carte di cuori: sample(1:13) #4 9 12 3 13 7 10 11 5 6 2 1 8 # Posso estrarre (del tutto a caso) 4 numeri tra 1 e 13: sample(1:13, 4) #6 10 11 12 # # I numeri generati con runif e sample appaiono distribuiti in modo del tutto casuale. # In realtà posso ottenere la stessa sequenza di numeri se batto set.seed(n) ("poni # il seme eguale a n") con il medesimo numero intero n. Esempio: sample(5:10); set.seed(273); sample(5:10); set.seed(273); sample(5:10) # 6 10 7 5 8 9 # 6 7 5 9 10 8 # 6 7 5 9 10 8 set.seed(273); sample(5:10,3) # 6 7 5 # La cosa è ad es. utile per controllare dei procedimenti che si vogliono realizzare # in cui si impiega il generatore di numeri casuali. # Per approfondimenti su come "è fatto" il generatore di numeri casuali vedi qui. # # Con questi comandi sono stati realizzati anche Codifica(seme,testo) e Decodifica( , ) # per realizzare codifiche "segrete"; occorre mettere in "seme" un numero intero. # Ecco qualche esempio (vedi qui se vuoi vedere come sono stati realizzati): frase = "Ci vediamo alle 17:30 davanti a casa tua. Ciao!" Codifica(314, frase) # "\"m%rfgmcew%c~~f%5(*-?%gcrcdqm%c%vcac%q|c6%\"mcw/" Decodifica(314, "\"m%rfgmcew%c~~f%5(*-?%gcrcdqm%c%vcac%q|c6%\"mcw/") # "Ci vediamo alle 17:30 davanti a casa tua. Ciao!" frase2 = "195/3 = 65" Codifica(-10, frase2) # "F!I74TXTSI" Decodifica(-10, "F!I74TXTSI") # "195/3 = 65" # # # NOTA tecnica su come "salvare" i grafici e i comandi usati. # Per memorizzare un'immagine posso cliccare su di essa e dal menu che si # apre salvarla come BitMap e poi o # -- incollarla direttamente in un documento del programma con cui si sta # scrivendo (OpenOffice, LibreOffice, Wordpad, Word, ...) o # -- incollarla in Paint (o un altro programma di grafica), aggiungere # quello che si vuole, spostarla e ridurre i margini, e poi o # - copiarla nell'eventuale documento in cui si sta scrivendo o, # - da Paint, salvarla in formato PNG o GIF, che occupano poco spazio. # Per memorizzare i comandi usati posso, dal menu File, azionare il comando # Save History che consente di salvare in un file (in formato testo) tutti # i comandi eseguiti (devo dare un nome al file; posso lasciargli la # estensione "Rhistory" o mettere l'estensione "txt"). Poi, dopo, posso # [ o caricare automaticamente i comandi dal file usando il comando Load # History, che può essere usato anche per richiamare automaticamente # comandi caricati in altri file, o ] # aprire il file con BloccoNote (NotePad) o un altro editor e copiare e # incollare in R la sequenza di comandi che voglio. # # Per un formulario completo di tutti i comandi vai qui. # # Per usi più avanzati (scarto quadratico medio, gaussiana, altre distribuzioni, # tabelle di contingenza, ) vai qui. Ricordiamo solo che accanto ai comandi presenti # nel software var(dati) e sd(dati), che calcolano varianza e deviazione standard (o # scarto quadratico medio) "sperimentali", sono presenti, se hai caricato: # source("http://macosa.dima.unige.it/r.R") # i comandi Var(dati), Sd(dati) e SdM(dati) che calcolano la varianza, la deviazione # standard (o scarto quadratico medio) "teorici" e la deviazione standard della media. #