library(rmf)
data(studenti)

# distribuzioni di frequenza e grafici
frequenze(df$esami,cumul=TRUE)
bastoncini(df$esami) # Figura 3.1
frequenze(df$altezza,cumul=TRUE)
frq(df,cumul=TRUE)
pippo<-cbind(df$peso,df$altezza)
frq(pippo,cumul=TRUE)
frq(df$altezza,cumul=TRUE)     # frq e frequenze producono stessi risultati su una variabile
stem(df$altezza)
stem(df$altezza, scale=2)

# Le classi
x <- cutpoints(df$eta,c(19,20,22,29,49))
summary(x)
x <- cutpoints(df$eta,c(19,20,22,29,49), destra=FALSE)
summary(x)

# L'istogramma
hist(df$altezza) # Figura 3.2
istogramma(df$altezza) # Figura 3.3
cumulata(df$altezza)
set.seed(123456)
x <- c(runif(10),sample(c(1:6),20,replace=TRUE))
hist(x,breaks=c(0,1,6)) # Figura 3.4

# EDA
fps(df$altezza)
median(df$crediti,na.rm=TRUE)
mediana(df$crediti)
quantile(df$peso, probs=c(0.05,0.15,1/3,2/3,0.7,0.85,0.95),na.rm=TRUE)
quantile(df$altezza, probs = seq(0, 1, 0.25), na.rm = TRUE)
quartili(df$altezza)
quintili(df$peso)
decili(df$altezza)
boxplot(df$altezza,horizontal = TRUE,range=0) # Figura 3.5
boxplot(df$altezza,horizontal = TRUE)         # Figura 3.6

# Figura 3.7
# Il paradigma della simmetria e le sue violazioni
x <- rbeta(100000,2,5)
y <- rbeta(100000,5,5)
z <- rbeta(100000,5,2)
par(mfrow=c(2,3))
hist(x)
hist(y)
hist(z)
boxplot(x,range=0,horizontal=TRUE)
boxplot(y,range=0,horizontal=TRUE)
boxplot(z,range=0,horizontal=TRUE)
par(mfrow=c(1,1))

# EDA: analisi stratificate
boxplot(df$peso ~ df$sex) # Figura 3.8
boxplot(df$peso ~ df$altezza)
hcat <- cutpoints(df$altezza,seq(150,200,10))
boxplot(df$peso ~ hcat)

curve(dnorm,-3.5,3.5); abline(h=0); box() # Figura 3.9

# Il Q-Q plot
qqn(df$bmi) # Figura 3.10
fps(df$bmi)


# Figura 3.12
# Diagrammi temporali
par(mfrow=c(2,2))

# Tabella 1.3 e figura 1.8 - dal libro di Moore (2005)
privati <- c(7851,7870,7572,7255,7272,7664,7652,7665,7374,7411,7758,
             8389,8882,9324,9984,10502,10799,11723,12110,12380,12601,
             13012,13362,13830,14035,14514,15128,15881,16289,16456,17123)
anni <- c(1971:2001)
cbind(anni, privati)
plot(anni,privati,ylab="Media delle tasse pagate",
                  xlab="Anno accademico",
                  main="Tasse pagate nei college privati")
lines(anni,privati)

# Numero indice elementare (base 1971)
percento <- privati*100/privati[1]
plot(anni,percento,ylab="Numero indice elementare (base 1971)",
                   xlab="Anno accademico",
                   main="Tasse pagate nei college privati")
lines(anni,percento)

# Differenza prima
differenza <- diff(privati) # N.B. c'è un valore in meno
anni <- c(1972:2001)
plot(anni,differenza,ylab="Differenza prima nelle tasse pagate",
                     xlab="Anno accademico",
                     main="Tasse pagate nei college privati")
lines(anni,differenza)

# Differenza relativa o tasso di variazione [x(t)-x(t-1)]/x(t-1)
differenza <- diff(privati)
diffperc <- differenza/privati[1:30]*100 # devo togliere l'ultimo
plot(anni,diffperc,ylab="Tasso di variazione",
                   xlab="Anno accademico",
                   main="Tasse pagate nei college privati")
lines(anni,diffperc)

par(mfrow=c(1,1))


# La media aritmetica
x <- c(1,1,8,6)
mean(x)
x - mean(x)
cbind(x,(x-mean(x)))
sum(x-mean(x))
mean(df$eta)
mean(df$eta,na.rm=TRUE)
media(df$eta)
mean(df$peso,na.rm=TRUE)
by(df$peso,df$sex,mean,na.rm=TRUE)
by(cbind(df$peso,df$altezza,df$bmi),df$sex,mean,na.rm=TRUE)


# Figura 3.13
# La proprietà di minimo della media
set.seed(123456)
par(mfrow=c(2,2))
x <- runif(10)
x
(m <- mean(x))
var(x)*(9/10)

# primo passo
xx <- seq(0,1,0.1)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum((x - a)^2)
}

plot(xx,yy,,xlab='c',ylab='somma degli scarti al quadrato')
abline(v=m,lty=3)

# zoom 1
xx <- seq(0.4,0.6,0.01)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum((x - a)^2)
}

plot(xx,yy,,xlab='c',ylab='somma degli scarti al quadrato')
abline(v=m,lty=3)

# zoom 2
xx <- seq(0.4,0.5,0.001)
yy <- numeric(length(xx))

k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum((x - a)^2)
}

