Calcolare la media mobile e la mediana mobile

Data una sequenza di valori numerici la media mobile è la media calcolata in corrispondenza di ciascun dato su un numero fisso di dati, determinato da una finestra di ampiezza prestabilita che scorre i dati e nella quale il dato più vecchio viene ogni volta sostituito con un nuovo dato.

Così ad esempio per una finestra con l'ampiezza di 5 dati, in corrispondenza dei primi 4 dati la media mobile non può essere calcolata, in corrispondenza del dato numero 5 viene calcolata la media dei dati da 1 a 5, in corrispondenza del dato numero 6 viene calcolata la media dei dati da 2 a 6, in corrispondenza del dato numero 7 viene calcolata la media dei dati da 3 a 7 e così via. La mediana mobile a sua volta è la mediana calcolata in corrispondenza di ciascun dato nello stesso modo.

Per calcolare la media mobile e la mediana mobile facciamo ricorso al pacchetto zoo, che deve essere scaricato dal CRAN.

Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# CALCOLARE LA MEDIA MOBILE E LA MEDIANA MOBILE 

#

library(zoo) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediamob <- rollmean(valore_y, k=7, fill=NA, align='right') # calcola la media mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata$medianamob <- rollmedian(valore_y, k=7, fill=NA, align='right') # calcola la mediana mobile (dato corrente + 6 precedenti) e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media mobile e mediana mobile", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone il grafico della media mobile

