Applicazioni pratiche di machine learning/Previsioni sugli incendi
Caricamento librerie
modifica library(dplyr)
library(ggplot2)
library(h2o)
Parte 1: Dati
modificaIl dataset "forestfires.csv", scaricabile da qui: https://archive.ics.uci.edu/ml/datasets/forest+fires è costituito da 513 righe e 13 variabili e contiene i dati sugli incendi boschivi provenienti dal parco naturale Montesinho, dal Trás-os-Montes regione nord-orientale del Portogallo . Le 13 variabili in questione del dataset sono le seguenti :
- X - Coordinata spaziale dell'asse X all'interno della mappa del parco Montesinho con valori compresi tra 1 e 9
- Y - Coordinata spaziale dell'asse Y all'interno della mappa del parco Montesinho con valori compresi tra 2 e 9
- month - mese dell'anno: da "gen" a "dec"
- day - giorno della settimana: da "mon" a "sun"
- FFMC - Indice FFMC del sistema FWI: 18.7 a 96.20
- DMC - Indice DMC del sistema FWI: da 1.1 a 291.3
- DC - Indice DC del sistema FWI: da 7,9 a 860,6
- ISI - Indice ISI del sistema FWI: da 0,0 a 56,10
- temp - temperatura in gradi Celsius: da 2,2 a 33,30
- RH - umidità relativa in %: da 15,0 a 100
- wind - velocità del vento in km/h: da 0.40 a 9.40
- rain - pioggia esterna in mm/m^2: da 0,0 a 6,4
- area - l'area bruciata della foresta (in ha-ettari): da 0,00 a 1090,84 (questa variabile di output è molto proiettata verso lo 0, quindi potrebbe avere senso modellare con la trasformazione del logaritmo log(X+1)).
Il Forest Fire Weather Index (FWI) è il sistema canadese per la classificazione del pericolo di incendio e comprende sei componenti : FFMC rappresenta il contenuto di umidità dei rifiuti di superficie e influenza l'accensione e la propagazione del fuoco, DMC e DC rappresentano il contenuto di umidità degli strati organici superficiali e profondi che influenzano l'intensità del fuoco. L'ISI è un punteggio correlato alla diffusione della velocità del fuoco.
Caricamento dei dati:
forestfires <- read.csv("forestfires.csv", stringsAsFactors = TRUE)
Ricerca valori mancanti:
colSums(is.na(forestfires))
X Y month day FFMC DMC DC ISI temp 0 0 0 0 0 0 0 0 0 RH wind rain area 0 0 0 0
Ricerca valori anomali o outliers:
forestfires %>%
ggplot(aes(1:517,area)) +
geom_boxplot()
Dal grafico si vede che ci sono almeno 3 outliers relativi ad un'area superiore a 250 ettari, ma comunque non è necessario cancellarli.
forestfires %>%
filter(area>250)
X Y month day FFMC DMC DC ISI temp RH wind rain area 1 6 5 sep sat 92.5 121.1 674.4 8.6 25.1 27 4.0 0 1090.84 2 8 6 aug thu 94.8 222.4 698.6 13.9 27.5 27 4.9 0 746.28 3 7 4 jul mon 89.2 103.9 431.6 6.4 22.6 57 4.9 0 278.53
Parte 2: Domanda di ricerca
modificaSi vuole prevedere l'area bruciata in ettari dagli incendi boschivi. Tale conoscenza è particolarmente utile per la gestione delle risorse antincendio e quindi per definire gli obiettivi prioritari delle autobotti e degli equipaggiamenti di terra (pianificazione delle risorse).
Parte 3: Esplorazione dei dati
modificaCome si vede dal seguente grafico la maggiorparte degli incendi avviene nei mesi di agosto e settembre, mentre a novembre e gennaio non ce ne sono.
forestfires %>%
group_by(month) %>%
summarise(total=sum(area)) %>%
mutate(month=reorder(month,total)) %>%
ggplot(aes(month,total)) +
geom_bar(stat="identity", fill="red") +
coord_flip()
Creo una nuova variabile per poter valutare il tipo di danno:
danno <- rep(NA,nrow(forestfires))
for (i in 1:nrow(forestfires)) {
if (forestfires[i,13]==0)
danno[i]="Nessun danno"
else if (forestfires[i,13]<1)
danno[i]="Danno basso"
else if (forestfires[i,13]<25)
danno[i]="Danno moderato"
else if (forestfires[i,13]<100)
danno[i]="Danno alto"
else danno[i]="Danno molto elevato"
}
forestfires <- cbind(forestfires,danno)
Come si vede dal seguente grafico i danni alti o molto elevati avvengono nei mesi di luglio, agosto e settembre.
forestfires %>%
ggplot(aes(month,area, fill=danno)) +
geom_bar(stat="identity") +
coord_flip()
Parte 4: Modellizzazione e previsione
modificaPer creare un modello previsionale si utilizza la libreria "H2o" che consente con la sua funzione autoML di selezionare automaticamente l'algoritmo di machine learning migliore, cioè quello nel caso specifico che da un più basso valore dell'errore previsionale RMSE.
h2o.init()
Starting H2O JVM and connecting: ....... Connection successful!
R is connected to the H2O cluster: H2O cluster uptime: 5 seconds 344 milliseconds H2O cluster timezone: Europe/Rome H2O data parsing timezone: UTC H2O cluster version: 3.38.0.1 H2O cluster version age: 1 month and 9 days H2O cluster name: H2O_started_from_R_gian_pqa647 H2O cluster total nodes: 1 H2O cluster total memory: 1.91 GB H2O cluster total cores: 2 H2O cluster allowed cores: 2 H2O cluster healthy: TRUE H2O Connection ip: localhost H2O Connection port: 54321 H2O Connection proxy: NA H2O Internal Security: FALSE R Version: R version 4.2.1 (2022-06-23)
Si divide il dataset df in un training set fatto dal 75% delle prime osservazioni e su cui si addestra il modello ed il rimanente 25% costituisce il testing set su cui verrà testato il modello :
n <-nrow(forestfires)
training <- forestfires[1:round(n*0.75),]
testing <- forestfires[round(n*0.75):n,]
Si addestra il modello:
train <- as.h2o(training)
y <- "area"
x <- setdiff(names(train), y)
aml <- h2o.automl(x = x, y = y,
training_frame = train,
max_runtime_secs = 300)
Di seguito in ordine i migliori algoritmi selezionati da H2o, tra cui il primo è GBM_lr_annealing_selection_AutoML_3_20221029_184630_select_model :
lb <- aml@leaderboard
lb
model_id rmse 1 GBM_lr_annealing_selection_AutoML_3_20221029_184630_select_model 57.32336 2 GBM_grid_1_AutoML_3_20221029_184630_model_10 57.54180 3 GBM_grid_1_AutoML_3_20221029_184630_model_3 58.86349 4 GBM_grid_1_AutoML_3_20221029_184630_model_9 59.18868 5 GBM_5_AutoML_3_20221029_184630 59.20226 6 GBM_grid_1_AutoML_3_20221029_184630_model_12 59.28968
L'errore previsionale RMSE dell'algoritmo di Boosting selezionato sul testing set è pari a 73,14 ettari...
model <- aml@leader
test <- as.h2o(testing)
p1 = h2o.predict(model, newdata=test)
p1_df <- as.data.frame(p1)
cat("RMSE_test=",sqrt(mean((p1_df$predict-testing$area)^2)))
RMSE_test= 73.13926
Parte 5: Previsioni in Python
modificaimport pandas as pd
# Caricamento dati
data = pd.read_csv('forestfires.csv')
# Separo la variabile da predire area dai predictors
y = data.area
X = data.drop(['area'], axis=1)
# Divido i dati in training e validation subsets, il primo con il 75% delle partite e il secondo con il 25%
X_train_full=X.iloc[1:round(0.75*X.shape[0])+1,]
X_valid_full=X.iloc[round(0.75*X.shape[0])+1:,]
y_train=y.iloc[1:round(0.75*X.shape[0])+1]
y_valid=y.iloc[round(0.75*X.shape[0])+1:]
# Seleziono le colonne con valori categoriali in numero minore di 12
categorical_cols = [cname for cname in X_train_full.columns if X_train_full[cname].nunique() < 12 and X_train_full[cname].dtype == "object"]
# Seleziono le colonne numeriche
numerical_cols = [cname for cname in X_train_full.columns if X_train_full[cname].dtype in ['int64', 'float64']]
# Mantengo soltanto le colonne selezionate
my_cols = categorical_cols + numerical_cols
X_train = X_train_full[my_cols].copy()
X_valid = X_valid_full[my_cols].copy()
Le colonne numeriche e categoriche nel dataset che hanno valori mancanti NA vengono trasformate con valori approssimativi nei NA tramite SimpleImputer, mentre le variabili categoriche vengono trasformate in numeriche tramite OneHotEncoder :
from sklearn.compose import ColumnTransformer
from sklearn.pipeline import Pipeline
from sklearn.impute import SimpleImputer
from sklearn.preprocessing import OneHotEncoder
# Preprocessing per i dati numerici
numerical_transformer = SimpleImputer(strategy='constant')
# Preprocessing per i dati categoriali
categorical_transformer = Pipeline(steps=[
('imputer', SimpleImputer(strategy='most_frequent')),
('onehot', OneHotEncoder(handle_unknown='ignore'))
])
preprocessor = ColumnTransformer(
transformers=[
('num', numerical_transformer, numerical_cols),
('cat', categorical_transformer, categorical_cols)
])
Per la previsione utilizzo l'algoritmo RandomForest e addestro il modello sul training set dopo avere preprocessato i dati con il preprocessor precedentemente creato. Infine prevedo i dati sul validation set.
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(random_state=0)
my_pipeline = Pipeline(steps=[('preprocessor', preprocessor),
('model', model)
])
# addestro il modello
my_pipeline.fit(X_train, y_train)
# ottengo le previsioni sul validation set
preds_valid = my_pipeline.predict(X_valid)
Sul validation set ottengo un errore previsionale RMSE di 79.08 ettari... in R ho ottenuto il risultato migliore pari a 73,14 ettari...
from sklearn.metrics import mean_squared_error
print("RMSE_valid=",mean_squared_error(y_valid,preds_valid, squared=False))
RMSE_valid= 79.07977319002202