Esercizi

Autore/Autrice

Paolo Bosetti

Data di Pubblicazione

11 ottobre 2025

Data di modifica

4 dicembre 2025

Sommario

Gli esercizi sono riportati di seguito, prima elencando le domande, poi (a fine documento), gli svolgimenti. I numeri delle domande della prima sezione corrispondeono ai numeri delle risposte nella seconda.

1 Domande

1.1 Test di Student

Si hanno due campioni da 9 e 12 elementi. Il primo ha i valori 12, 11.7, 9.7, 11, 12.5, 12.7, 9.9, 11.4, 9.2; il secondo 12.5, 15.1, 14.4, 12.5, 14.9, 14.4, 13.2, 11.2, 12.6, 14.8, 13.9, 11.9. Assumendo che i campioni abbiano la stessa varianza, effettuare un test di Student a due lati senza utilizzare la funzione t.test().

1.2 Potenza del test di Student

Consideriamo il secondo campione dell’esercizio precedente. Risulta \(\bar y_2 = 13.45\). Verificare con la funzione t.test() l’ipotesi che il valore atteso della popolazione d’origine sia diverso da \(\mu_2 = 12.7\). Dato che nella fattispecie risulta \(\bar y_2 > \mu_2\), effettuare anche un test a un lato e discutere il risultato in termini di potenza del test.

Qual è il valore minimo di \(\mu_2\) che ci porta a rigettare \(H_0\) con una probabilità d’errore minore del 1%?

1.3 Manipolazione dati

Considerare il dataset ggplot2::diamonds che riporta il prezzo di vendita di diamanti in funzione delle loro caratteristiche morfologiche.

Raggruppare la dimensione in carati (carat) in cinque gruppi di caratura. Costruire una tabella che riporti, per ciascun gruppo e per ciascuna categoria di taglio (cut), il minimo, il massimo e i tre quartili del prezzo di vendita, in colonne separate.

1.4 Dati di produzione

La tabella dati prod_defects.csv sul sito http://paolobosetti.quarto.pub/data.html#list-of-csv-files contiene dei dati di difettosità (come numero difetti per settimana) su 50 settimane di un impianto di produzione con tre macchine, M1, M2 e M3. Per le tre macchine i limiti massimi di difettosità sono, rispettivamente, 27, 32 e 100.

Creare una tabella che riporti, per ognuna delle 50 settimane, la difettosità in percentuale sul massimo per ciascuna macchina e realizzare il corrispondente grafico.

Per ogni macchina, aggiungere anche il calcolo della media progressiva, cioè il valor medio dei campioni dal primo all’i-esimo.

Nota

Ricordare che tutte le tabelle presenti in http://paolobosetti.quarto.pub/data.html possono essere caricate direttamente dalle funzioni read_table() o read_csv() passando l’URL generato come adas.utils::examples_url("nome_del_file.csv").

2 Risposte

2.1 Test di Student

c1 <- c(12, 11.7, 9.7, 11, 12.5, 12.7, 9.9, 11.4, 9.2)
c2 <- c(12.5, 15.1, 14.4, 12.5, 14.9, 14.4, 13.2, 11.2, 12.6, 14.8, 13.9, 11.9)

c1.m <- mean(c1)
c2.m <- mean(c2)
n1 <- length(c1)
n2 <- length(c2)
Sp <- sqrt(((n1-1)*var(c1) + (n2-1)*var(c2)) / (n1 + n2 -2))

t0 <- abs(c1.m - c2.m) / (Sp * sqrt(1/n1 + 1/n2))
p.value <- pt(t0, n1 + n2 -2, lower.tail=FALSE) * 2
cat(paste("t0:", t0, "\n"))
t0: 4.09946306243968 
cat(paste("p-value:", p.value))
p-value: 0.00061049441430277

Per verifica:

t.test(c1, c2, var.equal = TRUE)

    Two Sample t-test

data:  c1 and c2
t = -4.0995, df = 19, p-value = 0.0006105
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
 -3.516249 -1.139306
