Questo blog ha come obiettivo quello di fornire gli strumenti per "rompere il ghiaccio" con R, lasciando poi agli eventuali interessi e all'iniziativa personale l'approfondimento del linguaggio di programmazione, che ciascuno potrà realizzare in base ai proprio bisogni. Tuttavia occasionalmente può essere interessante affrontare qualche argomento un poco più tecnico, come quello che vediamo oggi.
Il punto di partenza è rappresentato dal calcolo delle statistiche elementari parametriche. Per un campione che include n dati:
a) la media, la misura di posizione dei dati, viene calcolata come
b) la deviazione standard, la misura di dispersione dei dati, viene calcolata come
essendo
Invece di calcolare le statistiche elementari parametriche servendoci di qualcuna delle funzioni disponibili in R le vogliamo calcolare seguendo alla lettera la loro definizione matematica. Come si vede i calcoli sono elementari, e di fatto l'unico problema da risolvere in termini di programmazione è trovare un metodo per il calcolo della sommatoria Σ che ci serve per ottenere:
1) la somma di tutti i dati
da dividere poi per il numero n dei dati al fine di calcolare la media;
2) la somma dei quadrati delle differenze tra ciascun dato e la media, ovvero la devianza
da dividere poi per il numero dei dati meno uno per calcolare la varianza dalla quale otteniamo infine la deviazione standard.
Il primo metodo per il calcolo di una sommatoria che impieghiamo fa ricorso ad una iterazione (o ciclo o loop), una struttura di controllo che consente di eseguire ripetutamente una o più istruzioni, fino al verificarsi di una condizione specificata.
Accertatevi innanzitutto di avere installato il pacchetto psychTools che contiene i dati che ci servono, o in alternativa procedete come indicato nel post Il set di dati galton nel quale trovate anche illustrato il significato dei dati impiegati. Quindi copiate questo script, incollatelo nella Console di R e premete ↵ Invio:
# STATISTICHE ELEMENTARI PARAMETRICHE con loop
#
library(psychTools) # carica il pacchetto che include il set di dati galton
#
library(psychTools) # carica il pacchetto che include il set di dati galton
x <- galton$parent # salva nel vettore x le altezze dei padri
#
n <- length(x) # calcola il numero dei dati
#
n <- length(x) # calcola il numero dei dati
#
### (a)
#
sx <- 0 # inizializza la sommatoria dei dati
for (i in 1:n) {
sx <- sx+x[i] # calcola la sommatoria dei dati
}
sx <- 0 # inizializza la sommatoria dei dati
for (i in 1:n) {
sx <- sx+x[i] # calcola la sommatoria dei dati
}
#
media <- sx/n # calcola la media
media <- sx/n # calcola la media
#
### (b)
### (b)
#
devianza <- 0 # inizializza la sommatoria degli scarti quadratici (devianza)
#
for (i in 1:n) {
devianza <- devianza+(x[i]-media)^2 # calcola la sommatoria degli scarti quadratici
}
#
varianza <- devianza/(n-1) # calcola la varianza
ds <- sqrt(varianza) # calcola la deviazione standard
devianza <- 0 # inizializza la sommatoria degli scarti quadratici (devianza)
#
for (i in 1:n) {
devianza <- devianza+(x[i]-media)^2 # calcola la sommatoria degli scarti quadratici
}
#
varianza <- devianza/(n-1) # calcola la varianza
ds <- sqrt(varianza) # calcola la deviazione standard
#
### (c)
### (c)
#
statistica <- c("numero dei dati =", "media =", "varianza =", "deviazione standard =") # predispone i nomi delle statistiche calcolate
risultato <- c(n, media, varianza, ds) # predispone i risultati delle statistiche calcolate
print(data.frame(statistica, risultato), row.names=FALSE) # mostra la tabella dei risultati
statistica <- c("numero dei dati =", "media =", "varianza =", "deviazione standard =") # predispone i nomi delle statistiche calcolate
risultato <- c(n, media, varianza, ds) # predispone i risultati delle statistiche calcolate
print(data.frame(statistica, risultato), row.names=FALSE) # mostra la tabella dei risultati
#
Dopo avere caricato con la prima riga di codice il pacchetto che include i dati, con la seconda riga di codice al vettore x sono assegnati (<-) i valori della variabile galton$parent che vogliamo utilizzare (si rimanda nuovamente al post Il set di dati galton per la loro storia e il loro significato).
Ora che abbiamo in x i dati che dobbiamo elaborare impieghiamo (terza riga di codice) la funzione length() per calcolarne il numero n.
Il primo loop (a) ci serve per calcolare la sommatoria dei dati che chiamiamo sx e che come prima cosa (quarta riga) inizializziamo ponendola uguale a 0 (sx <- 0).
Per calcolare la sommatoria dei dati impieghiamo la funzione for() che esegue il codice compreso tra la parentesi graffa aperta { e la parentesi graffa chiusa }. Il loop for (i in 1:n) procede facendo variare l'indice i da 1 a n e, per ognuno dei valori assunti da i, somma in sx il valore del dato i-esimo x[i]. Alla prima iterazione verrà sommato in sx il valore del dato x1, alla seconda iterazione verrà aggiunto il valore del dato x2, alla terza iterazione verrà aggiunto il valore del dato x3, e così via fino al dato xn, dopo il quale il loop si interrompe: a questo punto avremo in sx la sommatoria dei valori degli n dati. Questa prima parte si conclude con il calcolo della media dei dati ottenuta dividendo la sommatoria dei dati per il loro numero (media <- sx/n).
Il secondo loop (b) calcola la somma dei quadrati delle differenze tra ciascun dato e la media che chiamiamo devianza e che viene come prima cosa inizializzata ponendola uguale a 0 (devianza <- 0).
Quindi con un loop identico al precedente, facendo variare l'indice i da 1 a n, per ognuno dei valori assunti da i viene sommata nella variabile devianza la differenza tra ciascun dato e la media dei dati elevata al quadrato cioè (x[i]-media)^2 ottenendo in tal modo nella variabile devianza la sommatoria desiderata, che divisa per n - 1 ci consente di ottenere la varianza, la cui radice quadrata è infine la deviazione standard, come dalle formule riportate all'inizio.
La terza e ultima parte dello script (c) si limita a generare e a riportare nella Console di R una tabella con la sintesi dei risultati ottenuti:
statistica risultato
numero dei dati = 928.000000
media = 68.308190
varianza = 3.194561
deviazione standard = 1.787333
Il secondo metodo per il calcolo di una sommatoria impiega la vettorializzazione. Molte funzioni in R sono vettorializzate, il che significa che le operazioni specificate all'interno della funzione sono effettuate in parallelo su tutti gli elementi del vettore contenente i dati.
Copiate questo script, incollatelo nella Console di R e premete ↵ Invio:
# STATISTICHE ELEMENTARI PARAMETRICHE con vettorializzazione
#
library(psychTools) # carica il pacchetto che include il set di dati galton
x <- galton$parent # salva la variabile nel vettore x
#
n <- length(x) # calcola il numero dei dati
#
### (a)
#
media <- sum(x)/n # calcola la media
#
### (b)
#
varianza <- sum((x-media)^2)/(n-1) # calcola la varianza
ds <- sqrt(varianza) # calcola la deviazione standard
#
### (c)
#
statistica <- c("numero dei dati =", "media =", "varianza =", "deviazione standard =") # predispone i nomi delle statistiche calcolate
risultato <- c(n, media, varianza, ds) # predispone i risultati delle statistiche calcolate
print(data.frame(statistica, risultato), row.names=FALSE) # mostra la tabella dei risultati
#
Si vede chiaramente in che modo R realizza la vettorializzazione nella quinta riga di codice: la funzione sum() calcola la sommatoria dei quadrati delle differenze tra ciascun dato e la media (la devianza), necessaria a sua volta per il calcolo della varianza, impiegando direttamente come argomento il vettore x contenente gli n dati, ragion per cui l'espressione (x-media)^2 si intende applicata, e viene applicata, a tutti gli n dati dati del vettore x. Confrontato con quello dello script precedente il codice della parte (a) e della parte (b) è molto più conciso e più facile da leggere, oltre ad essere più efficiente di quello impiegato da un linguaggio non vettorializzato, mentre la parte (c) dello script è semplicemente identica alla precedente (risultati inclusi).
Se volete infine controllare la correttezza dei calcoli effettuati con loop e con vettorializzazione potete confrontarne i risultati con quelli ottenuti con queste tre righe di codice che fanno direttamente ricorso alle funzioni statistiche standard della versione base di R (che si assumono fornire risultati corretti):
#
mean(galton$parent)
var(galton$parent)
sd(galton$parent)
#
In conclusione abbiamo visto attraverso un breve esempio come con R sia possibile effettuare la ripetizione ciclica di un calcolo, o più in generale di una qualsiasi operazione, sia ricorrendo alla classica iterazione (loop), sia impiegando funzioni che prevedono la vettorializzazione [1]. E questo ci ha consentito di illustrare due tecniche di programmazione semplici ma interessanti, che possono risultare utili nel corso dell'utilizzo di R.
----------
[1] R-bloggers. Vectorization in R: Why? Nell'articolo sono riportati anche i riferimenti ad alcune risorse sull'argomento. URL consultato il 20/09/2021: https://bit.ly/3nOsxuS
