% I dati sono su 3 colonne: la prima e` il giorno di rilevazione,
% la seconda e` la durata dell'eruzione del geyser "Old Faithful"
% la terza e` il tempo che intercorre tra l'eruzione e la successiva
%
% Estraiamo i dati delle colonne 2 e 3
du=geyser(:,2);
int=geyser(:,3);

% Calcoliamo gli indici di posizione per du e int
[mean(du), median(du)]
[mean(int), median(int)]

% Calcoliamo il massimo ed il minimo
[max(du), min(du)]
[max(int), min(int)]

% Creiamo una matrice con i dati du e int in colonna
data=[du,int];

% Oppure con questo metodo alternativo
data=[du';int']';

% Calcoliamo varianza, deviazione standard e differenze interquartili
% Nella prima colonna i dati relativi a du e nella seconda quelli relativi a int
[var(data);var(data,1); std(data);std(data,1);range(data);iqr(data)]

% Calcoliamo i quartili dei dati
prctile(data, [25:25:75])

% Calcoliamo i percentili con passo 10
prctile(data, [10:10:100])

% Immettiamo il numero di rilevazioni nella variabile n
n=length(du);

% Creiamo l'istogramma di du e di
hist(du)

% Creiamo una nuova figura e miglioriamo l'istogramma di du
hist(du, 10*range(du)+1)


% Facciamo lo stesso per int
hist(int)
hist(du,range(int)+1)

% Notiamo che in entrambi i casi sembrano esserci due
% comportamenti separati (esamineremo la questione meglio con lo scatterplot)
% Calcoliamo ora frequenze assolute e relative (semplici e cumulative)
fadu=hist(du,10*range(du)+1);
faint=hist(int,range(int)+1);
frdu=fadu/n
frint=faint/n


Frdu(1)=frdu(1);
for i=2:length(frdu)
Frdu(i)=Frdu(i-1)+frdu(i)
end

Frint(1)=frint(1);
for i=2:length(frint)
Frint(i)=Frint(i-1)+frint(i)
end


% Plottiamo le due frequenze semplici in due grafici differenti
plot(min(du):.1:max(du), frdu)
bar(min(du):.1:max(du), frdu)
plot(min(int):max(int), frint)
bar(min(int):max(int), frint)

% Plottiamo le due frequenze cumulative in due grafici differenti
plot(min(du):.1:max(du), Frdu)
bar(min(du):.1:max(du), Frdu)
plot(min(int):max(int), Frint)
bar(min(du):.1:max(du), Frdu)

% Calcoliamo la matrice di covarianza dei dati
cov(data)

% Calcoliamo i coefficienti di correlazione dei dati
corrcoef(data)

% Per renderci conto della bonta` della nostra osservazione sui
% due gruppi differenti di comportamenti e sulla correlazione
% facciamo lo scatterplot
scatter(du,int)
% o equivalentemente
plot(du,int,'o')

% Calcoliamo la retta di regressione
% coefficienti della retta
c=polyfit(du,int,1);
% Valutiamo i valori di du sulla retta
z=polyval(c, du);
% Confrontiamo i due grafici
plot(du,int,'o',du,z)

% Valutiamo la bonta` dell'approssimazione lineare
% Calcoliamo il vettore dei residui e la sua norma2
r=z-int;
plot(z,r,'o')
norma2=norm(r) %=sqrt(sum(r.^2))

intmedio=mean(int)
DT=norm((int-intmedio))^2 %=sum((int-intmedio).^2)  Devianza Totale
DS=norm(z-intmedio)^2 %=sum((z-intmedio).^2)        Devianza spiegata
DR=norma2^2 %=sum(r.^2)                             Devianza residua
DS+DR-DT  %deve essere 0
R2=DS/DT      % e` buono quanto piu` e` prossimo a 1

% Cerchiamo ora la moda delle durate
mat1=[min(du):.1:max(du);frdu]';
mat2=sortrows(mat1,2);
modadu=mat2(length(frdu),1)

% Cerchiamo ora la moda dei tempi di intereruzione
mat1=[min(int):max(int);frint]';
mat2=sortrows(mat1,2);
modadu=mat2(length(frint),1)


% Dividiamo le osservazioni in due gruppi poiche` dall'esame dei due istogrammi 
% e dallo scatterplot sembrano esserci due comportamenti distinti
% (i dati sullo scatterplot sembrano essere contenuti in due rettangoli 
% opposti)
% Scegliamo di dividerli guardando la durata intereruzione
% calcoleremo quindi gli indici di dispersione e posizione per i due 
% gruppi. Osservata una durata T di una eruzione,
% la previsione del tempo che intercorrera` fino alla prossima
% eruzione sara`, in prima approssimazione, dato dal valore medio 
% della durata dei tempi di intereruzione del gruppo cui T appartiene

% Dividiamo le osservazioni in due gruppi
% du<=3, du>3
sdata=sortrows(data,1);
i=1;
while sdata(i,1)<=3
i=i+1;
end

% Costruiamo le due matrici di osservazioni relative ai due gruppi
data1=sdata(1:(i-1),1:2);
data2=sdata(i:n,1:2);

% Calcoliamo gli indici di posizione e dispersione del gruppo 1
[mean(data1); median(data1)]
[var(data1);var(data1,1)]

% Calcoliamo gli indici di posizione e dispersione del gruppo 2
[mean(data2); median(data2)]
[var(data2);var(data2,1)]

% Boxplot tempi intereruzione
subplot(2,2,1), boxplot(data(:,2))        %subplot(M,N,K) crea M*N spazi
subplot(2,2,2), boxplot(data1(:,2))       %il cui k-esimo e' occupato dal plot
subplot(2,2,3), boxplot(data2(:,2))

% Boxplot durate eruzioni
subplot(1,3,1), boxplot(data(:,1))
subplot(1,3,2), boxplot(data1(:,1))
subplot(1,3,3), boxplot(data2(:,1))