giovedì 3 giugno 2021

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). URL consultato il 02/01/2023.

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

Nessun commento:

Posta un commento