Statistica_6alezione

Author

Alessandro Fusta Moro

RIPASSO

Probabilità condizionata

\(P(A|B) = \frac{P(A \cap B)}{P(B)} \qquad \text{cioè} \qquad P(A \cap B) = P(A|B)P(B)\)

Principio delle probabilità totali

\(P(A) = P(B)P(A|B) + P(\bar{B})P(A|\bar{B}) = P(A \cap B) + P(A \cap \bar{B})\)

Teorema di Bayes

\(P(A|B) = \frac{P(A)P(B|A)}{P(B)} \qquad \text{da cui} \qquad P(A) = \frac{P(A|B)P(B)}{P(B|A)}\)

Esercizio

Una macchina produce pezzi difettosi (evento B) con probabilità \(P(B)= 0.001\). Una procedura di controllo scarta pezzi difettosi (evento A condizionato da B) con probabilità \(P(A|B) = 0.9\) mentre scarta pezzi non difettosi con probabilità \(P(A|\bar{B})= 0.01\). Qual’è la probabilità \(P(B|A)\) che un pezzo sia non difettoso dato che è stato scartato?

Dal principio delle probabilità totali (sostituisco P(A)):

\(P(B)P(A|B) + P(\bar{B})P(A|\bar{B}) = \frac{P(A|B)P(B)}{P(B|A)}\)

\(P(B|A) = \frac{P(A|B)P(B)}{P(B)P(A|B) + P(\bar{B})P(A|\bar{B})}\)

rm(list = ls())
P_B = 0.001 #P(B) probabilità di produrre pezzi difettosi
P_AcondB = 0.9 #P(A|B) la prob di scartarlo quando il pezzo è difettoso
P_AcondcompB = 0.01 #P(A|compB) la prob di scartare non difettosi

#probabilità che il pezzo scartato fosse difettoso = P(B|A)
P_BcondA <-
  P_AcondB * P_B / ((P_B * P_AcondB) + ((1 - P_B) * P_AcondcompB))
P_BcondA # =~ 8 %
[1] 0.08264463

Indipendenza in probabilità

\(P(A \cap B) = P(A)P(B)\)

\(P(A|B) = \frac{P(A \cap B)}{P(B)} = \frac{P(A)P(B)}{P(B)} = P(A)\)

Modelli statistici

\({p(x,θ): θ ∈ Θ,x ∈ S}\) per variabili casuali discrete \({f(x,θ): θ ∈ Θ,x ∈ S}\) per variabili casuali continue

Esempio di un modello noto: binomiale

bin <- function(n, k, p) {
  bin <- choose(n,k)*
    p ^ (k) * (1 - p) ^ (n - k) # con 0 < p < 1 e x=(0,1)
  return(bin)
}
bin(10,5,0.2)
[1] 0.02642412
dbinom(5,10,0.2)
[1] 0.02642412

E così anche per Poisson, Normale, esponenziale (ed altre)

Campione casuale

Come generare un campione casuale dato un modello. Le variabili casuali \(X_i \text{ con } i = 1,\ldots,n\) servono per modellare le osservazioni campionarie \(x_i \text{ con } i = 1,\ldots,n\). L’inferenza statistica consiste nel fare congetture sull’incognito vettore di parametri \(\theta = [\theta_1,\ldots,\theta_p]\) basandosi su una realizzazione campionaria \([x_1,\ldots,x_n]\).

La definizione di campione casuale da un modello statistico formalizza che lo stesso fenomeno ripetibile con risultato incerto viene osservato n volte (le n variabili casuali campionarie con la stessa funzione di ripartizione) e che le n osservazioni sono ottenute in modo in- dipendente (la condizione di indipendenza in probabilità delle n variabili casuali campionarie).

set.seed(1)
x <- sample(1:10, 100, replace = T) #da un uniforme discreta

Una funzione delle sole variabili casuali campionarie è chiamata statistica. Le statistiche sono uno degli strumenti di base utilizzati nell’ineferenza statistica. Se usate in un problema di stima puntuale le statistiche prendono il nome di stimatori se usate in un problema di verifica di ipotesi vengono chiamate statistiche test.

Stima puntuale

Media campionaria

somma <- sum(x)
somma_cumulata <- cumsum(x)
n <- length(x)
media_campionaria <- somma / n
mu <- media_campionaria
mu
[1] 6.06
mean(x)
[1] 6.06

Varianza campionaria

n ^ -1 * sum((x - mu) ^ 2)
[1] 8.0564
var(x)
[1] 8.137778

Varianza campionaria corretta

(n - 1) ^ -1 * sum((x - mu) ^ 2)
[1] 8.137778
var(x)
[1] 8.137778

frequenza cumulata relativa campionaria

X <- seq(1, 10, 1)
Cn <- c()
for (i in X) {
  Cn <- c(Cn, n ^ -1 * sum(x <= i))
}
plot(Cn,type = "h",xlab="x")

