#
#   - - -   INCONTRO N. 3 - Statistica/probabilità (continuazione)
#                            e (avvio alla) Geometria                - - -
#
# Completiamo l'attività su statistica e probabilità, vediamo qualche esempio
# proposto da chi segue "a distanza" e cerchiamo di avviare il tema geometria.
#
# Vediamo come ottenere alcune animazioni. Vediamo prima come si può usare
# il generatore di numeri casuali per studiare un dado equo:
n <- 50; r <- runif(n); for (i in 1:10) {n <- n*2; r <- c(r,runif(n));
                                         hist(r,seq(0,1,1/6))}
# Facciamo qualche variante: coloriamo di grigio le colonne; stampiamo i valori
# di n, visulizziamo le frequenze relative, e ... introduciamo degli intervalli
# di tempo di 2 secondi tra un'immagine e la successiva. Lo possiamo fare con
# la seguente riga, che possiamo usare a scatola nera (vediamo, poi, come
# "caricarla" da rete)
tic <- function(x) {sec <-proc.time()[3]; while(proc.time()[3] < sec+x) sec<-sec}
n <- 50; r <- runif(n); for (i in 1:10) { tic(2); n <- n*2; r <- c(r,runif(n));
                  hist(r,seq(0,1,1/6), probability=TRUE, col="grey"); print(n) }
# All'aumentare del numero delle prove l'istogramma tende ad appiattirsi.
#
# Un altro esempio, a scatola nera (ma che da' un'idea di come si possono fare
# cose di "matematica" - relativamente - divertenti)
for (k in 0:20) { tic(0.3);
  plot(c(0,100),c(0,100),type="n",xlab="", ylab="", asp=1);
  symbols(50,80, circles=5, inches=FALSE, add=TRUE); #testa,occhi
  symbols(52,82, circles=1, inches=FALSE, add=TRUE);
  lines(c(55,50,50),c(0,0,75)); tic(0.3);  #corpo con gambe allineate
  plot(c(0,100),c(0,100),type="n",xlab="", ylab="", asp=1);
  symbols(50,80, circles=5, inches=FALSE, add=TRUE);
  symbols(52,82, circles=1, inches=FALSE, add=TRUE);
  # corpo con gambe non allineate (prima le x e poi le y)
  lines(c(35,30,50,70,75,70,50,50,35,50,65),
        c(1.5,2,40,2,2.5, 2,40,75,45,75,45)) }