plot(xx,yy,xlab='c',ylab='somma degli scarti al quadrato')
abline(v=m,lty=3)

# zoom finale (centrato sulla media)

xx <- seq((m-0.001),(m+0.001),0.0001)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum((x - a)^2)
}

plot(xx,yy,,xlab='c',ylab='somma degli scarti al quadrato')
abline(v=m,lty=3)

# Media di una distribuzione di frequenze
f <- c(2,1,2,1,1,3,0,1,1,1,0,1,1,0,1,2,0,1,0,1,1,1,0,2)
table(f)
mean(f)
x <- c(0,1,2,3)
n <- c(6,13,4,1)
x*n
cbind(x,n,x*n)
sum(x*n)/sum(n)

# Media ponderata: CFU
wt <- c(8,12,8)
x <- c(27,30,28)
round(mean(x),1)
round(weighted.mean(x, wt),1)
wt <- c(6,13,4,1)
x <- c(0,1,2,3)
weighted.mean(x, wt)

# La media è la stessa ma le distribuzioni sono molto diverse
x <- rep(1:5,40)
y <- rep(1:5,c(10,16,148,16,10))
par(mfrow=c(1,2))
bastoncini(x)
bastoncini(y)
par(mfrow=c(1,1))

# La mediana
x <- c(20,23,22,19,27,30,25)
x
sort(x)
median(x)
round(mean(x),2)
x <- c(x,29)
x
sort(x)
median(x)
x <- c(3,4,7,2,3,1,8,12,1,3,5,6,9)
median(x)
mean(x)
y <- c(x,1000)
median(y)
mean(y)
round(mean(df$bmi,na.rm=T),2)


# Figura 3.14
# La proprietà di minimo della mediana
x <- c(20,23,22,19,27,30,25)
(m <- median(x))
par(mfrow=c(2,2))

# primo passo
xx <- seq(18,30,0.1)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom 1
xx <- seq(20,26,0.01)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom 2
xx <- seq(21,25,0.001)
yy <- numeric(length(xx))

k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom finale (centrato sulla mediana)

xx <- seq((m-0.001),(m+0.001),0.0001)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)
par(mfrow=c(1,1))


# Figura 3.15
# La proprietà di minimo della mediana (N pari)
x <- c(20,23,22,19,27,30,25,29)
(m <- median(x))
par(mfrow=c(2,2))

# primo passo
xx <- seq(18,30,0.1)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom 1
xx <- seq(21,27,0.01)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom 2
xx <- seq(22,26,0.001)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)

# zoom finale (centrato sulla mediana)

xx <- seq((m-1.1),(m+1.1),0.0001)
yy <- numeric(length(xx))
k <- 0
for (a in xx)
{
    k <- k+1
    yy[k] <- sum(abs(x - a))
}
plot(xx,yy,type='l',xlab='c',ylab='somma degli scarti in valore assoluto')
abline(v=m,lty=3)
par(mfrow=c(1,1))



# Misure di variabilità
x <- c(2,4,6)
y <- c(4,4,4)
mean(x);mean(y)
bastoncini(x)
bastoncini(y)
wt  <- c(40,40,40,40,40)
wt1 <- c(10,16,148,16,10)
x   <- c(1,2,3,4,5)
weighted.mean(x, wt);weighted.mean(x, wt1)

# Il Range risente dei dati anomali
x <- c(1,2,3,4,5,6,6,6,6,6,6,7,7)
y <- c(1,2,3,3,4,4,4,4,4,5,5,6,7)
z <- c(1,5,6,6,6,6,6,7,7,7,7,7,7)
max(x)-min(x)
max(y)-min(y)
max(z)-min(z)

# Figura 3.17
par(mfrow=c(1,3))
bastoncini(x)
bastoncini(y)
bastoncini(z)
par(mfrow=c(1,1))


# ...e allora calcoliamo IQR (range interquartilico)
quantile(x, probs = seq(0, 1, 0.25))
quantile(y, probs = seq(0, 1, 0.25))
quantile(z, probs = seq(0, 1, 0.25))
# oppure, con rmf, usiamo la funzione quartili
quartili(x)
quartili(y)
quartili(z)

# La varianza
x <- c(228,226,194,125,180,215,219,218,241,300,292,336,296,319,324,212,241,329,373,177)
mean(x)
x-mean(x)                  # gli scarti dalla media
sum(x-mean(x))             # sommano a zero
(x-mean(x))^2              # scarti dalla media al quadrato
cbind(x,(x-mean(x)),(x-mean(x))^2)
d  <- sum((x-mean(x))^2)   # la devianza
(v <- d/length(x))         # la varianza
sqrt(v)                    # lo scarto quadratico medio
sd(x)                      # la deviazione standard di R
ds(x)                      # la deviazione standard di rmf
var(x)                     # la varianza di R
varianza(x)                # la varianza di rmf
ds(x)*sqrt((length(x)-1)/length(x))     # una riconciliazione
var(x)*(length(x)-1)/length(x)

# Implementiamo le due formule di calcolo della varianza
x <- c(6,2,7,17)
(sum((x-mean(x))^2)/length(x))
(sum(x^2)/length(x)-mean(x)^2)

# Il coefficiente di variazione
(m <- by(df$altezza,df$sex,mean,na.rm=TRUE))
(s <- by(df$altezza,df$sex,sd,na.rm=TRUE))
(cv_F <- s[1]/m[1])
(cv_M <- s[2]/m[2])