frequenza relativa campionaria

Fn <- c()
for (i in X) {
  Fn <- c(Fn, n ^ -1 * sum(x == i))
}
plot(Fn,ylim = c(0,0.3),type = "h",xlab="x")

Funzione di verosimiglianza

La funzione di probabilità congiunta delle variabili casuali campionarie è definitita come:

\(V(x_1,x_2,\dots,x_n) = p(x_1,\theta) \times p(x_2,\theta) \times \ldots \times p(x_n,\theta) = \prod_{i=1}^n p(x_i,\theta)\)

ovvero il prodotto delle singole probabilità. Vediamo un esempio:

V<-c() #initialization
theta.s <- seq(.1,.9,.01)
for (theta in theta.s) {
  V <- c(V,prod(rep(bin(1,1,theta),50),#successi
                rep(bin(1,0,theta),50)#insuccessi
                ))
  # print(V[length(V)])
}
plot(y=V,x=theta.s, xlab="theta",type = "l")

Aumentiamo il numero dei successi

V<-c() #initialization
theta.s <- seq(.1,.9,.01)
for (theta in theta.s) {
  V <- c(V,prod(rep(bin(1,1,theta),90),#successi
                rep(bin(1,0,theta),10)#insuccessi
                ))
  # print(V[length(V)])
}
plot(y=V,x=theta.s, xlab="theta",type = "l")

90/100
[1] 0.9

Diminuiamo la numerosità del campione

V<-c() #initialization
theta.s <- seq(.1,.9,.01)
for (theta in theta.s) {
  V <- c(V,prod(rep(bin(1,1,theta),10),#successi
                rep(bin(1,0,theta),10)#insuccessi
                ))
  # print(V[length(V)])
}
plot(y=V,x=theta.s, xlab="theta",type = "l")

Log-verosimiglianza