#
# La riga "tic <- ..." definisci una funzione a cui si può dare come input un
# numero di secondi di attesa (in questo caso ne abbiamo dato 2). La spegazione
# della funzione è semplice, ma sulla sua ideazione non vale la pena perdersi ogni
# volta. Vediamo come possiamo caricarla da rete. Cose analoghe le possiamo fare
# per ogni argomento. Possiamo caricare una serie di comandi da usare con i nostri
# alunni e farli richiamare quando serve. Vediamo come.
#   [ lavorando non in rete conviene procedere così: con Change Dir
#     scegliere .../R 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/ - ad es. al posto del comando
#     seguente considerare:  source("tic.txt") ]
#
# Proviamo a introdurre:
source("http://macosa.dima.unige.it/om/prg/r/tic.txt")
# Abbiamo visto che cosa viene stampato. Vediamo il file. Possiamo usare vai
# comani per farlo. Un modo particolarmente semplice è ricorrere a Display file
# dal menu file e mettere nella finestra di dialogo che si apre
# http://macosa.dima.unige.it/om/prg/r/tic.txt.  Vedete il contenuto del file:
#  # commento stampato automaticamente
#  print(c(
#    "viene caricata la funzione tic(...) tale che se metto",
#    "tic(x) vengono attesi x secondi"
#    ), quote=FALSE)
#  # funzione definita
#  tic <- function(x) {sec <-proc.time()[3]; while(proc.time()[3] < sec+x) sec<-sec}
# Tra  print(c(  e  ), quote=FALSE)  possiamo inserire tutte le righe di commento
# che vogliamo. Nelle righe successive scriviamo i comandi che gli alunni trovano
# caricati automaticamente (se vogliamo, poi, possiamo fare degli altri esempi).
# Questo consente diversi livelli d'uso del programma. Uno può preparare una banca
# di comandi che si possono usare a scatola nera, consentendo un uso facile del
# programma anche a livello di scuola di base!!!
# Nota: se non mi ricordo "source" posso ricorrere al primo comando del menu
# a cascata in alto a sinistra, File (vedi)
#
# Ulteriore esempio (che fa?):
freq <- c(2,3,4,5); interv <- c(1,5,10,20,50)
source("http://macosa.dima.unige.it/om/prg/r/daticlas.txt")
griglia <- seq(0.01,0.04,0.01); abline(h = griglia, lty=3)
# Si tratta del caricamento automatico (a "scatola nera", come fosse un comando)
# del procedimento per analizzare dati classificati visto la vostra scorsa, qui).
#
# Abbiamo visto come può essere caricata una procedura complessa.
#
# Finora abbiamo considerato singole variabili casuali. Passiamo ora, rapidamente,
# a considerare SISTEMI di VARIABILI CASUALI.
#
# Esaminiamo tre file, memorizzati in diversi formati, e analizziamoli:
readLines("http://macosa.dima.unige.it/om/prg/r/battit1.txt",n=12)
readLines("http://macosa.dima.unige.it/om/prg/r/battit2.txt",n=11)
readLines("http://macosa.dima.unige.it/om/prg/r/battit3.txt",n=11)
# Consideriamo prima il secondo formato, che sembra "standard"
# È una tabella, proviamo ad esaminarne due righe, dopo i commenti
# (salto le 8 righe di commenti)
read.table("http://macosa.dima.unige.it/om/prg/r/battit2.txt",skip=8, nrows=2)
# OK. Leggo "tutta" la tabella e la metto, ad es., nella variabile  dati
# Poiché le righe da saltare iniziano con # (simbolo che in R indica un commento)
# potrei battere  read.table("http://macosa.dima.unige.it/om/prg/r/battit2.txt")
# senza skip=8 poichè le righe verrebbero saltate automaticamente.
dati <- read.table("http://macosa.dima.unige.it/om/prg/r/battit2.txt",skip=8)
# In questo caso, come posso controllare con nrow, sono poche righe di dati
nrow(dati)
# e posso visualizzarle tutte:
dati
# Con edit li avrei anche visti in formato "foglio di calcolo":
edit(dati)
# ma, in generale, conviene più in breve con str (struttura)
str(dati)
# Per un'idea rapida del file posso usare (come già fatto nel caso univariato) "summary"
summary(dati)
# Nota: il programma "ha capito", per come avevamo inserito i dati (qualche spazio
# iniziale nella prima riga) che 01, 02, ... erano i numeri delle righe, non dati.
#
# Prima di andare avanti vediamo che cosa si poteva fare negli altri casi.
read.table("http://macosa.dima.unige.it/om/prg/r/battit3.txt",skip=8, nrows=2)
# Vedo che mi prende BatPrima BatDopo ... come riga di dati.
# Allora la salto:
read.table("http://macosa.dima.unige.it/om/prg/r/battit3.txt",skip=9, nrows=2)
dati2 <- read.table("http://macosa.dima.unige.it/om/prg/r/battit3.txt",skip=9)
summary(dati2)
# Nel primo caso capisco che devo saltare 10 righe:
read.table("http://macosa.dima.unige.it/om/prg/r/battit1.txt",skip=10, nrows=2)
# Vedo che mi mette insieme più dati: devo spiegare sono separati da ";"
read.table("http://macosa.dima.unige.it/om/prg/r/battit1.txt",skip=10,sep=";",nrows=2)
# OK. A questo punto proseguo come sopra:
dati3 <- read.table("http://macosa.dima.unige.it/om/prg/r/battit1.txt",skip=10,sep=";")
summary(dati3)
# Analizziamo, dunque, questa tabella di dati
str(dati)
# Posso esaminare le singole variabili. Vediamo come:
altezze <- dati$Alt
summary(altezze)
# Se avessi usato dati2 avrei dovuto chiamare V6 le altezze
altezze2 <- dati2$V6; summary(altezze2)
#
hist(altezze)
# Ma ho due sottopopolazioni diverse per altezza:
datiF <- subset(dati,dati$Sesso==2); datiM <- subset(dati,dati$Sesso==1)
nrow(datiF); nrow(datiM)
str(datiF); str(datiM)
summary(datiF); summary(datiM)
hist(datiF$Alt,probability=TRUE); dev.new(); hist(datiM$Alt,probability=TRUE)
# Per confrontare i dati può essere comodo averli su una stessa scala
# (per abbreviare, memorizzo alcuni dati con nuove variabili):
M <- datiM$Alt; F <- datiF$Alt; int <- seq(150,200,5)
hist(F,int,right=FALSE,angle=45,density=14,xlab="",main="",probability=TRUE,ylim=c(0,0.06),col="red")
hist(M,int,right=FALSE,add=TRUE,angle=135,density=7,xlab="",main="",probability=TRUE,col="blue")
#
# i coefficienti di correlazione tra le diverse variabili (vedi)
cor(dati)
# Si vede che c'è una forte correlazione tra altezza e peso, però se mi restringo
# ad uno dei due sessi:
cor(datiM)
# vedo che questa scende decisamente. Elaborare dati statistici non è banale
# (spesso, purtroppo, le statistiche elaborate dai medici costituiscono dei buoni
# esempi di come non andrebbero fatte). Rappresentiamo graficamente (pesi,altezze)
# dei due sessi. Ricordiamo dove variano:
summary(dati$Peso); summary(dati$Alt)
# Scelgo di conseguenza la scala e traccio i due insiemi di punti
# Prima traccio tutti i punti:
plot(c(150,200),c(40,100),type="n",xlab="", ylab="")
points(dati$Alt,dati$Peso,col="blue")
# che si dispongono in forma allungata, e poi quelli femminili:
points(datiF$Alt,datiF$Peso,col="red")
#
# Finiamo questa rapida panoramica esaminando la retta di regressione
# tra peso e altezza nel caso dei maschi:
x <- datiM$Alt; y <- datiM$Peso; plot(x,y)
# la retta di regressione si ottiene così:
lm(y ~ x)
# Volendo una scrittura migliore posso fare:
mod = lm(y ~ x); mod$coefficients; abline(mod$coefficients)
# y = -62.70275 + 0.74876*x. Verifichiamo:
f <- function(x) -62.70275 + 0.74876*x
plot(f,150,200,add=TRUE,col="red")
#
# OK. Ci fermiamo qui. Il discorso si farebbe troppo lungo.
# Vediamo solo come, nel caso delle altezze (non nel caso dei pesi!!!)
# potremmo approssimare i dati con una gaussiana:
hist(datiM$Alt,seq(150,210,4),probability=TRUE,,col="yellow")
m <- mean(datiM$Alt); s <- sd(datiM$Alt); m; s
z <- function(x) dnorm(x,mean=m,sd=s)
plot(z,150,210,add=TRUE,col="blue")
#
# Sulla gaussiana, vediamo come si possono verificare/congetturare
# facilmente alcune proprietà:
media <- 3; stdev <- 5
z <- function(x) dnorm(x, mean=media, sd=stdev)
integrate(z,-Inf,Inf)
integrate(z,3,Inf)
integrate(z, media-stdev, media+stdev)
integrate(z, media-3*stdev, media+3*stdev)
media <- 9; stdev <- 5
integrate(z,9,Inf)
integrate(z, media-3*stdev, media+3*stdev)
#
# E vediamo un ultimo esempio, legato ad un ampio campione di dati.
readLines("http://macosa.dima.unige.it/om/prg/stf/leva.tab",n=11)
dati <- read.delim("http://macosa.dima.unige.it/om/prg/stf/leva.tab",header=FALSE,skip=9,sep =";")
summary(dati)# Aggiungo 1/2 ai dati troncati agli interi, in modo da non dover poi
# aggiungere 1/2 alle medie
dati$V1 <- dati$V1+1/2; dati$V2 <- dati$V2+1/2; dati$V4 <- dati$V4+1/2
summary(dati)
# quanti sono i dati:
length(dati$V1)
# i coefficienti di correlazione tra le diverse variabili
cor(dati)
# Vediamo, ad es., la relazione tra V4 e V2 (statura-100 e peso):
plot(c(40,100),c(45,130),type="n",xlab="", ylab="")
abline(v=seq(40,100,10),h=seq(50,130,10),lty=3)
# Con densCols posso colorare diversamente le uscite in base alle frequenze
dcols <- densCols(dati$V4,dati$V2)
points(dati$V4,dati$V2,col=dcols,pch=".", cex=5)
# cex dimensiona i punti
mod = lm(dati$V2 ~ dati$V4); mod$coefficients; abline(mod$coefficients,col="red")
lm(dati$V4 ~ dati$V2)$coefficients
#  x = a y + b  ->  ay = x-b  ->  y = (x-b)/a
g <- function(x) (x-55.4836084)/0.2738877
plot(g,40,100,add=TRUE,col="green4")
#
# FINE delle considerazioni statistico/probabilistiche
#
#
# Nel poco tempo che ci rimane avviamo alcune considerazioni di GEOMETRIA
#
# Esempio introduttivo:
ita <- read.table("http://macosa.dima.unige.it/om/prg/R/rmacosa/penisola.txt",header=FALSE,sep =",")
plot(ita,pch=".",asp=1)
#
lines(ita,col="blue")
abline(h=seq(37,47,1),v=seq(7,18,1),lty=3,col="grey")
P <- mean(ita); P; points(P[1],P[2],pch=20,col="red")
# Posso tracciare figure di ogni tipo.
# Questo era un esempio di tracciamento per punti.
#
# Un fascio di curve
x1 <- -5; x2 <- 5; y1 <- -10; y2 <- 10
f <- function(x) a*x^2+2*x-2
plot(c(x1,x2),c(y1,y2),type="n",xlab="", ylab="")
abline(h=seq(y1,y2,2), v=seq(x1,x2,1), lty=3, col="grey")
axis(1,pos=0,label=FALSE,col="brown"); axis(2,pos=0,label=FALSE,col="brown")
for (i in 1:10) {a <- -3+i; print(a); curve(f, add=TRUE,col=i)}
# L'aggiunta, sul piano cartesiano, di qualche frase
text(3,2,"retta"); text(-1/4,1,"O")
# I colori sono stati dati non con "nomi", come fatto altre volte,
# ma sotto forma di numeri, che qui richiamo:
dev.new(); plot(c(0,21),c(-1,1),type="n",xlab="", ylab="")
for (c in 1:20) points(c,0,col=c,pch=19)
#
# Su questi aspetti ritorneremo nella parte iniziale del prossimo incontro.
# Purtroppo non ce l'ho fatta a elaborare i dati proposti nei messaggi inviatimi
# da chi segue a distanza. Lo farò per il prossimo incontro ...
#
# Torniamo indietro (<- clicca).