Introduzione a R - Università di Bologna

Introduzione a R
Cinzia Viroli
Dipartimento di Scienze Statistiche
Università di Bologna
[email protected]
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
1
0. Premessa
R è uno dei software statistici più utilizzati dalla comunità scientifica. La sua grande diffusione nasce dal
fatto che si tratta di un ambiente di tipo open source ovvero disponibile gratuitamente sotto i vincoli della
GPL (General Public Licence) per diversi sistemi operativi. Si tratta di un ambiente molto versatile e
potente che permette di gestire basi di dati, di effettuare analisi statistiche e di produrre grafici di alto
livello. Inoltre esso rappresenta un vero e proprio linguaggio di programmazione orientato ad oggetti che
consente l'uso di strutture condizionali e cicliche, nonché di funzioni create dall'utente. Per informazioni
sull'installazione del programma:
http://cran.r-project.org/
Queste semplici note rappresentano uno schema molto sintetico e non esaustivo degli argomenti trattati a
lezione. Per maggiori dettagli sull'uso di R si faccia riferimento al volume:
A first course in Statistical programming with R, Braun and Murdoch, Cambridge University Press, 2008
(in inglese)
oppure ai tanti manuali gratuiti scaricabili dal sito di R http://cran.r-project.org sotto la voce Contributed
Documentation.
1. Introduzione al linguaggio R
R è un linguaggio di programmazione ad oggetti. I possibili oggetti sono vettori, matrici, liste, basi di dati
(dataframe), funzioni ecc. Sia in windows che in linux non esiste una interfaccia grafica al programma.
Quando si apre R per la prima volta appare una semplice finestra bianca (chiamata console) con un
cursore che indica la posizione in cui inserire un qualche comando. L'area di lavoro iniziale di R
(chiamata workspace) è vuota. Per verificarlo si digiti:
> ls()
oppure
> objects()
Se si costruisce un semplice oggetto di nome pluto a cui si assegna il valore 3:
> pluto=3
l'area di lavoro ora conterrà il nuovo oggetto.
Si immagini ora di lavorare in R e alla fine della sessione di lavoro si vogliano salvare gli oggetti
costruiti. Occorre quindi salvare il workspace in una determinata directory di lavoro. Alcuni comandi
utili:
> getwd() per conoscere la directory di lavoro
> setwd(“directory”) per impostare la directory di lavoro
> save.image(“nome.RData”) per salvare il workspace
> load(“nome.RData”) per caricare il workspace
> savehistory(“nome”) e loadhistory(“nome”) per salvare
l'elenco dei comandi digitati nella console.
e
caricare
Per uscire da R: q()
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
2
1.1 Lavorare con gli oggetti
> x=2+(3-4*5)/2 #costruisci un oggetto
> x #guardalo
[1] -6.5
> x-2 #...utilizzalo
[1] -8.5
> x=5 #...sovrascrivilo
> y=2.5 #crea un altro oggetto
> x/y #fai qualsiasi operazione con gli oggetti
[1] 2
Il simbolo # è il commento e viene ignorato dall’ambiente.
1.1 Operatori logici
> 2*2==4
[1] TRUE
> 2*3==4
[1] FALSE
> 2*2>=4
[1] TRUE
> 2*3>4
[1] TRUE
> 2*3>=4
[1] TRUE
> 2*3!=4
[1] TRUE
1.2 Vettori
> x=c(2,5,9.5,-3) #costruisci un vettore
> x[2] #seleziona il suo secondo elemento
[1] 5
> x[c(2,4)] #seleziona i suoi elementi nella posizione 2 e 4
[1] 5 -3
> x[x>0] #seleziona i sui elementi positivi
[1] 2.0 5.0 9.5
> x[x>0]-1
[1] 1.0 4.0 8.5
> x[x>0][2]
[1] 5
> x=1:10
> x
[1] 1 2 3 4 5 6 7 8 9 10
> x1=seq(1,1000, length=10) #vettore da 1 a 1000 di ampiezza 10
> x1
[1] 1 112 223 334 445 556 667 778 889 1000
> x2=rep(2,times=10) #ripeti 2 10 volte
> x2
[1] 2 2 2 2 2 2 2 2 2 2
> rep(c(1,3),times=4) #ripeti (1,3) 4 volte
[1] 1 3 1 3 1 3 1 3
> rep(c(1,9),c(3,1)) #ripeti (1,9) 3 e 1 volta rispettivamente
[1] 1 1 1 9
> length(c(x,x1,x2,3))
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
3
[1] 31
1.3 Matrici
> x=matrix(1:10,ncol=5) #costruisci una matrice
> x
[,1] [,2] [,3] [,4] [,5]
[1,] 1 3 5 7 9
[2,] 2 4 6 8 10
> xx=t(x) #trasponila
> x[,1] #seleziona la prima colonna
[1] 1 2
> x[2,] #seleziona la seconda riga
[1] 2 4 6 8 10
> x[3,2] #seleziona l’elemento [3,2]
Error: subscript out of bounds
> x[2,3] #...e quello [2,3]
[1] 6
> x[,4:5] #seleziona solo le colonne 4 e 5
[,1] [,2]
[1,] 7 9
[2,] 8 10
> x[,-c(2,4)] #seleziona le colonne 1, 3 e 5
[,1] [,2] [,3]
[1,] 1 5 9
[2,] 2 6 10
> cbind(1:2,c(1,-2),c(0,9)) #dispone per colonne
[,1] [,2] [,3]
[1,] 1 1 0
[2,] 2 -2 9
> rbind(1:4,c(0,5,-4,6)) #dispone per righe
[,1] [,2] [,3] [,4]
[1,] 1 2 3 4
[2,] 0 5 -4 6
> diag(x[,4:5]) #estrae la diagonale principale
[1] 7 10
> D=diag(1:3) #costruisce una matrice diagonale
> A=matrix(0,3,3) # crea una matrice di zeri
> A=edit(A) # inserisci dei valori nella matrice
1.4 Dataframe
> X=data.frame(a=1:4, sesso=c("M","F","F","M")) #crea un dataframe
#con una variabile ‘quantitativa’ e una ‘qualitativa’.
> X
a sesso
1 1 M
2 2 F
3 3 F
4 4 M
> dim(X)
[1] 4 2
> X$eta=c(2.5,3,5,6.2) #aggiungi una variabile
> X
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
4
a
1
2
3
4
sesso eta
1 M 2.5
2 F 3.0
3 F 5.0
4 M 6.2
> X[X$sesso=="M","eta"] #guarda eta per i maschi
[1] 2.5 6.2
> X$sesso[X$eta<=3] #guarda sesso per cui eta<=3
1.4 Array
> a=array(1:24,c(3,4,2)) #crea un array
> a #guardalo
1.4 Liste
> lista=list(oggetto1=a,oggetto2=X) #crea una lista
> lista$oggetto2 #guarda un oggetto della lista
Esercizio
Calcolare la media aritmetica e la varianza
1 n
( x i − x )2
s2 =
∑
n − 1 i=1
delle vendite di un certo quotidiano di un edicolante in una settimana:
lunedì
martedì
mercoledì
giovedì
venerdì
sabato
50
84
44
43
40
65
2. Le funzioni
In R esistono tantissime funzioni diverse. Ad esempio
mean()
summary()
sum()
var()
Altre possono essere costruite e caricate. Per fare questo è comodo utilizzare un editor. In linux esistono
diverse possibilità fra cui KATE e Emacs Speack Statistics (ESS).
Proviamo a costruire la funzione che calcola la varianza.
varianza=function(vettore)
{
numobs=length(vettore)
varianza=sum((vettore-mean(vettore))^2)
varianza=varianza/(numobs-1)
return(varianza)
}
Una volta che la funzione è stata scritta su un editor per compilarla e quindi renderla disponibile in R si
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
5
può procedere in due modi: o con copia ed incolla, oppure occorre salvarla in un file chiamato ad esempio
“varianza.R” e digitare nella console:
source(“varianza.R”)
Si può costruire anche una funzione che riporti in uscita più oggetti attraverso le liste. Vogliamo costruire
una funzione che restituisca in output media e varianza:
media.var=function(vettore)
{
numobs=length(vettore)
media=mean(vettore)
varianza=sum((vettore-media)^2)
varianza=varianza/(numobs-1)
return(list(media=media,var=varianza))
}
Oppure si può costruire una funzione che prenda in entrata più oggetti. Ad esempio dato un certo capitale
e un certo tasso di interesse mensile si vuole calcolare il montante dopo 6 mesi (capitalizzazione
semplice) M=C(1+6i):
calcola.montante=function(capitale,interesse,mesi=6)
{
montante=capitale*(1+interesse*mesi)
return(montante)
}
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
6
3. La programmazione in R
Il ciclo for:
for (name in vector) {commands}
Esempio: i numeri di Fibonacci
Calcolare i primi 12 numeri di Fibonacci.
Fibonacci=rep(0,12)
Fibonacci[1]=1
Fibonacci[2]=1
for (i in 3:12) Fibonacci[i]=Fibonacci[i-2]+Fibonacci[i-1]
La condizione if:
if (condition) {commands when TRUE}
if (condition) {commands when TRUE} else {commands when FALSE}
L'istruzione while:
while (condition) {statements}
L'istruzione repeat:
repeat {statements}
if (condition) break
Esempio: il metodo di Newton per trovare le radici di una funzione
Si tratta di un algoritmo per determinare le radici di una funzione
f(x)=0
sulla base dell'approssimazione di Taylor. Se f'(x) è la derivata prima di f(x) allora l'algoritmo consiste
nell'iterare
(x0 scelta iniziale)
xn=xn-1
–
f(xn-1)/f'(xn-1)
sino a convergenza.
Esempio
Si applichi l’algoritmo di Newton per determinare le radici della funzione f ( x) = x 3 + 2 x 2 − 7
### Newton con while###
newton.1=function(x0,tolerance=0.000001)
{
x=x0
f=x^3+2*x^2-7
while (abs(f)>tolerance) {
f.prime=3*x^2+4*x
x=x-(f/f.prime)
f=x^3+2*x^2-7
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
7
}
return(x)
}
### Newton con repeat###
newton.2=function(x0,tolerance=0.000001)
{
x=x0
repeat {
f=x^3+2*x^2-7
if (abs(f) < tolerance) break
f.prime=3*x^2+4*x
x=x-(f/f.prime)
}
return(x)
}
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
8
4. Simulazione
4.1 Generazione di numeri casuali
Generazione di numeri pseudo-casuali
Esistono diversi modi per generare una sequenza di numeri nell’intervallo [0,1] che ha tutte le
caratteristiche di una sequenza casuale (ma non lo è perciò la si definisce pseudo-casuale). Un algoritmo
molto semplice è il seguente. Si fissi un numero intero m e un numero intero b<m, non uno fra i divisori
di m. Si scelga un valore di partenza indicato con x0 chiamato seed. Per generare una sequenza di numeri
pseudo-casuali (u1,u2,…,un) si può procedere come segue:
x1 = b x0 mod m
u1 = x1/m
dove u1 è il primo numero pseudo-casuale. Il secondo lo si ottiene come:
x2 = b x1 mod m
u2 = x2/m
e così via.
Esempio
Costruire il programma che generi una sequenza di numeri casuali dati un seme (seed) e i valori per b ed
m. Eseguire il programma con m=7 e b=3.
genera.random=function(n,seme=1,m=100000,b=trunc(sqrt(m)))
{
num.casuali=rep(0,n)
for (i in 1:n) {seme=(b*seme)%%m
num.casuali[i]=seme/m
}
return(num.casuali)
}
Cosa succede se b è un divisore di m? Quale problema si può riscontrare quando m è un numero non
elevato?
Generazione di numeri (pseudo-)casuali da una uniforme
runif(n,min=a,max=b)
Esempio
Generare 1000 osservazioni da due variabili aleatorie X1 e X2 distribuite secondo una uniforme in [0,1].
Farne uno scatterplot. Sono indipendenti? Costruirne poi la variabile somma e la variabile differenza e
farne uno scatterplot. Misurare la correlazione fra le due coppie di variabili.
set.seed(1) ## per fissare un punto di partenza comune
x1=runif(1000)
x2=runif(1000)
plot(x1~x2)
S=x1+x2
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
9
D=x1-x2
plot(S~D)
cor(x1,x2)
cor(S,D)
Generazione di numeri casuali da una bernoulli
Una variabile aleatoria di Bernoulli ha due esiti possibili con probabilità rispettivamente p e (1-p).
Esempio
Una prova di esame consiste di 20 domande alle quali la mia preparazione mi consente di rispondere
correttamente con una probabilità del 20%. Generare uno fra i possibili risultati all'esame.
risposte=runif(20)
risposte.corrette=(risposte<0.2)
table(risposte.corrette)
pie(table(risposte.corrette),col=c("red","green"))
La variabile aleatoria binomiale
La somma X di n prove indipendenti di Bernoulli ciascuna con probabilità p è la variabile binomiale. X
rappresenta il numero di successi in n prove.
dbinom(x,size,prob) # calcola la probabilità di avere “x” successi in
# “size” prove ciascuna con probabilità “prob”
pbinom(x,size,prob) # calcola la probabilità di avere al massimo “x”
# successi ovvero P(X≤x)
rbinom(n,size,prob) # genera n variabili binomiali
Esempio
Quale è la probabilità che in 100 lanci di una moneta non truccata si ottengano esattamente 50 teste e 50
croci? Sorprendentemente questa probabilità è alquanto bassa:
dbinom(50,100,0.5)
[1] 0.07958924
E in 20 lanci quale è la probabilità di ottenere 10 teste e 10 croci? Costruire il grafico che associa ad ogni
possibile numero di teste in 20 lanci la sua probabilità.
x=0:20
prob=rep(0,21)
for (i in 1:21) prob[i]=dbinom(x[i],20,0.5)
barplot(prob,names.arg=x)
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
10
Esempio
Si sa che in un processo industriale il 10% delle tubi prodotte da una vecchia macchina risultano difettose.
Si sa inoltre che vengono prodotti 15 tubi ogni ora. Si definisce un processo fuori controllo quando più di
4 tubi escono difettosi ogni ora. Simulare il numero di pezzi difettosi prodotti ogni ora nell'arco di 24 ore
e contare i processi fuori controllo.
pezzi.dif=rbinom(24,15,0.1)
sum(pezzi.dif>4)
La variabile aleatoria normale
La variabile aleatoria normale è caratterizzata da due parametri: la media e la standard deviation.
dnorm(x,mean,sd) # calcola la densità di probabilità
pnorm(x,mean,sd) # calcola P(X≤x) chiamata funzione di ripartizione
rnorm(n,mean,sd) # genera n variabili normali
qnorm(p,mean,sd) # restituisce il valore di ascissa (quantile) in
# corrispondenza del valore dell'integrale p
Esempio
Simulare 10000 realizzazioni da una distribuzione normale con media 50 e deviazione standard 5.2 e
costruirne l'istogramma.
x=rnorm(10000,50,5.2)
hist(x,probability=TRUE)
lines(density(x))
Esempio
Calcolare la probabilità che le realizzazioni di una variabile aleatoria 'standardizzata' (ovvero con media
zero e standard deviation 1) siano comprese nell'intervallo [-1.946, 1.946].
## Primo modo ##
prob=pnorm(-1.946)
1-2*prob ## in virtù della simmetria della distribuzione
## Secondo modo ##
n=10000000
x=rnorm(n)
compresi=(x>-1.946 & x<1.946)
sum(compresi)/n
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
11
4.2 Integrazione di Monte Carlo
Le realizzazioni dalle varie leggi distributive possono essere utilizzate per approssimare gli integrali del
tipo:
∫
b
a
g (u )du
Siano u1,u2,…,un n realizzazioni indipendenti da variabili aleatorie uniformi nell’intervallo [a,b]. Allora
per la definizione di valore atteso di una variabile aleatoria si sa
E [g (u )] = ∫ g (u )
b
a
1
du
b−a
da cui, per la legge dei grandi numeri
∫
b
a
1 n 
g (u )du = E [g (u )](b − a ) ≅  ∑ u i (b − a )
 n i =1 