sample estimates:
mean of x mean of y 
 11.12222  13.45000

2.2 Potenza del test di Student

A due lati:

t.test(c2, mu=12.7)

    One Sample t-test

data:  c2
t = 1.9921, df = 11, p-value = 0.07177
alternative hypothesis: true mean is not equal to 12.7
95 percent confidence interval:
 12.62136 14.27864
sample estimates:
mean of x 
    13.45

A un lato:

(t <- t.test(c2, mu=12.7, alternative="greater", conf.level = 0.99))

    One Sample t-test

data:  c2
t = 1.9921, df = 11, p-value = 0.03588
alternative hypothesis: true mean is greater than 12.7
99 percent confidence interval:
 12.42668      Inf
sample estimates:
mean of x 
    13.45

Adottando il test a un lato, il p-value risulta inferiore. Se ci fossimo dati una soglia per la probabilità di un falso positivo pari a \(\alpha=0.05\), in particolare, nel primo caso non avremmo potuto rigettare \(H_0\), mentre nel secondo sì.

Osservando l’intervallo di confidenza al 99%, valori di \(mu_2\) inferiori a 12.4266787 consentirebbero di rifiutare l’ipotesi nulla con una probabilità d’errore inferiore all’1%.

Nota

Provare a modificare il valore del parametro conf.level e osservare come l’unico risultato nell’output di t.test() è l’ampiezza dell’intervallo di confidenza.

2.3 Manipolazione dati

diamonds %>% 
  mutate(size = cut(carat, breaks = 5), .after=carat) %>% 
  select(size:clarity, price) %>% 
  group_by(size,cut) %>% 
  summarise(
    probs = seq(0, 100, length.out=5),
    price = quantile(price),
    .groups="keep"
  ) %>% 
  mutate(probs = paste0(probs, "%")) %>% 
  pivot_wider(names_from = probs, values_from = price) %>% 
  knitr::kable()
size cut 0% 25% 50% 75% 100%
(0.195,1.16] Fair 337 1715.00 2684.5 3785.25 10752
(0.195,1.16] Good 327 896.50 2258.0 4009.50 17499
(0.195,1.16] Very Good 336 773.00 1975.0 3937.00 18542
(0.195,1.16] Premium 326 882.00 1826.5 4033.00 18279
(0.195,1.16] Ideal 326 827.00 1389.0 2947.00 17590
(1.16,2.12] Fair 2336 5334.00 7388.0 10913.00 18574
(1.16,2.12] Good 3098 6686.50 8820.0 12236.00 18707
(1.16,2.12] Very Good 2774 6945.25 9526.0 12762.75 18818
(1.16,2.12] Premium 2360 6604.75 9384.5 12751.50 18795
(1.16,2.12] Ideal 3168 7155.50 9761.0 12735.50 18806
(2.12,3.09] Fair 5405 7069.00 12500.0 14975.00 18308
(2.12,3.09] Good 9068 13751.50 14613.5 16894.25 18788
(2.12,3.09] Very Good 6512 13782.00 15753.0 17161.00 18692
(2.12,3.09] Premium 6357 13717.00 15857.5 17262.75 18823
(2.12,3.09] Ideal 8709 14548.75 16427.0 17513.00 18791
(3.09,4.05] Fair 9823 10745.50 11668.0 13816.00 15964
(3.09,4.05] Very Good 15984 15984.00 15984.0 15984.00 15984
(3.09,4.05] Premium 12300 15223.00 15223.0 16193.00 18701
(3.09,4.05] Ideal 12545 12555.50 12566.0 12576.50 12587
(4.05,5.01] Fair 17329 17673.50 18018.0 18274.50 18531
NotaOsservazioni
  • size:clarity: in dplyr è possibile selezionare più colonne dando un intervallo di nomi di colonne.
  • summarise() può dare un warning suggerendo di usare reframe: durante il corso e durante l’esame NON utilizzate reframe(), perché essendo sperimentale non ne è garantita la disponibilità sui PC usati per l’esame.
  • Quando le funzioni sommario utilizzate in summarise() ritornano più di un valore, il sommario riporta altrettante righe. In questo caso come sommario utilizziamo quantile(), che ritorna 5 valori, e otteniamo quindi i valori desiderati, sebbene in un’unica colonna. Aggiungiamo la colonna perc per avere le corrispondenti percentuali.
  • Ridistribuiamo poi i valori in colonne mediante pivot_wider().