lines(mydata$rileva_x, mydata$medianamob, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone il grafico della mediana mobile

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media mobile dei valori osservati", "Mediana mobile dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto zoo i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media mobile e la mediana mobile;

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna/vettore mydata$mediamob nella quale è riportata la media mobile calcolata in corrispondenza di ogni dato;

→ la colonna/vettore mydata$medianamob nella quale è riportata la mediana mobile calcolata in corrispondenza di ogni dato.

La media mobile e la mediana mobile sono calcolate rispettivamente con le funzioni rollmean() e rollmedian() che impiegano i seguenti argomenti:

→ la variabile (valore_y) sulla quale calcolare la media mobile;

→ l'ampiezza della finestra (k=7) che specifica che la media (o mediana) mobile va calcolata, in corrispondenza di ogni dato, su quel dato e sui sei dati che lo precedono;

→ fill=NA che specifica di riportare NA in corrispondenza dei valori non computabili;

→ align='right' che - in quanto con k=7 la media e la mediana possono essere calcolate solamente a partire dal settimo dato - specifica di riportare i valori non computabili NA nelle prime sei posizioni.

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto zoo [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media mobilea (mediamob) e mediana mobile (medianamob)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y mediamob medianamob

1         1        5       NA         NA

2         2        9       NA         NA

3         3       12       NA         NA

4         4        7       NA         NA

5         5        8       NA         NA

6         6       13       NA         NA

7         7       26 11.42857          9

8         8        7 11.71429          9

9         9        8 11.57143          8

10       10       19 12.57143          8

11       11       10 13.00000         10

12       12        7 12.85714         10

13       13        5 11.71429          8

14       14       36 13.14286          8

15       15       12 13.85714         10

16       16       18 15.28571         12

17       17       31 17.00000         12

18       18       17 18.00000         17

19       19       16 19.28571         17

20       20       43 24.71429         18

21       21        7 20.57143         17

22       22       12 20.57143         17

23       23        9 19.28571         16

24       24       27 18.71429         16

25       25       15 18.42857         15

26       26       18 18.71429         15

27       27       24 16.00000         15

28       28       29 19.14286         18

29       29       32 22.00000         24

30       30       23 24.00000         24

31       31       29 24.28571         24

32       32       12 23.85714         24

33       33        7 22.28571         24

34       34       48 25.71429         29

35       35       34 26.42857         29

36       36       45 28.28571         29

37       37       19 27.71429         29

38       38       12 25.28571         19

39       39       16 25.85714         19

40       40       18 27.42857         19

41       41       22 23.71429         19

42       42       31 23.28571         19

43       43       38 22.28571         19

44       44        7 20.57143         18

45       45       11 20.42857         18

46       46        4 18.71429         18

47       47       38 21.57143         22

48       48       36 23.57143         31

49       49       43 25.28571         36

50       50       15 22.00000         15

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media mobile con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana mobile con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplicativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto zoo: S3 Infrastructure for Regular and Irregular Time Series (Z's Ordered Observations).

https://cran.r-project.org/web/packages/zoo/index.html

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.

Calcolare la media cumulativa e la mediana cumulativa

Data una sequenza di valori numerici la media cumulativa è la media, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Così ad esempio in corrispondenza del dato numero 5 la media cumulativa è la media dei dati da 1 a 5, in corrispondenza del dato numero 6 la media cumulativa è la media dei dati da 1 a 6, in corrispondenza del dato numero 7 la media cumulativa è la media dei dati da 1 a 7 e così via. Lo stesso principio vale per la mediana cumulativa che a sua volta è la mediana, calcolata in corrispondenza di ciascun dato, di tutti i dati dal primo fino a quello corrente.

Per il calcolo della media cumulativa e della mediana cumulativa impieghiamo le funzioni contenute nel pacchetto cumstats, che deve essere scaricato dal CRAN.

Copiate e incollate questo script nella Console di R e premete ↵ Invio:

# CALCOLARE LA MEDIA CUMULATIVA E LA MEDIANA CUMULATIVA 

#

library(cumstats) # carica il pacchetto

rileva_x <- seq(1:50) # numero progressivo della rilevazione in ascisse

valore_y <- c(5, 9, 12, 7, 8, 13, 26, 7, 8, 19, 10, 7, 5, 36, 12, 18, 31, 17, 16, 43, 7, 12, 9, 27, 15, 18, 24, 29, 32, 23, 29, 12, 7, 48, 34, 45, 19, 12, 16, 18, 22, 31, 38, 7, 11, 4, 38, 36 ,43, 15) # valore osservato corrispondente in ordinate

mydata <- data.frame(rileva_x, valore_y) # combina i vettori in una matrice

#

mydata$mediacum <- cummean(mydata$valore_y) # calcola la media cumulativa e la aggiunge in una nuova colonna

mydata$medianacum <- cummedian(mydata$valore_y) # calcola la mediana cumulativa e la aggiunge in una nuova colonna

mydata # mostra la matrice con i dati definitivi

#

plot(mydata$rileva_x, mydata$valore_y, xlim=c(0,50), ylim=c(0,60), type="l", lty=1, lwd=1, col="black", main="Valori osservati, media cumulativa e mediana cumulativa", xlab="Rilevazioni", ylab="Valori osservati") # grafico dei dati originali

lines(mydata$rileva_x, mydata$mediacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="red") # sovrappone grafico della media cumulativa

lines(mydata$rileva_x, mydata$medianacum, xlim=c(0,50), ylim=c(0,60), type="l", lty=2, lwd=1, col="blue") # sovrappone grafico della mediana cumulativa

#

legend(0,60, legend=c("Valori osservati nelle rilevazioni", "Media cumulativa dei valori osservati", "Mediana cumulativa dei valori osservati"), col=c("black", "red", "blue"), lty=c(1,2,2), lwd=c(1,1,1), cex=0.8) # aggiunge la legenda

#

Dopo avere caricato il pacchetto cumstats i dati impiegati come esempio sono generati in questo modo:

→ nel vettore rileva_x viene riportato il numero progressivo della rilevazione;

→ nel vettore valore_y viene riportato il valore numerico osservato in corrispondenza di ciascuna rilevazione, valore numerico sul quale calcoleremo la media cumulativa e la mediana cumulativa.

→ con la funzione dataframe() i due vettori sono combinati nella matrice mydata [1].

A questo punto alla matrice mydata sono aggiunte due colonne:

→ la colonna (vettore) mediacum nella quale è riportata la media cumulativa calcolata in corrispondenza di ogni dato con la funzione cummmean();

→ la colonna (vettore) medianacum nella quale è riportata la mediana cumulativa calcolata in corrispondenza di ogni dato con la funzione cummedian().

Per i dettagli delle funzioni impiegate digitare help(nomedellafunzione) nella Console di R. Al bisogno potete consultare il manuale del pacchetto cumstats [2].

Non resta quindi che mostrare il contenuto della matrice mydata nella versione definitiva contenente il numero della rilevazione (rileva_x), il valore osservato (valore_y) e le rispettive media cumulativa (mediacum) e mediana cumulativa (medianacum)

> mydata # mostra la matrice con i dati definitivi

   rileva_x valore_y  mediacum medianacum

1         1        5  5.000000        5.0

2         2        9  7.000000        7.0

3         3       12  8.666667        9.0

4         4        7  8.250000        8.0

5         5        8  8.200000        8.0

6         6       13  9.000000        8.5

7         7       26 11.428571        9.0

8         8        7 10.875000        8.5

9         9        8 10.555556        8.0

10       10       19 11.400000        8.5

11       11       10 11.272727        9.0

12       12        7 10.916667        8.5

13       13        5 10.461538        8.0

14       14       36 12.285714        8.5

15       15       12 12.266667        9.0

16       16       18 12.625000        9.5

17       17       31 13.705882       10.0

18       18       17 13.888889       11.0

19       19       16 14.000000       12.0

20       20       43 15.450000       12.0

21       21        7 15.047619       12.0

22       22       12 14.909091       12.0

23       23        9 14.652174       12.0

24       24       27 15.166667       12.0

25       25       15 15.160000       12.0

26       26       18 15.269231       12.0

27       27       24 15.592593       12.0

28       28       29 16.071429       12.5

29       29       32 16.620690       13.0

30       30       23 16.833333       14.0

31       31       29 17.225806       15.0

32       32       12 17.062500       14.0

33       33        7 16.757576       13.0

34       34       48 17.676471       14.0

35       35       34 18.142857       15.0

36       36       45 18.888889       15.5

37       37       19 18.891892       16.0

38       38       12 18.710526       15.5

39       39       16 18.641026       16.0

40       40       18 18.625000       16.0

41       41       22 18.707317       16.0

42       42       31 19.000000       16.5

43       43       38 19.441860       17.0

44       44        7 19.159091       16.5

45       45       11 18.977778       16.0

46       46        4 18.652174       16.0

47       47       38 19.063830       16.0

48       48       36 19.416667       16.5

49       49       43 19.897959       17.0

50       50       15 19.800000       16.5

che sono impiegati per la successiva rappresentazione grafica

che viene realizzata con tre sole funzioni:

→ la funzione plot() che traccia il grafico del valore osservato mydata$valore_y in corrispondenza di ciascuna rilevazione mydata$rileva_x impiegando una linea (type="l") continua (lty=1) di larghezza unitaria (lwd=1) e di colore nero (col="black");

→ la funzione lines() che viene impiegata una prima volta per sovrapporre il grafico della media cumulativa con una linea tratteggiata (lty=2) di colore rosso (col="red") e una seconda volta per sovrapporre il grafico della mediana cumulativa con una linea tratteggiata (lty=2) di colore blu (col="blue") 

→ la funzione legend() che riporta la legenda esplocativa dei grafici [3].

Ovviamente lo script può essere riadattato per rappresentare dati differenti ma anche per riportare solamente l'uno o l'altro dei grafici.

----------

[1] Parliamo di array o vettore nel caso di dati numerici monodimensionali, disposti su una sola riga, 

8

6

11

7

di matrice nel caso di dati numerici disposti su più righe e più colonne

8	9	15	14
6	7	18	12
11	8	17	13
7	4	19	17

e di tabella nei casi in cui il contenuto, disposto su più righe e più colonne, è rappresentato oltre che da dati numerici, anche da testo e/o operatori logici

M	7	9	VERO
F	3	12	VERO
F	5	10	FALSO

[2] Vedere la voce Reference manual del pacchetto cumstats: Cumulative Descriptive Statistics.

https://cran.r-project.org/web/packages/cumstats/index.html

[3] Per un approfondimento sull'impiego della funzione legend() rimando al post Aggiungere la legenda a un grafico.

Statistiche elementari parametriche

Siamo al terzo passo delle valutazioni che è sempre necessario effettuare nell'ambito di un percorso logico che prevede, per il calcolo delle statistiche elementari di una singola variabile (analisi univariata):

→ analisi esplorativa dei dati;

→ esecuzione dei test di normalità (gaussianità) per valutare se i dati seguono una distribuzione gaussiana;

→ calcolo delle statistiche elementari parametriche (media, deviazione standard, varianza, quantili parametrici) se i dati seguono una distribuzione gaussiana;

→ calcolo delle statistiche elementari non parametriche (mediana, deviazione assoluta mediana o MAD, quartili e altri quantili non parametrici) se i dati non seguono una distribuzione gaussiana.

Si tratta del passo conclusivo per la variabile altezza (espressa in cm) rilevata in 202 atleti australiani, in quanto con i test di normalità abbiamo stabilito che la variabile è distribuita in modo gaussiano: questa era la sintesi grafica dei risultati.

Ora copiate questo script, incollatelo nella Console di R e premete ↵ Invio:

# STATISTICHE ELEMENTARI PARAMETRICHE 

#

library(DAAG) # carica il pacchetto DAAG che include il set di dati ais

mydata <- as.matrix(ais[c(10)]) # ci interessa la sola colonna con l'altezza

#

# statistiche elementari parametriche 

#

min(mydata) # valore minimo

mean(mydata) # media 

max(mydata) # valore massimo

range(mydata) # range 

max(mydata)-min(mydata) # campo di variazione

var(mydata) # varianza

sd(mydata) # deviazione standard

(sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

#

# confronta media e mediana

#

media <- mean(mydata) # calcola la media

mediana <- median(mydata) # calcola la mediana

data.frame(media, mediana) # le mette a confronto

#

# confronta deviazione standard e MAD

#

ds <- sd(mydata) # calcola la deviazione standard

mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

data.frame(ds, mad) # le mette a confronto

#

# confronta quantili parametrici e non parametrici

#

qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

data.frame(qpar, qnon, dper) # li mette a confronto

#

Non avrebbe senso ripetere l'analisi esplorativa dei dati, i test di normalità e la sintesi grafica dei risultati - per questi si rimanda agli script riportati nel percorso logico indicato all'inizio. Per cui qui ci limitiamo a calcolare le usuali statistiche elementari parametriche, riportando per confronto i risultati delle statistiche elementari non parametriche corrispondenti. 

Dopo avere caricato il pacchetto DAAG [1] con la tabella ais, la variabile che ci interessa, contenuta nella colonna 10 della tabella ais viene assegnata (<-) all'oggetto mydata, sul quale sono calcolate le statistiche elementari parametriche con le funzioni:

→ min() per riportare il valore minimo osservato;

→ mean() per calcolare la media;

→ max() per riportare il valore massimo osservato;

→ range() per calcolare il range, che qui è correttamente;

→ max(mydata)-min(mydata) per calcolare il campo di variazione, talora indicato (impropriamente) come "range";

→ var() per calcolare la varianza,

→ sd() per calcolare la deviazione standard;

→ e con l'espressione (sd(mydata)/mean(mydata))*100 che calcola il coefficiente di variazione percentuale.

> min(mydata) # valore minimo

[1] 148.9

> mean(mydata) # media 

[1] 180.104

> max(mydata) # valore massimo

[1] 209.4

> range(mydata) # range 

[1] 148.9 209.4

> max(mydata)-min(mydata) # campo di variazione

[1] 60.5

> var(mydata) # varianza

         ht

ht 94.76038

> sd(mydata) # deviazione standard

[1] 9.734494

> (sd(mydata)/mean(mydata))*100 # coefficiente di variazione CV %

[1] 5.404931

Sono poi messe a confronto la media con la mediana

> media <- mean(mydata) # calcola la media

> mediana <- median(mydata) # calcola la mediana

> data.frame(media, mediana) # le mette a confronto

    media mediana

1 180.104   179.7

la deviazione standard con la MAD [2]

> ds <- sd(mydata) # calcola la deviazione standard

> mad <- mad(mydata) # calcola la Median Absolute Deviation (about the median) o MAD

> data.frame(ds, mad) # le mette a confronto

        ds     mad

1 9.734494 8.96973

e i quantili parametrici con i quantili non parametrici, riportando nell'ultima colonna la differenza percentuale tra i due

> qpar <- round(qnorm(c(seq(0.025, 0.975, 0.025)), mean=mean(mydata), sd=sd(mydata)), digits=2) # calcola i quantili parametrici

> qnon <- round(quantile(mydata, probs=seq(0.025, 0.975, 0.025)), digits=2) # calcola i quantili non parametrici

> dper <- abs(round(100*(qpar-qnon)/((qpar+qnon)/2), digits=2)) # calcola la differenza in percentuale

> data.frame(qpar, qnon, dper) # li mette a confronto

        qpar   qnon dper

2.5%  161.02 158.98 1.28

5%    164.09 163.96 0.08

7.5%  166.09 167.35 0.76

10%   167.63 169.17 0.91

12.5% 168.91 170.36 0.85

15%   170.01 171.40 0.81

17.5% 171.01 172.22 0.71

20%   171.91 172.76 0.49

22.5% 172.75 173.52 0.44

25%   173.54 174.00 0.26

27.5% 174.29 174.18 0.06

30%   175.00 175.00 0.00

32.5% 175.69 175.73 0.02

35%   176.35 176.07 0.16

37.5% 177.00 177.30 0.17

40%   177.64 177.94 0.17

42.5% 178.26 178.44 0.10

45%   178.88 178.99 0.06

47.5% 179.49 179.60 0.06

50%   180.10 179.70 0.22

52.5% 180.71 180.10 0.34

55%   181.33 180.36 0.54

57.5% 181.94 181.00 0.52

60%   182.57 182.66 0.05

62.5% 183.21 183.06 0.08

65%   183.85 183.90 0.03

67.5% 184.52 184.67 0.08

70%   185.21 185.17 0.02

72.5% 185.92 185.60 0.17

75%   186.67 186.18 0.26

77.5% 187.46 187.28 0.10

80%   188.30 188.62 0.17

82.5% 189.20 189.18 0.01

85%   190.19 190.67 0.25

87.5% 191.30 191.94 0.33

90%   192.58 192.79 0.11

92.5% 194.12 193.86 0.13

95%   196.12 195.17 0.49

97.5% 199.18 197.48 0.86

Le differenze trascurabili tra statistiche elementari parametriche e statistiche elementari non parametriche riflettono il fatto che nel caso di una distribuzione perfettamente gaussiana [3] le due coincidono - è quindi assolutamente lecito impiegare le seconde al posto delle prime. Mentre il contrario è sbagliato - nel caso di una distribuzione non gaussiana le statistiche elementari parametriche non devono essere impiegate.

Lo script può essere riutilizzato molto semplicemente, assegnando (<-) all'oggetto mydata i propri dati.

----------

[1] Vedere il post Il set di dati ais nel quale trovate anche come caricare i dati della tabella senza impiegare il pacchetto DAAG. 

[2] La Median Absolute Deviation about median o MAD ovvero la deviazione assoluta mediana (dalla mediana) è l'equivalente non parametrico della deviazione standard. Vedere: Rousseeuw PJ, Croux C. Alternatives to the Median Absolute Deviation. Journal of the American Statistical Association 88 (424), 1273-1283, 1993.

https://www.jstor.org/stable/2291267?seq=1#page_scan_tab_contents

[3] Un esempio di dati distribuiti in modo perfettamente gaussiano è riportato nel post La distribuzione gaussiana.