Spesso si usa la versione logaritmica per rendere il calcolo più agevole, chiamata funzione di log-verosimiglianza: \(L(\theta) = ln(p(x_1,\theta) \times \ldots \times p(x_n,\theta)) = ln(p(x_1,\theta) + \ldots + ln(p(x_n,\theta) = \sum_{i=1}^n ln(p(x_i,n))\)

L<-c() #initialization
theta.s <- seq(.1,.9,.01)
for (theta in theta.s) {
  L <- c(L,sum(log(rep(bin(1,1,theta),50)),#successi
                log(rep(bin(1,0,theta),50))#insuccessi
                ))
  # print(V[length(V)])
}
plot(y=L,x=theta.s, xlab="theta",type = "l")

log-verosimiglianza

Binomiale

\(L(p)=\sum_{i=1}^n x_i ln(p) + (n - \sum_{i=1}^n x_i) ln(1-p)\)

L<-c() #initialization
theta.s <- seq(.1,.9,.01)
suc <- 50 #successi
N <- 100
for (theta in theta.s) {
L <- c(L,sum(rep(log(theta),suc))+sum(rep(log(1-theta),N-suc)))
}
plot(y=L,x=theta.s, xlab="theta",type = "l")

Esponenziale

L<-c() #initialization
theta.s <- seq(.1,5,.1)
x <- rexp(100)
for (theta in theta.s) {
L <- c(L,length(x)*log(theta) - (theta*sum(x)))
}
plot(y=L,x=theta.s, xlab="theta",type = "l")

e così via per Poisson, Normale ecc.

Proprietà delle Stime Puntuale

Correttezza

Il valore atteso deve coincidere con la stima ### Consistenza Ovvero la probabilità di un errore di stima grossolano deve tendere a zero al di- vergere della numerosità campionaria n ### Efficienza La varianza dev’essere la minore possibile (per approf Teorema di Cramer-Rao)

La media campionaria è uno stimatore corretto e consistente!

(vedi inizio script 5a lezione)

Regressione lineare

Caricare un dataset in R

data(mtcars)
head(mtcars)
                   mpg cyl disp  hp drat    wt  qsec vs am gear carb
Mazda RX4         21.0   6  160 110 3.90 2.620 16.46  0  1    4    4
Mazda RX4 Wag     21.0   6  160 110 3.90 2.875 17.02  0  1    4    4
Datsun 710        22.8   4  108  93 3.85 2.320 18.61  1  1    4    1
Hornet 4 Drive    21.4   6  258 110 3.08 3.215 19.44  1  0    3    1
Hornet Sportabout 18.7   8  360 175 3.15 3.440 17.02  0  0    3    2
Valiant           18.1   6  225 105 2.76 3.460 20.22  1  0    3    1

Analisi Esplorativa dei Dati

# Codice
summary(mtcars)
      mpg             cyl             disp             hp       
 Min.   :10.40   Min.   :4.000   Min.   : 71.1   Min.   : 52.0  
 1st Qu.:15.43   1st Qu.:4.000   1st Qu.:120.8   1st Qu.: 96.5  
 Median :19.20   Median :6.000   Median :196.3   Median :123.0  
 Mean   :20.09   Mean   :6.188   Mean   :230.7   Mean   :146.7  
 3rd Qu.:22.80   3rd Qu.:8.000   3rd Qu.:326.0   3rd Qu.:180.0  
 Max.   :33.90   Max.   :8.000   Max.   :472.0   Max.   :335.0  
      drat             wt             qsec             vs        
 Min.   :2.760   Min.   :1.513   Min.   :14.50   Min.   :0.0000  
 1st Qu.:3.080   1st Qu.:2.581   1st Qu.:16.89   1st Qu.:0.0000  
 Median :3.695   Median :3.325   Median :17.71   Median :0.0000  
 Mean   :3.597   Mean   :3.217   Mean   :17.85   Mean   :0.4375  
 3rd Qu.:3.920   3rd Qu.:3.610   3rd Qu.:18.90   3rd Qu.:1.0000  
 Max.   :4.930   Max.   :5.424   Max.   :22.90   Max.   :1.0000  
       am              gear            carb      
 Min.   :0.0000   Min.   :3.000   Min.   :1.000  
 1st Qu.:0.0000   1st Qu.:3.000   1st Qu.:2.000  
 Median :0.0000   Median :4.000   Median :2.000  
 Mean   :0.4062   Mean   :3.688   Mean   :2.812  
 3rd Qu.:1.0000   3rd Qu.:4.000   3rd Qu.:4.000  
 Max.   :1.0000   Max.   :5.000   Max.   :8.000  
str(mtcars)
'data.frame':   32 obs. of  11 variables:
 $ mpg : num  21 21 22.8 21.4 18.7 18.1 14.3 24.4 22.8 19.2 ...
 $ cyl : num  6 6 4 6 8 6 8 4 4 6 ...
 $ disp: num  160 160 108 258 360 ...
 $ hp  : num  110 110 93 110 175 105 245 62 95 123 ...
 $ drat: num  3.9 3.9 3.85 3.08 3.15 2.76 3.21 3.69 3.92 3.92 ...
 $ wt  : num  2.62 2.88 2.32 3.21 3.44 ...
 $ qsec: num  16.5 17 18.6 19.4 17 ...
 $ vs  : num  0 0 1 1 0 1 0 1 1 1 ...
 $ am  : num  1 1 1 0 0 0 0 0 0 0 ...
 $ gear: num  4 4 4 3 3 3 3 4 4 4 ...
 $ carb: num  4 4 1 1 2 1 4 2 2 4 ...
plot(mtcars$hp, mtcars$mpg, main="HP vs MPG", xlab="Horsepower", ylab="Miles per Gallon")

plot(mtcars$wt, mtcars$mpg, main="Weight vs MPG", xlab="Weight", ylab="Miles per Gallon")

Creazione del Modello di Regressione Lineare

model <- lm(mpg ~ hp, data=mtcars)
summary(model)

Call:
lm(formula = mpg ~ hp, data = mtcars)

Residuals:
    Min      1Q  Median      3Q     Max 
-5.7121 -2.1122 -0.8854  1.5819  8.2360 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 30.09886    1.63392  18.421  < 2e-16 ***
hp          -0.06823    0.01012  -6.742 1.79e-07 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3.863 on 30 degrees of freedom
Multiple R-squared:  0.6024,    Adjusted R-squared:  0.5892 
F-statistic: 45.46 on 1 and 30 DF,  p-value: 1.788e-07

Diagnostica dei residui

plot(model)

Regressione Lineare Multipla

model2 <- lm(mpg ~ hp + wt, data=mtcars)
summary(model2)

Call:
lm(formula = mpg ~ hp + wt, data = mtcars)

Residuals:
   Min     1Q Median     3Q    Max 
-3.941 -1.600 -0.182  1.050  5.854 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 37.22727    1.59879  23.285  < 2e-16 ***
hp          -0.03177    0.00903  -3.519  0.00145 ** 
wt          -3.87783    0.63273  -6.129 1.12e-06 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 2.593 on 29 degrees of freedom
Multiple R-squared:  0.8268,    Adjusted R-squared:  0.8148 
F-statistic: 69.21 on 2 and 29 DF,  p-value: 9.109e-12

Diagnostica dei residui

plot(model2)

Predizione con il Modello

newdata <- data.frame(hp=150, wt=3.2)
predict(model2, newdata)
       1 
20.05227 

Validazione del Modello

set.seed(123)  # Per riproducibilità
sample_index <- sample(seq_len(nrow(mtcars)), size = 0.7*nrow(mtcars))
train <- mtcars[sample_index, ]
test <- mtcars[-sample_index, ]

model_train <- lm(mpg ~ hp + wt, data=train)
predictions <- predict(model_train, test)

# Calcolo dell'errore medio quadratico (MSE)
mse <- mean((test$mpg - predictions)^2)
mse
[1] 4.242533