2.4 Dati di produzione

Osserviamo il file di dati originale:

examples_url("prod_defects.csv") %>% 
  read_file() %>% 
  substr(1,250) %>% # solo i primi 250 caratteri
  cat()
# Difetti in produzione su 50 settimane
# S: numero settimana
# M1: conteggio difetti per macchina 1
# M2: conteggio difetti per macchina 2
# M3: conteggio difetti per macchina 3
S,M1,M2,M3
1,3,1,2
2,1,0,2
3,1,1,1
4,2,0,1
5,3,0,5
6,0,0,2
7,3,0,0
8,4,

È evidente che il file contiene dei commenti inizianti con # e che i dati non sono tidy.

Importiamo quindi la tabella e trasformiamola in formato tidy:

prod <- examples_url("prod_defects.csv") %>% 
  read_csv(comment = "#") %>% 
  pivot_longer(-S, names_to="Macchina", values_to="Difetti") %>% 
  rename(Settimana=S)

Visualizziamo i dati:

prod %>% 
  ggplot(aes(x=Settimana,y=Difetti, color=Macchina)) + 
  geom_point()

Creiamo una tabella di appoggio per i limiti:

limits <- tribble(
  ~Macchina, ~max,
  "M1", 27,
  "M2", 32,
  "M3", 100
)

Uniamo le due tabelle e calcoliamo il rapporto:

prod %>% 
  left_join(limits, by = join_by(Macchina)) %>% 
  mutate(`Dif.rel.`=(Difetti/max*100) %>% round(2)) %>% 
  ggplot(aes(x=Settimana, y=`Dif.rel.`, color=Macchina)) +
  geom_line()

Per il calcolo della media progressiva possiamo usare una mappa, calcolando la media dei valori da 1 a \(i\):

examples_url("prod_defects.csv") %>% 
  read_csv(comment = "#") %>%
  mutate(
    M1.m = map_dbl(1:n(), \(i) mean(M1[1:i])),
    M2.m = map_dbl(1:n(), \(i) mean(M2[1:i])),
    M3.m = map_dbl(1:n(), \(i) mean(M3[1:i]))
  ) %>% 
  slice_head(n=10)
Rows: 50 Columns: 4
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
dbl (4): S, M1, M2, M3

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
# A tibble: 10 × 7
       S    M1    M2    M3  M1.m  M2.m  M3.m
   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
 1     1     3     1     2  3    1      2   
 2     2     1     0     2  2    0.5    2   
 3     3     1     1     1  1.67 0.667  1.67
 4     4     2     0     1  1.75 0.5    1.5 
 5     5     3     0     5  2    0.4    2.2 
 6     6     0     0     2  1.67 0.333  2.17
 7     7     3     0     0  1.86 0.286  1.86
 8     8     4     0     0  2.12 0.25   1.62
 9     9     2     0     1  2.11 0.222  1.56
10    10     2     1     3  2.1  0.3    1.7

Operare per colonne, però, non è molto pulito: supponiamo che in un futuro l’impianto cresca a 20 macchine, sarebbe necessario ripetere la stessa operazione più volte. È meglio risolvere il problema in formato tidy,

Per farlo partiamo dalla tabella tidy (prod). Non possiamo applicare direttamente una mappa, perché otterremo la media progressiva di tutti i difetti mescolando le varie macchine.

Possiamo però sfruttare group_by, che consente di operare su una tabella per gruppi:

