\(P(A|B) = \frac{P(A \cap B)}{P(B)} \qquad \text{cioè} \qquad P(A \cap B) = P(A|B)P(B)\)
\(P(A) = P(B)P(A|B) + P(\bar{B})P(A|\bar{B}) = P(A \cap B) + P(A \cap \bar{B})\)
\(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)}\)
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\(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)\)
\({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.02642412dbinom(5,10,0.2)[1] 0.02642412E così anche per Poisson, Normale, esponenziale (ed altre)
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 discretaUna 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.
somma <- sum(x)
somma_cumulata <- cumsum(x)
n <- length(x)
media_campionaria <- somma / n
mu <- media_campionaria
mu[1] 6.06mean(x)[1] 6.06n ^ -1 * sum((x - mu) ^ 2)[1] 8.0564var(x)[1] 8.137778(n - 1) ^ -1 * sum((x - mu) ^ 2)[1] 8.137778var(x)[1] 8.137778X <- seq(1, 10, 1)
Cn <- c()
for (i in X) {
Cn <- c(Cn, n ^ -1 * sum(x <= i))
}
plot(Cn,type = "h",xlab="x")
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")
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.9Diminuiamo 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")
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")
\(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")
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.
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)
(vedi inizio script 5a lezione)
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 1Analisi 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-07Diagnostica 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-12Diagnostica 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