Esempio
Approssimare l’integrale
5
∫ sin( x)dx (vero valore –0.7)
2
u=runif(10000,2,5)
mean(sin(u))*(5-2)
Esempio
Approssimare l’integrale doppio
10
7
∫ ∫ sin( x − y)dxdy
3
1
u=runif(10000,1,7)
v=runif(10000,3,10)
mean(sin(u-v))*(10-3)*(7-1)
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
12
5. Alcuni elementi di algebra lineare
Prendiamo una matrice quadrata A:
A=matrix(c(1,2,1.5,0.5),2,2)
Alcune proprietà:
det(A)
# calcola il determinante di A
dim(A)
# restituisce le dimensioni della matrice
diag(A) # restituisce gli elementi sulla diagonale principale di A
sum(diag(A)) # calcola la traccia di A
t(A)
# restituisce la trasposta di A
diag(diag(A)) # restituisce una matrice diagonale con gli elementi
della diagonale principale di A
eigen(A)
# calcola autovalori e autovettori di A
solve(A)
# calcola l’inversa della matrice
Prendiamo 2 matrici con le stesse dimensioni:
A=matrix(c(1,2,1.5,0.5),2,2)
B=matrix(c(0.3,1,1,-0.5),2,2)
A+B
A*B
A%*%B
# somma gli elementi delle due matrici
# moltiplica a due a due fra gli elementi delle due matrici
# fa il prodotto matriciale
Esempio
La funzione solve(A,b) restituisce la soluzione del sistema lineare Ax=b. Determinare x tale che H3
x=b con b=[1 2 3]T e H3 la matrice di Hilbert di ordine 3:
H3=1/cbind(seq(1,3),seq(2,4),seq(3,5))
b=c(1,2,3)
x=solve(H3,b)
Esempio
Trovare la scomposizione di Choleski e la scomposizione in valori singolari della matrice di Hilbert di
ordine 3.
H3.chol=chol(H3)
t(H3.chol)%*%H3.chol
H3.svd=svd(H3)
H3.svd$u%*%diag(H3.svd$d)%*%t(H3.svd$v)
Cinzia Viroli, Dipartimento di Scienze Statistiche, [email protected]
13