prod %>% 
  group_by(Macchina) %>% 
  mutate(
    mean = map_dbl(1:n(), ~ mean(Difetti[1:.]))
  ) %>% 
  slice_head(n=3)
# A tibble: 9 × 4
# Groups:   Macchina [3]
  Settimana Macchina Difetti  mean
      <dbl> <chr>      <dbl> <dbl>
1         1 M1             3 3    
2         2 M1             1 2    
3         3 M1             1 1.67 
4         1 M2             1 1    
5         2 M2             0 0.5  
6         3 M2             1 0.667
7         1 M3             2 2    
8         2 M3             2 2    
9         3 M3             1 1.67
Importante

Nonostante slice_head(n=3) otteniamo 9 righe: come mai? è dovuto al fatto che la tabella è ancora raggruppata quindi ogni operazione si applica un gruppo alla volta, ottenendo tre righe per ognuno dei tre gruppi.

Il raggruppamento può essere disabilitato con ungroup().

Ora mettiamo tutto assieme, con l’attenzione di calcolare la media progressiva per la difettosità relativa (e non per il conteggio):

prod %>% 
  # Calcolo difettosità relativa
  left_join(limits, by = join_by(Macchina)) %>% 
  mutate(`Dif.rel.`=(Difetti/max*100) %>% round(2)) %>% 
  # Calcolo le medie progressive per gruppi
  group_by(Macchina) %>% 
  mutate(
    mean = map_dbl(1:n(), ~ mean(`Dif.rel.`[1:.]))
  ) %>% 
  ungroup() %>% 
  # Grafico
  ggplot(aes(x=Settimana, color=Macchina)) +
  geom_line(aes(y=`Dif.rel.`)) +
  geom_line(aes(y=mean), linetype=2)

3 Regressione

Carichiamo i dati e dividiamoli già in due gruppi training e validation (80:20):

set.seed(0)
df <- examples_url("regression.csv") %>% 
  read_csv(show_col_types = F) %>% 
  mutate(train = sample(c(T,F), size=n(), replace=T, prob=c(0.8, 0.2)))

df %>% 
  ggplot(aes(x=x, y=y)) +
  geom_point(aes(color=train))

Confrontiamo due regressioni lineari, una col modello \(y=a + b\log(x)\), una col modello \(y=a + b/x^{0.1}\):

df %>% filter(train) %>% 
  ggplot(aes(x=x, y=y)) +
  geom_point() +
  geom_smooth(
    method = "lm",
    formula = y ~ I(x^0.1)
  ) +
  geom_smooth(
    method = "lm",
    formula = y ~ log(x),
    color="red"
  )

summary(lm(y~log(x), data=df))

Call:
lm(formula = y ~ log(x), data = df)

Residuals:
    Min      1Q  Median      3Q     Max 
-1.3905 -0.5088 -0.1289  0.5283  1.8303 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 12.78895    0.24923   51.31   <2e-16 ***
log(x)       2.88642    0.09549   30.23   <2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.73 on 48 degrees of freedom
Multiple R-squared:  0.9501,    Adjusted R-squared:  0.949 
F-statistic: 913.6 on 1 and 48 DF,  p-value: < 2.2e-16
summary(lm(y~I(x^0.1), data=df))

Call:
lm(formula = y ~ I(x^0.1), data = df)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.42760 -0.53774 -0.04605  0.49055  2.13253 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -12.3808     1.0906  -11.35  3.4e-15 ***
I(x^0.1)     25.1179     0.8513   29.50  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.7469 on 48 degrees of freedom
Multiple R-squared:  0.9477,    Adjusted R-squared:  0.9467 
F-statistic: 870.5 on 1 and 48 DF,  p-value: < 2.2e-16

Difficile dire quale sia meglio, anche se \(R^2\) del modello logaritmico è leggermente più elevato. Proviamo con la cross-validazione, confrontando diverse potenze oltre al logaritmo:

