Stat101 - parte 1
Valerio Calà
12/10/2020
Statistics & Probability
In questo documento markdown vengono svolti esercizi base di Statistica e Calcolo delle Probabilità. Programma del documento:
Distribuzioni di frequenza, indici di posizione, variabilità, forma e concentrazione
Distribuzioni di frequenza multiple per caratteri qualitativi (indici di connessione)
Distribuzioni di frequenza multiple per caratteri quantitativi (indici di correlazione)
Indici di posizione, variabilità, forma e concentrazione
In questa sezione carichiamo dei dati provenienti da una popolazione \(X\) (carattere univariato, quantitativo trasferibile). L’aggettivo ‘trasferibile’ indica il fatto che la totalità delle modalità di \(X\) è prefissata (somma pari a \(T\)). \[\sum_{i=1}^{n}{X_i}=T\]
# fenomeno X quantitativo; si osserva X=x
x = c(0.5, 2.1, 3.4, 4.1, 4.6, 5.7, 6.2, 6.6, 7.8, 8.1, 8.3,
8.4, 8.6, 8.9, 9.2, 9.8, 10, 10.3, 10.5, 10.6, 10.8, 11.2,
11.3, 11.6, 11.7, 12.1, 12.5, 12.6, 12.8, 12.9, 12.9, 13.2,
14.4, 15, 15.5, 16.3, 17, 17.3, 18.5, 23.5)
n = length(x)
# dividi x in classi d'intervallo d'ampiezza 3
estremi_classi = seq(3, 24, by=3)
intervalli = c("0-3","3-6","6-9","9-12","12-15","15-18","18-21","21-24")
frequenze_cumulate = c()
f_ripartizione = function(estremi, x, h=3) {
for (i in 1:length(estremi)) {
frequenze_cumulate[i] = as.data.frame(x) %>% filter(x<=estremi[i]) %>% count()
}
as.numeric(frequenze_cumulate)
}
frequenze_cumulate = f_ripartizione(estremi=estremi_classi, x=x)
frequenze_marginali = c()
f_marginale = function(fdr) {
frequenze_marginali[1] = fdr[1]
for (i in 2:length(fdr)) {
frequenze_marginali[i] = fdr[i] - fdr[i-1]
}
as.numeric(frequenze_marginali)
}
frequenze_marginali = f_marginale(fdr=frequenze_cumulate)
p_x = frequenze_marginali/n
F_x = frequenze_cumulate/n
m_x = c()
media_estremi = function(estremi) {
m_x[1] = estremi[1]/2
for (i in 2:length(estremi)) {
m_x[i] = m_x[i-1] + 3
}
as.numeric(m_x)
}
m_x = media_estremi(estremi_classi)
crea_distribuzione = function(intervalli, estremi, n, N, p, Effe, m) {
data.frame(intervalli, estremi, n, N, p, Effe, m)
}
DF = crea_distribuzione(intervalli, estremi_classi, n=frequenze_marginali, N=frequenze_cumulate, p_x, F_x, m_x)
DF## intervalli estremi n N p Effe m
## 1 0-3 3 2 2 0.050 0.050 1.5
## 2 3-6 6 4 6 0.100 0.150 4.5
## 3 6-9 9 8 14 0.200 0.350 7.5
## 4 9-12 12 11 25 0.275 0.625 10.5
## 5 12-15 15 9 34 0.225 0.850 13.5
## 6 15-18 18 4 38 0.100 0.950 16.5
## 7 18-21 21 1 39 0.025 0.975 19.5
## 8 21-24 24 1 40 0.025 1.000 22.5Fin qui abbiamo caricato i dati x, creato le classi d’intervallo d’ampiezza 3 e costruito la tabella DF contenente i dati raggruppati con le relative frequenze marginali e cumulate, assolute e relative, nonché il valore centrale di ogni classe d’intervallo.
Adesso utilizziamo la tabella creata, DF, per calcolare gli indici di posizione, variabilità, forma e concentrazione. In particolare, per ciascuna classe di indici, calcoleremo:
posizione: media (aritmetica, geometrica, armonica), mediana, moda, primo e terzo quartile
dispersione: devianza, varianza, scarto quadratico medio, campo di variazione, differenza interquartile, coefficiente di variazione (cv)
forma: indice di asimmetria di Fisher ed indici di curtosi (Pearson e Fisher)
concentrazione: rapporto di concentrazione di Gini, curva di Lorenz e superficie S
# visualizzo frequenze marginali e f. di ripartizione
X = seq(0, 24, by=0.1)
n_x = case_when(
X>=0 & X<=3 ~ p_x[1],
X>3 & X<=6 ~ p_x[2],
X>6 & X<=9 ~ p_x[3],
X>9 & X<=12 ~ p_x[4],
X>12 & X<=15 ~ p_x[5],
X>15 & X<=18 ~ p_x[6],
X>18 & X<=21 ~ p_x[7],
X>21 & X<=24 ~ p_x[8]
)
N_x = case_when(
X>=0 & X<=3 ~ F_x[1],
X>3 & X<=6 ~ F_x[2],
X>6 & X<=9 ~ F_x[3],
X>9 & X<=12 ~ F_x[4],
X>12 & X<=15 ~ F_x[5],
X>15 & X<=18 ~ F_x[6],
X>18 & X<=21 ~ F_x[7],
X>21 & X<=24 ~ F_x[8]
)
df = data.frame(X, n_x, N_x)
ggplot(data=df, aes(x=X, y=n_x)) + geom_line() + labs(x='', y='',title='Funzione di Probabilità') + theme_light()
ggplot(data=df, aes(x=X, y=N_x)) + geom_line() + labs(x='', y='',title='Funzione di Ripartizione') + theme_light()
# INDICI DI POSIZIONE con dati disaggregati
media = mean(x)
media_geometrica = (prod(x))^(1/n)
media_armonica = n/sum(1/x)
quartile1 = quantile(x, 0.25)
mediana = (x[n/2] + x[n/2 +1])/2
quartile3 = quantile(x, 0.75)
# INDICI DI POSIZIONE con dati aggregati in DF
posizione = DF %>% summarise(media = sum(m * n)/sum(n),
media_geom = (prod(m^n))^(1/sum(n)),
media_arm = sum(n)/sum(n/m),
#1quartile appartiene alla classe (6,9)
quartile_1 = 6 + (9-6)*(0.25 - Effe[2])/(Effe[3] - Effe[2]),
#2quartile appartiene alla classe (9,12)
median = 9 + (12-9)*(0.5 - Effe[3])/(Effe[4] - Effe[3]),
#3quartile appartiene alla classe (12,15)
quartile_3 = 12 + 3*(0.75 - Effe[4])/(Effe[5] - Effe[4])
) %>% t()
# confronto posizione
data.frame(row.names=c("media","media geom","media arm", "1°quartile", "mediana","3°quartile"),
disaggregati=c(media, media_geometrica, media_armonica, quartile1, mediana, quartile3), aggregati=posizione)## disaggregati aggregati
## media 10.670000 10.650000
## media geom 9.214724 9.377651
## media arm 6.148860 7.488778
## 1°quartile 8.250000 7.500000
## mediana 10.700000 10.636364
## 3°quartile 12.900000 13.666667# INDICI DI DISPERSIONE con dati disaggregati
devianza = sum((x - mean(x))^2)
varianza = devianza/(n-1)
sd = sqrt(varianza)
campo_variazione = max(x) - min(x)
diff_interquartile = quartile3 - quartile1
cv = sd/abs(media)
# INDICI DI DISPERSIONE con dati aggregati
dispersione = DF %>% summarise(
devianza = sum(n*(m - posizione[1])^2),
varianza = devianza/sum(n),
sd = sqrt(varianza),
campo_variazione = max(m) - min(m),
diff_interquartile = posizione[6] - posizione[4],
cv = sd/abs(posizione[1])
) %>% t()
# confronto dispersione
data.frame(row.names = c("devianza","varianza","scarto quadratico medio", "campo di variazione", "differenza interquartile","coeff variazione"), disaggregati=c(devianza, varianza, sd, campo_variazione, diff_interquartile, cv), aggregati=dispersione)## disaggregati aggregati
## devianza 837.8040000 827.100000
## varianza 21.4821538 20.677500
## scarto quadratico medio 4.6348844 4.547252
## campo di variazione 23.0000000 21.000000
## differenza interquartile 4.6500000 6.166667
## coeff variazione 0.4343847 0.426972# INDICI DI FORMA con dati disaggregati
asimmetria = sum(((x-mean(x))/(sqrt(var(x))))^3)/n
curtosi_beta = sum(((x-mean(x))/(sqrt(var(x))))^4)/n
curtosi_fisher = curtosi_beta - 3 # se >0 asimmetria positiva
# INDICI DI FORMA con dati aggregati
forma = DF %>% summarise(
asimmetria = sum(n*((m - posizione[1])/dispersione[3])^3)/sum(n),
curtosi_beta = sum(n*((m - posizione[1])/dispersione[3])^4)/sum(n),
curtosi_fisher = curtosi_beta - 3
) %>% t()
# confronto forma
data.frame(row.names = c("asimmetria", "curtosi di Pearson", "curtosi di Fisher"),
disaggregati = c(asimmetria, curtosi_beta, curtosi_fisher), aggregati=forma)## disaggregati aggregati
## asimmetria 0.1735845 0.17380070
## curtosi di Pearson 3.2624104 3.02064268
## curtosi di Fisher 0.2624104 0.02064268# INDICI DI CONCENTRAZIONE con dati disaggregati
Q = c()
P = c()
for (i in 1:n) {
Q[i] = sum(x[1:i])/sum(x)
P[i] = i/n
}
# indice di Gini (=0 equadistribuzione, =1 massima concentrazione)
RG = 1- 2*sum(Q[-40])/(n-1)
# oppure in alternativa
RG2 = (sum(P) - sum(Q))/(sum(P)-1)
RG==RG2## [1] FALSE# curva di Lorenz
L = data.frame(P=c(0,P), Q=c(0,Q))
ggplot(L, aes(x=P, y=Q)) + geom_line() + geom_abline(slope=1, intercept=0) + theme_light()
# superficie S compresa tra la bisettrice e la curva di lorenz
# se S=0 equadistribuzione; se S=0.5 massima concentrazione
S = 0.5 - AUC(x=L$P, y=L$Q)
# INDICI DI CONCENTRAZIONE con dati aggregati
Tot = sum(DF$m*DF$n)
Qs = c()
for (i in 1:8) {
Qs[i] = sum(DF$m[1:i]*DF$n[1:i])/Tot
}
sumQ = sum(Qs) - Qs[8]
Ps = c()
for (i in 1:8) {
Ps[i] = sum(DF$n[1:i])/40
}
sumP = sum(Ps) - Ps[8]
RG_aggr = (sumP - sumQ)/sumP
# curva di Lorenz con dati aggregati
L_aggr = data.frame(P=c(0,Ps), Q=c(0,Qs))
ggplot(L_aggr, aes(x=P, y=Q)) + geom_line() + geom_abline(slope=1, intercept=0) + theme_light()
S_aggr = 0.5 - AUC(L_aggr$P, L_aggr$Q)
# confronto concentrazione
data.frame(row.names=c("Gini Ratio", "Area"), disaggregati=c(RG, S), aggregati=c(RG_aggr, S_aggr))## disaggregati aggregati
## Gini Ratio 0.2437339 0.1638438
## Area 0.1188203 0.1177817(X,Y) qualitativi (indice di connessione)
In questa sezione carichiamo i dati relativi a due fenomeni qualitativi, rappresentati congiuntamento con le frequenze multiple (tabelle di contingenza) assolute e relative. Dopo aver caricato i dati, si condurrà un’analisi del grado di associazione-connessione tra queste due variabili. Se non c’è dipendenza statistica, otterremo delle frequenze congiunte (osservate) molto vicine al prodotto delle frequenze marginali (frequenze teoriche, in ipotesi di indipendenza statistica).
# fenomeni X e Y qualitativi (nominali); si osserva congiuntamente (X,Y) = (x,y)
# X = sesso ; Y = churn ; n = 100 individui
n = 100
X = c(rep(1, n/2), rep(0,n/2))
Y = c(rep(1, n/4), rep(0, 3*n/4))
df = data.frame(X, Y)
# frequenza congiunta
Fxy = table(X,Y)
# freq.assoluta marginale di X
mX = c(sum(X==0),sum(X==1))
# freq.assoluta marginale di Y
mY = c(sum(Y==0), sum(Y==1))
# display Fxy, mX e mY insieme
XY = as.matrix(rbind(c(Fxy[1,1],Fxy[1,2],mX[1]), c(Fxy[2,1],Fxy[2,2],mX[2]),c(mY[1],mY[2],n)))
colnames(XY) = c("Y=0","Y=1","mX")
rownames(XY) = c("X=0","X=1","mY")
XY## Y=0 Y=1 mX
## X=0 50 0 50
## X=1 25 25 50
## mY 75 25 100# frequenze relative
fxy = table(X,Y)/n
# freq.relativa di X
fX = mX/n
# freq.relativa di Y
fY = mY/n
# TEST ESATTO DI FISHER per calcolo prob. (a,b,c,d)=(50,0,25,25) | (dato)
# (a+c)=75,(a+b)=50,(b+d)=25,(c+d)=50, n=100
p_Fisher = function(a,b,c,d) {
(factorial(a+b)*factorial(a+c)*factorial(b+d)*factorial(c+d))/(factorial(n)*factorial(a)*factorial(b)*factorial(c)*factorial(d))
}
# si somma la prob. di quanto osservato e le prob. di tutti i casi più estremi
# le marginali NON devono cambiare --> 25,25 | X=0 e 50,0 | X=1
(p_Fisher(50,0,25,25) + p_Fisher(25,25,50,0))*100 #0.000000104% i.e. solo 1 caso ogni 10Mln## [1] 1.042479e-07# poiché la prob. calcolata è bassa, rifiuto l'H0 di X indipendente da Y
# vediamo la stessa ipotesi (X e Y indipendenti) con il TEST CHI-QUADRATO
# frequenze osservate O_ij = fxy_ij
# frequenze attese E_ij = fX_i*fY_j
E = c(fX[1]*fY[1],fX[1]*fY[2],fX[2]*fY[1],fX[2]*fY[2])
O = c(fxy[1,1],fxy[1,2],fxy[2,1],fxy[2,2])
X2 = sum(((O-E)^2)/E) # stat. test ~ (sotto H0) X2(g-1)*(k-1) gdl, i.e. X2(2-1)*(2-1)=X2(1) gdl
# decision rule: se X2 osservato è più estremo del quantile di una X2(1) di livello alpha, rifiuto l'ipotesi nulla di indipendenza a livello 1-alpha.
# confidenza: 1-alpha = 0.95
qAlpha = qchisq(p=0.95, df=1)
# test quantile corretto: area da 0 a qAlpha pari a 0.95
integrale = integrate(function(x) dchisq(x, df=1), lower=0, upper=qAlpha)
# confronto X2 osservato con quantile teorico qAlpha
X2 > qAlpha## [1] FALSE# X2 è minore, quindi non rifiuto H0 con livello di confidenza 95%
# CONCLUSIONE: per Fisher X e Y NON sono indipendenti (rifiuto H0)
# per test Chi-Quadrato, invece, X e Y sono indipendenti (non rifiuto H0)
# le contingenze (differenze tra freq. osservate e teoriche) non sono abbastanza forti per concludere che la distribuzione congiunta non è il prodotto delle marginali
# adesso studiamo (X,Y) caratteri, rispettivamente, ordinali e nominali.
# NON posso usare il test esatto di Fisher (solo per X e Y dicotomiche), ma posso usare
# il test Chi-Quadrato per avere una misura di connessione (se misura bassa -> X ind. da Y)
# X = voti da 1 a 6 di 1000 studenti all'esame di matematica
# Y = laurea (5 anni dopo) in facoltà scientifica (=1) o umanistica (=0)
n=1000
X = c(rep(1,0.2*n),rep(2,0.2*n),rep(3,0.2*n),rep(4,0.2*n),rep(5,0.1*n),rep(6,0.1*n))
Y = c(rep(0,n/8), rep(1,n/8), rep(0,3*n/8), rep(1, 3*n/8))
table(X,Y)## Y
## X 0 1
## 1 125 75
## 2 150 50
## 3 200 0
## 4 25 175
## 5 0 100
## 6 0 100# calcoliamo frequenze relative, congiunte e marginali
fX = c()
for (i in 1:6) {
fX[i] = sum(X==i)/n
}
options(scipen=999) # disabilita notazione scientifica
fY = c()
for (i in 0:1) {
fY[i+1] = sum(Y==i)/n
}
fxy = table(X,Y)/n
XY_rel = rbind(c(fxy[1,1],fxy[1,2],fX[1]),c(fxy[2,1],fxy[2,2],fX[2]),c(fxy[3,1],fxy[3,2],fX[3]),
c(fxy[4,1],fxy[4,2],fX[4]),c(fxy[5,1],fxy[5,2],fX[5]),c(fxy[6,1],fxy[6,2],fX[6]),
c(fY[1],fY[2],1))
colnames(XY_rel) = c("Y=0","Y=1","mX")
rownames(XY_rel) = c("X=1","X=2","X=3","X=4","X=5","X=6","mY")
XY_rel## Y=0 Y=1 mX
## X=1 0.125 0.075 0.2
## X=2 0.150 0.050 0.2
## X=3 0.200 0.000 0.2
## X=4 0.025 0.175 0.2
## X=5 0.000 0.100 0.1
## X=6 0.000 0.100 0.1
## mY 0.500 0.500 1.0XY_ass = XY_rel*n
# costruisco frequenze teoriche E = prodotto delle freq. marginali relative
E = fX %*% t(fY)
O = XY_rel[1:6,1:2]
X2 = sum(((O-E)^2/E))
# gdl = (6-1)*(2-1)=5*1=5
# qAlpha da una X2(5)
X2 > qchisq(p=0.95, df=5)## [1] FALSE# la stat test è minore del quantile, quindi NON rifiuto l'H0 di indipendenza
# calcolo V di Cramer (normalizzo X2)
V = sqrt(X2/min(ncol(XY_rel) -1, nrow(XY_rel) -1))
#V=1 max dipendenza, V=0 indipendenza assoluta (join=prod marginali)
#V=0.5362 indica modesta dipendenza, abbastanza per rifiutare H0:X ind. Y
V## [1] 0.5361903(X,Y) quantitativi (indice di correlazione)
In questa sezione carichiamo i dati relativi a due fenomeni quantitativi (X,Y) e cerchiamo di misurare il grado di associazione tra questi due fenomeni, ricorrendo al concetto di correlazione. In generale due fenomeni quantitativi si dicono correlati se all’aumentare di un fenomeno ci aspettiamo in media un aumento o una diminuizione anche dell’altro fenomeno.
# modello di regressione lineare univariato: Y=f(X)+e con f() lineare, e errore white noise
n = 100
X = rnorm(n, 10, sd=1.5)
f = function(x) {4.8*X - 2}
e = rnorm(n, 0, 1)
Y = f(X) + e
# visualizzo scatter-plot di Y,X e sovraimpongo funzione di regressione
df = data.frame(X,Y)
ggplot(df, aes(x=X, y=Y)) + geom_point() + geom_smooth(method='lm', se=F) + theme_light()## `geom_smooth()` using formula 'y ~ x'
# misura analitica della correlazione tra X,Y e del coefficiente di corr.lineare
correlazione = cor(X,Y)
correlazione # 0.99 è molto vicino a +1; fortissima correlazione lineare positiva## [1] 0.9900239# dimostro che R2 è uguale a correlazione^2 nel caso di modello lineare univariato
mod = lm(Y ~ 1+X, data=df)
# def funzione r2; rss=residual sum of squares; tss=total sum of squares
rsquared = function(rss, tss) {1 - rss/tss}
summary(mod) # intercetta stimata -1.62; intercetta reale -2##
## Call:
## lm(formula = Y ~ 1 + X, data = df)
##
## Residuals:
## Min 1Q Median 3Q Max
## -3.07693 -0.75053 -0.01544 0.76566 1.75755
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -2.35444 0.69549 -3.385 0.00102 **
## X 4.82435 0.06936 69.558 < 0.0000000000000002 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.9922 on 98 degrees of freedom
## Multiple R-squared: 0.9801, Adjusted R-squared: 0.9799
## F-statistic: 4838 on 1 and 98 DF, p-value: < 0.00000000000000022 # coeffX stimato +4.76; reale +4.8
tss = sum((Y - mean(Y))^2)
rss = sum((mod$residuals)^2)
data.frame(r2=rsquared(rss=rss, tss=tss), rho2 = correlazione^2)## r2 rho2
## 1 0.9801473 0.9801473# calcolo manualmente i coefficienti di regressione beta_ols
N = cbind(rep(1,n),X)
beta_ols = solve(t(N)%*%N)%*%t(N)%*%Y