library(modelr)
models <- c(
  map(1:15, \(n) lm(y~I(x^(1/n)), data=filter(df, train))),
  list(lm(y~log(x), data=filter(df, train)))
)

tibble(
  n = seq_along(models),
  train = map_dbl(models, ~ rmse(., filter(df, train))),
  valid = map_dbl(models, ~ rmse(., filter(df, !train)))
) %>% 
  pivot_longer(-n, names_to = "type", values_to = "RMSE") %>% 
  ggplot(aes(x=n, y=RMSE, fill=type)) + 
  geom_col(position = "dodge")

L’ultima barra corrisponde al modello logaritmico. Si osserva che RMSE di tale modello per i dati di training e per quelli di valdazione è effettivamente il minimo, anche se l’esponente \(1/15\) produce risultati ancora migliori di \(1/10\). Tutavia, per \(n\gt10\) non si osservano sensibili miglioramenti per l’RMSE sui dati di validazione.

Osserviamo infine i residui per verificare l’assenza di pattern:

df %>% 
  add_residuals(lm(y~log(x), data=df), var="log") %>% 
  add_residuals(lm(y~I(x^0.1), data=df), var="pow") %>% 
  select(x, log, pow) %>% 
  pivot_longer(-x, names_to = "model", values_to = "residual") %>% 
  ggplot(aes(x=x, y=residual, color=model)) + 
  geom_point()

In entrambi i casi non si notano pattern significativi. Dal punto di vista fisico, è più frequente trovare relazioni logaritmiche che del tipo \(y\propto x^{0.1}\), quindi accettiamo il modello logaritmico.

Nota

Il grafico di RMSE in funzione del grado sopra riportato non mostra un minimo perché stiamo confrontando lo stesso modello con parametri (l’esponente) differenti e non diversi modelli. Quindi, il modello esponenziale non può **perdere di generalità*, cioè non rischiamo sovra-adattamento al diminuire dell’esponente.

Piuttosto, l’approccio corretto sarebbe quello di regredire anche l’esponente, trattandolo come il coefficiente di una regressione non-lineare, cioè regredire il modello \(y~a + bx^{1/c}\), dove \(c\) compare ovviamente in una relazione non lineare con gli altri coefficienti \(a\) e \(b\).

Proviamo con una regressione ai minimi quadrati inserendo anche l’esponente \(1/c\) come parametro di regressione:

f <- \(x, a, b, c) a + b*x^(1/c)
df.nls <- nls(y ~ f(x, a, b, c),
    data = df,
    start = list(
      a = 1, b=1, c=2
    ),
    trace=T
)
11457.91    (1.99e+01): par = (1 1 2)
1102.181    (6.55e+00): par = (4.304996 7.951322 7.461343)
440.4336    (4.07e+00): par = (-21.86335 34.57408 20.48135)
73.30484    (1.39e+00): par = (-60.8398 73.58786 30.53211)
25.17886    (5.97e-02): par = (-68.75292 81.50195 29.32612)
25.08932    (3.82e-04): par = (-68.82486 81.57389 29.5062)
25.08932    (1.71e-06): par = (-68.83131 81.58033 29.50943)

Si noti come il valore ottimale di \(c\) risulta in effetti moto più alto:

df.nls
Nonlinear regression model
  model: y ~ f(x, a, b, c)
   data: df
     a      b      c 
-68.83  81.58  29.51 
 residual sum-of-squares: 25.09

Number of iterations to convergence: 6 
Achieved convergence tolerance: 1.71e-06

Confrontando la regressione non-lineare con il modello logaritmico si osserva come stanno entrambi nella stessa banda di confidenza, quindi continuiamo a preferire il modello logaritmico:

df %>% 
  add_predictions(df.nls) %>% 
  ggplot(aes(x=x)) +
  geom_point(aes(y=y)) +
  geom_line(aes(y=pred), color="red") + 
  geom_smooth(
    aes(y=y),
    method = "lm",
    formula = y ~ log(x),
    linetype=2
  )