venerdì 23 aprile 2021

L’oracolo del cancro al seno

Il tasso di morte per carcinoma al seno nella popolazione femminile e’ molto alto (28% dei tumori femminili) e in Italia e’ la prima causa di morte nella fascia tra i 35 e 50 anni. Secondo la World Health Organization, il cancro al seno colpisce piu’ di 1.5 milioni di donne al mondo. La prima testimonianza di un cancro al seno e’ quella di una donna egizia del 1600 a.C. morta a causa di questa malattia. Contrariamente a quanto si pensa oggi il cancro non e’ una malattia del nostro secolo.  

 

 

Il rischio per una donna di ammalarsi di tumore al seno e’ del 10-12% e varia in base all’eta’:

·        fino a 49 anni à 2.4% (1 donna su 42)

·        tra i 50 e i 69 anni à 5.5% (1 donna su 18)

·        tra i 70 e gli 84 anni à 4.7% (1 donna su 21)

Data l’elevata incidenza ed il tasso di mortalita’ negli ultimi anni sono stati adottati degli screening ben precisi e definiti per identificare precocemente la patologia. Le evidenze scientifiche dimostrano come la diagnosi precoce sia l’azione piu’ efficace nella lotta contro il tumore della mammella. 

Le tecniche diagnostiche piu’ indicate per questa malattia sono l’ecografia, la mammografia, agoaspirato, TAC, PET e RNM. In caso di nodulo al seno, la diagnosi e’ quindi un primo fondamentale passaggio. Una tardiva o non corretta diagnosi della patologia potrebbe infatti portare, nei casi piu’ gravi al decesso della donna. Nonostante i progressi scientifici fatti in ambito medico purtroppo in assoluto non e’ possibile escludere che si possano verificare degli errori medico-sanitari anche durante la fase diagnostica. Questi errori per esempio possono derivare da difetti del macchinario, procedure non seguite scrupolosamente, scambio di analisi, esami mal letti o interpretati, mancato approfondimento con ulteriori visite o esami integrativi. Per questo motivo con l’avvento dell’intelligenza artificiale si e’ pensato di affiancare al medico degli oracoli digitali per riconoscere alcune forme di tumore. Numerosi ricercatori hanno pensato di applicare le tecniche proprie del machine learning per identificare il carcinoma al seno e non solo.

       

Immagine ecografica di un cancro al seno         Immagine mammografica di cancro al seno

 Negli ultimi anni, una grande quantita’ di dati relativi al tumore del seno sono stati collezionati e resi disponibili sul web in diveri siti (vedi come esempio UCI). Tra questi database c’e’ quello dell’Universita’ del Wisconsin utilizzato in questo post. E’ un classico database per l’applicazione di algoritmi di classificazione binaria per il machine learning. E’ composto da 569 esempi o istanze, 31 predittori o features e una risposta che in questo caso e’ una variabile binaria con due valori: benigno B e maligno M. Ci sono 357 casi di tumore benigno e 212 di tumore maligno. Le due classi sono abbastanza bilanciate con un rapporto del 60%.  

 

 Le variabili x (i predittori o features) sono state registrate utilizzando delle immagini digitali della massa tumorale catturata con la tecnica dell’ago aspirato. In effetti le variabili sono 10, ma diventano 30 in quanto e’ stato preso il valore medio, la standard deviation e il valore massimo. Si tratta di 10 misure tutte relative ai nuclei delle cellule tumorali (vedi immagine di seguito).  




Le features estratte dall’immagine sono: 

1.     Raggio di un nucleo (media di tutti i raggi, standard deviation e valore massimo)

2.     Perimetro di ogni nucleo mediato

3.     Area del nucleo (valore medio, std e valore max)

4.  Compattezza ottenuta combinando opportunamente l’area e il perimetro (P^2/A)

5.     Smoothness ottenuta facendo la differenza tra il raggio in un certo punto e il valore          medio dei primi vicini

6.     Concavita’ 

7. Punti di concavita’. Diversamente dalla concavita’ questo parametro determina quanti  punti di concavita’ ci sono.

8.     Simmetria misurata tracciando l’asse maggiore passante per il centro e poi facendo la differenza della lunghezza tra linee perpendicolari all’asse maggiore in entrambe le direzioni.

  9.  Dimensione frattale ottenuta approssimando la linea chiusa con la “linea di costa”            descritta da Mandelbrot

10. Texture ottenuta trovando la varianza dell’intensita’ della scala dei grigi nei pixels


Per questo studio e’ stato utilizzato JMP Pro 16 della SAS che contiene una serie di algoritmi di machine learning utilizzabili per scopi di classificazione supervisionata. Tra i modelli disponibili abbiamo: Decision tree, Random Forest, Boosted tree, Logical regression, Support Vector Machine, Naive Bayes, XGboost, kNN, Neural Network (per maggiori dettagli sui diversi  modelli di machine learning fare riferimeno al post “Alla scoperta delle leggi della fisica con il machine learning” su questo blog). Prima di iniziare, controlliamo le distribuzioni delle nostre features che risultano tutte leggermente skewed sulla destra (mostrano cioe’ una leggera coda verso gli alti valori). Non ci sono valori mancanti nel database e solo alcuni outliers.

 


Il training per questo dataset e’ stato fatto su 8 classificatori usando l’opzione 5-fold cross validation. Per stabilire il grado di accuratezza e l’efficacia di qualsiasi modello di machine learning e’ necessario eseguire una o piu’ valutazioni sugli errori che si ottengono con ogni previsione. In genere, dopo il training viene effettuata una stima dell’errore per il modello, meglio nota come valutazione dei residui. Si tratta della stima numerica della differenza tra la risposta prevista e quella originale, chiamata anche errore di addestramento (training error). Tuttavia questo non basta in quanto non da’ un idea di quanto il nostro modello riuscira’ a generalizzare in caso di nuovi dati. Per questo motivo viene utilizzata la valutazione incrociata. Essa consiste nella suddivisione dell’insieme di dati in k parti (5 nel nostro caso) di uguale numerosita’ e a ogni passo la k-esima parte dei dati viene usata come convalida, mentre la restante parte costituisce l’insieme di training (addestramento). In questo modo si allena il modello per ognuna delle k parti evitando problemi di overfitting (sovradattamento) ma anche di campionamento asimmetrico (distorsione) tipico della suddivisione dei dati in due sole parti. 

Per stabilire quanto un modello e’ bravo nella classificazione esistono diverse metriche. Quelle piu’ usate sono legate al concetto di matrice di confusione. La nostra risposta puo’ avere solo due possibili valori M e B. Quindi accoppiando le risposte del modello e quelle effettive possiamo avere la seguente matrice: 

 

FN rappresenta il numero di casi in cui il modello ha previsto Malignant e invece nella realta’ il tumore e’ Benign (si parla di falsi negativi). Viceversa per il caso FP (falsi positivi) che rappresenta i casi missing. Il rapporto tra la somma FN+FP e il numero totale di esempi TP+TN+FP+FN (dove TP e TN indica i totali positivi e negativi rispettivamente) e’ chiamato misclassification error ME da cui si ricava immediatamente l’accuratezza data da:

A=1-ME=1-((FN+FP)/(FN+FP+TN+TP))

Osserviamo che come ci si aspetta se l’errore e’ nullo allora l’accuratezza sara’ massima e cioe’ 1. Nel caso invece in cui l’errore e’ massimo e cioe’ 1 allora l’accuratezza del modello e’ pari a zero. Un’altra metrica e’ la Precisione P del modello che indica l’abilita’ di un classificatore di non etichettare un’istanza positiva che in realta’ e’ negativa.

P=TP/(TP+FN)

Altre metriche utilizzate da JMP Pro per i classificatori sono l’R2, il RASE e l’AUC. Trattandosi di un classificatore che lavora su dati categorici per la risposta ci si chiede come e’ possibile che JMP restituisca il valore del R2 che tutti associamo alle regressioni numeriche? E la stessa domanda nasce anche quando le X sono categoriche. Un algoritmo dopo tutto e’ una sequenza di operazioni matematiche. Ma come si fa a sommare due variabili nominali quali Benigno e Maligno?


Nell‘ambito del Machine learning si usa il cosiddetto algoritmo “one hot encoder”.  Nell’esempio dei colori della tabella qui sotto, la variabile ha 3 stati (red, blue e green) e quindi si creano altrettante variabili dummy con valori di 0 e 1.    

In questo modo e’ possibile calcolare il valore del R2 essendo adesso le variabili tutte numeriche. Allo stesso modo e’ possibile calcolare il cosiddetto RASE, cioe’ la radice quadrata del valore medio dell’errore (y-ya) al quadrato (ya e’ il valore effettivo e y quello previsto dal modello). L’AUC indica invece l’area al di sotto della curva ROC. Vediamo di cosa si tratta. In qualsiasi problema di classificazione, il modello genera per ogni record un valore di probabilita’ di appartenenza ad uno dei 2 gruppi. Per decidere quindi se un record appartiene ad una classe piu’ che ad un’altra, e’ importante stabilire un valore di soglia per la probabilita’. Ovviamente al variare di tale soglia cambiera’ l’output del modello e la relativa matrice di confusione. Per stabilire se un modello sta lavorando bene si puo’ utilizzare la curva di ROC. Sull’asse della x viene riportato il numero di falsi positivi (FP) o 1-SPE dove SPE=TN/(TN+FP)  e’ la cosiddetta Specificita’, mentre sull’asse delle y viene riportata la quantita’  TP/(TP+FN) che in modo indiretto e’ legata al numero di falsi negativi (quando FN e’ zero infatti questo rapporto vale 1). Il migliore modello di classificazione sara’ quello che ottiene FP=0 e FN=0, che nel piano della curva ROC corrisponde al punto x=0 e y=1. I due plot ROC seguenti si riferiscono al caso della rete neurale per il set di Training e di Validation. I colori rosso e blu indicano i due valori della classe di risposta: Benign e Malignant. L’AUC come indica l’acronimo e’ il valore dell’area al di sotto della curva ROC. Piu’ e’ prossimo ad 1 e meglio e’ la prestazione del modello.    




 

L’accuratezza di tutti i modelli testati con il comando Model Screening di JMP e’ ben al di sopra del 90%, il che significa che il massimo errore per mis-classification e’ minore del 10%. L’algoritmo del Neural Boosted addirittura mostra un errore del solo 1%. Non c’e’ che dire, l’oracolo digitale sembra fare delle ottime previsioni. Attenzione. Anche se nella lingua inglese non c’e’ una grossa differenza tra prediction e prevision (o forecast) questo non e’ vero in italiano. Nella nostra lingua la predizione e’ quella fatta da un mago che cerca di capire cosa accadra’ nel futuro aiutandosi con strumenti che di scientifico non hanno praticamente nulla (e quindi fallira’ miseramente). Al contrario la previsione e’ quella che cerca di capire cosa puo’ succedere utilizzando metodi scientifici (matematici, fisici, statistici) come la previsione del ritorno della cometa di Halley o la probabilita’ che il costo delle azioni Apple crescera’ nei prossimi giorni.   

Ritorniamo ai modelli testati. Il migliore e’ la rete Neurale Boosted. Ma cosa significa boosted? E’ una classe di modelli nati nel 1988 con l’idea che mettendo insieme piu’ modelli di apprendimento deboli si possa creare un modello piu’ forte (della serie che l’unione fa la forza). Si tratta di un modello iterativo (lavora in seriale) che stabilisce come collegare tra loro un insieme di weak learner per crearne uno strong. Anche se l’accuratezza raggiunta da questo modello e’ molto alta, il fatto che ci siano alcuni casi in cui abbiamo predetto che il cancro e’ benigno mentre invece e’ maligno non ci piace affatto, visto che si ha a che fare con le vite delle persone. Meglio caso mai  avere un Falso negativo (diciamo che e’ maligno ma in realta’ e’ benigno) che oltre alla paura non fara’ altri danni alla persona sottoposta alla diagnosi. C’e’ da dire comunque che nel Machine learning e’ possibile provare a penalizzare gli esempi che ricadono nella casella FN rispetto a quella FP. In JMP Pro questo puo’ essere fatto direttamente dal Model Screening utilizzando l’opzione Decision Thresholds. Questa permette di esplorare la soglia dei modelli per la classificazione binaria. C’e’ un report per ogni modello specificato dal metodo di validazione. 

 

Ogni report contiene un grafico della distribuzione delle probabilita’ previste, delle carte a barra per le diverse classificazioni e la matrice di confusione. Nelle barre la porzione in rosso e’ proporzionale all’errore di mis-classification. Spostando la linea nera al centro del grafico delle distribuzione si puo’ cambiare la soglia e cercare di diminuire il numero di falsi positivi rispetto a quelli negativi. Con la scelta operata nel nostro caso si e’ potuto ottenere un azzeramento dei Falsi positivi per le NN Boosted raggiungendo un’accuratezza del 100%.

 

Anche se in JMP le opzioni che vado a descrivere adesso vengono implementate automaticamente, in generale usando linguaggi come Python o R e le loro librerie, conviene prima di passare al training/test del modello di normalizzare le variabili X per esempio facendo in modo che tutti i predittori siano nel range 0-1 e che questi vengano trasformati con una funzione tipo logaritmo per cercare di eliminare la skewness  della distribuzione.  In definitiva i 5 steps piu’ importanti in ogni attivita’ di Machine learning sono:

1.   Data collection: si tratta dello step dove viene raccolto il materiale da dare in pasto agli algoritmi per trasformarlo in conoscenza utilizzabile. Nella maggior parte dei casi i dati devono essere combinati in una singola sorgente come un file testo, csv o excel. 

2. Data exploration/preparation: la qualita’ di qualsiasi progetto di machine learning dipende dalla qualita’ dei dati in ingresso. Quindi ogni qualvolta si parte col costruire un modello si devono pulire i dati dal rumore, eliminare quelli non necessari, e riempire le celle vuote del database (missing value). 

3. Model training: una volta che i dati sono stati prepararti si divide il set in training/validation/test e si fa partire la ricerca

4.   Model evaluation: poiche’ ogni machine learning tende ad essere biasato e’ importante valutare le prestazioni dell’algoritmo in termini di generalizzazione. Per fare questo si utilizzano diversi tipi di metriche a secondo che si tratta di un problema di regressione o di classificazione. 

5.  Model improvement: nel caso in cui siano necessarie prestazioni migliori si puo’ pensare di utilizzare delle strategie avanzate. Qualche volta basta cambiare il modello, o costruire dei nuovi predittori (feature engineering). Altre volte in caso di underfitting del sistema semplicemente raccogliere piu’ dati.             

 


sabato 3 ottobre 2020

Quando un diverso numero di neutroni può aiutare l’archeologo

 Nell’ultimo decennio, in tutti i campi delle Scienze Ambientali e dei Beni Culturali, è emersa una forte esigenza di metodologie innovative da affiancare a quelle convenzionali per una migliore comprensione dei processi ambientali e antropogenici. Per questo motivo le metodologie isotopiche hanno assunto un ruolo sempre più importante in un numero crescente di procedure di analisi e di controllo sia in campo industriale, che ambientale, biomedico e archeologico. Si tratta, in generale, di metodi di misura sviluppati nella ricerca di base in fisica e trasferiti ai laboratori in vari altri campi. Tra questi metodi la spettrometria di massa di isotopi stabili degli elementi più abbondanti in natura, costituisce un potente mezzo di indagine, largamente utilizzato nelle scienze applicate alla diagnostica dei beni culturali ed ambientali. L’ambiente è, per sua natura, in continua evoluzione ma l’uomo ne ha accelerato il passo in una direzione incognita. Ecco, allora, che la determinazione della composizione isotopica di elementi quali il Carbonio, l’Ossigeno e l’Azoto nell’atmosfera, nel suolo, nella vegetazione e nelle acque, si è dimostrata un indicatore molto sensibile delle piccole variazioni (spesso dovute proprio ad attività antropiche) dei flussi di questi elementi e per studiare meccanismi ancora poco chiari. Gli isotopi sono atomi di uno stesso elemento che posseggono un diverso numero di neutroni. Gli isotopi vengono definiti stabili quando il numero di protoni e neutroni non cambia nel tempo. Poiché i neutroni all’interno dell’atomo sono particelle elettricamente neutre, due isotopi dello stesso elemento mantengono la stessa carica ma hanno un diverso numero di massa (A). Il numero di massa è espresso dalla somma del numero di neutroni (N) e di protoni (Z, particelle cariche positivamente) ed è riportato come apice a sinistra del simbolo chimico. Qui di seguito un diagramma che mostra come la stabilità nucleare è più un’eccezione che una regola con solo 270 nuclidi stabili su 2500 conosciuti.



A causa della diversa massa, gli isotopi di un elemento hanno forze di legame diverse e di conseguenza possono avere comportamenti chimico-fisici diversi. La composizione in isotopi stabili di elementi di bassa massa atomica è riportata in termini di valori relativi ad uno standard di composizione isotopica nota. Questa notazione, il “delta per mille” (δ ‰) è calcolata secondo l’equazione:

dove R indica il rapporto tra l’isotopo più pesante e quello più leggero. In pratica, il valore δ esprime di quante parti per mille il rapporto isotopico si discosta dallo stesso rapporto in uno standard di riferimento con composizione isotopica nota. Il motivo dell’utilizzo di tale notazione risiede nel fatto che la differenza tra due isotopi di uno stesso elemento è molto piccola per cui utilizzare il rapporto R da solo implicherebbe il maneggiamento di numeri con molte cifre decimali e complicherebbe i calcoli. L’uso della notazione δ espande piccole differenze e le rende più facilmente visibili. Inoltre, l’uso di un rapporto di rapporti permette di confrontare elementi diversi che sarebbero molto più difficili da confrontare in valori assoluti.

Un valore δ positivo indica che il rapporto isotopico del campione considerato è più alto di quello dello standard, cioè che il campione è arricchito nell’isotopo pesante rispetto allo standard. Viceversa, un valore δ negativo indica che il rapporto isotopico del campione considerato è più basso di quello dello standard, e quindi che il campione è impoverito nell’isotopo pesante rispetto allo standard. Diversi termini vengono comunemente usati per comparare i valori δ: arricchito/impoverito in isotopi pesanti; più positivo/più negativo; maggiore/minore; più pesante/più leggero.


Esistono diversi standard internazionali di riferimento per vari isotopi (vedi tabella di seguito). Per quanto riguarda gli isotopi stabili dell’acqua (2H e 18O) per esempio, comunemente, lo standard di riferimento è rappresentato dallo VSMOW (Vienna Standard Mean Ocean Water), rimpiazzato nel 2006 dallo VSMOW2 che, per definizione, ha un valore di δ di 0‰.


Le proporzioni relative dei diversi isotopi non sono sempre le stesse a causa del cosiddetto frazionamento isotopico. Questo fenomeno è dovuto a differenze nella velocità di reazione degli isotopi leggeri e pesanti dello stesso elemento durante le reazioni termodinamiche. Nell’acqua, per esempio, il frazionamento isotopico avviene principalmente durante i cambiamenti di fase (evaporazione, condensazione, fusione). In questi casi, la velocità di diffusione varia a seconda della forza dei legami chimici tra gli isotopi pesanti e leggeri. I legami molecolari tra gli isotopi leggeri vengono spezzati più facilmente di quelli tra isotopi pesanti. La forma isotopica pesante dell’acqua (cioè quella contenente 2H e 18O, ossia 2H18O) richiede più energia per spezzare i legami chimici rispetto all’acqua che contiene isotopi più leggeri (1H16O), e di conseguenza la reazione avviene più lentamente. Nel caso del cambiamento di stato da liquido a gassoso, per esempio, le molecole di acqua che contengono 16O evaporano più facilmente rispetto a quelle che contengono 18O, a causa della differenza di massa tra i due isotopi e della facilità di rottura dei legami idrogeno che mantengono l’acqua in forma liquida. Inoltre, maggiore è la differenza di massa tra due isotopi, maggiore è il frazionamento isotopico. Esistono due meccanismi principali che determinano il frazionamento isotopico: frazionamento all’equilibrio e frazionamento non all’equilibrio (frazionamento cinetico). Il primo comporta la redistribuzione degli isotopi di un elemento tra vari composti, cosicché la velocità di reazione diretta e inversa di un certo isotopo è identica. Questo implica che il rapporto tra isotopi differenti in ogni composto o in ogni fase della reazione per una certa temperatura rimanga costante ma non implica necessariamente che la loro concentrazione sia uguale. Durante le reazioni all’equilibrio, l’isotopo più pesante si accumula nel composto col maggiore stato di ossidazione o nel composto più “denso”. Un tipico esempio di processo all’equilibrio è rappresentato dalla condensazione del vapore acqueo. In tale processo, gli isotopi pesanti della molecola d’acqua (18O and 2H) si concentrano nella fase liquida (che così viene definita arricchita in isotopi pesanti) mentre gli isotopi leggeri (16O and 1H) rimangono concentrati nella fase gassosa. In generale, il frazionamento è fortemente dipendente dalla temperatura: più alta è la temperatura, minore è la differenza tra la composizione isotopica di due specie all’equilibrio. Il frazionamento non all’equilibrio è tipico di sistemi aperti, dove il prodotto viene continuamente allontanato dal sistema, o nelle reazioni irreversibili. In questi casi, non è possibile raggiungere l’equilibrio fra reagenti e prodotti. La velocità di reazione è dipendente dal rapporto delle masse isotopiche e dalle loro energie vibrazionali. Come nel caso del frazionamento all’equilibrio, i legami tra gli isotopi leggeri vengono spezzati più facilmente di quelli tra isotopi pesanti. Questo comporta un accumulo di isotopi pesanti nel reagente mentre nelle reazioni reversibili all’equilibrio i prodotti della reazione possono essere isotopicamente sia più leggeri che più pesanti dei reagenti originari. Un tipico esempio di frazionamento isotopico cinetico è rappresentato dal processo di evaporazione dell’acqua dalla superficie degli oceani, in cui il vapore acqueo prodotto dall’evaporazione viene continuamente rimosso dai processi di turbolenza dell’atmosfera cosicché la reazione non può mai raggiungere l’equilibrio. Il vapore prodotto dalla reazione di evaporazione risulta impoverito in isotopi pesanti con una deviazione del rapporto di frazionamento rispetto a quello che avverrebbe all’equilibrio.

Gli isotopi stabili vengono utilizzati per risolvere diversi problemi. Una delle applicazioni più comuni è la geotermometria, cioè la misura di temperature che si sviluppano nel corso dei fenomeni geologici. Il frazionamento isotopico è inversamente proporzionale alla temperatura: il frazionamento è grande alle basse temperature e piccolo alle alte temperature. La geotermometria è inoltre un’analisi geochimica che consente di determinare, indirettamente e attraverso la composizione dei diversi minerali costituenti una roccia, i valori della temperatura alla quale la roccia si è andata formando. Questo metodo si applica se esistono condizioni di equilibrio termodinamico tra i componenti.

Un’altra applicazione degli isotopi stabili è l’identificazione di processi. Le piante, per esempio, che per fotosintesi producono catene di idrocarburi C4 (cioè idrocarburi con 4 atomi di Carbonio) frazionano il carbonio diversamente dalle piante che producono catene di idrocarburi C3. Questo frazionamento è mantenuto nella catena alimentare e quindi la misura di isotopi stabili nelle ossa di mammiferi fossili ci da’ informazioni sulla loro dieta. Le tecniche isotopiche giocano un ruolo importante anche nelle scienze archeologiche, dove sono state utilizzate per stabilire le diete del passato e la migrazione dei popoli, la provenienza di marmi usati in architettura e scultura, la provenienza di antiche ceramiche e metalli, la ricostruzione di paleoclimi e paleoambienti e l’utilizzo di oggetti di terracotta grazie all’analisi dei residui di cibo presenti sulla loro superficie.

Tra gli isotopi usati dagli archeologi quelli più importanti sono:

  • carbonio: 13C (o C-13) e 12C (o C-12). Il rapporto tra i due (13C/12C) viene indicato come δ13C
  • azoto: 15N (o N-15) e 14N (o N-14). Di nuovo il rapporto (15N/14N) è dato da δ15N.
  • ossigeno: 18O (O-18) e 16O (O-16), con il loro rapporto (18O/16O) indicato come δ18O.
  • stronzio: 87Sr (Sr-87) e 86Sr (Sr-86). Il rapporto è dato da 87Sr/86Sr.

Per il carbonio le due principali riserve in natura sono rappresentate dal carbonio organico e dai carbonati. Queste sono caratterizzate da differenti impronte isotopiche a causa dei diversi processi in cui sono coinvolte.

Il carbonio inorganico (carbonatico) interviene negli equilibri di scambio tra anidride carbonica atmosferica - ione bicarbonato disciolto - carbonato solido; queste reazioni di scambio portano ad un arricchimento dell’isotopo pesante nella forma carbonatica solida (δ13C pari a 0‰). Viceversa, le reazioni cinetiche in cui risulta principalmente coinvolto il carbonio organico, attraverso i processi fotosintetici, determinano una concentrazione dell’isotopo più leggero nel materiale organico prodotto (δ13C pari a circa −25‰). Il frazionamento del carbonio organico è principalmente legato al tipo di pianta che opera la fotosintesi. Le piante terrestri, classificate come C3 e C4, seguono due vie metaboliche differenti. Entrambe generano sostanza organica caratterizzata da valori di δ13C più negativi rispetto a quello dell’anidride carbonica (~-7‰), in quanto durante la fotosintesi la sostanza organica prodotta accumula l’isotopo leggero rispetto a quello pesante. Le piante C3, caratteristiche dei climi temperati, producono un composto a tre atomi di carbonio (ciclo di Calvin) con un valore medio di δ13C pari a -26,5‰. Le piante C4 generano un composto a 4 atomi di carbonio (ciclo di Hatch-Slack) caratterizzato da un valore di δ13C intorno a -12,5‰. Partendo da questo tipo di conoscenze il frazionamento del carbonio di origine organica è stato indagato per molteplici applicazioni, tra le quali lo studio delle reti trofiche. La composizione chimica dei tessuti animali è connessa alle risorse alimentari che essi assimilano e pertanto i tessuti riflettono la composizione isotopica della dieta. L’arricchimento tra i produttori e i consumatori primari (erbivori) è stato stimato intorno al +5‰ mentre andando ai livelli successivi della catena trofica l’arricchimento risulta meno marcato (+1‰). Il valore isotopico rilevato nei tessuti di un organismo, quindi, può essere utilizzato come indicatore della sua posizione trofica, ma dato che la variazione dei valori di δ13C associata ai passaggi trofici risulta relativamente modesta, questo dato è usato principalmente per rintracciare la fonte primaria di carbonio assunta. Attraverso l’analisi degli isotopi stabili del carbonio è possibile anche differenziare le reti trofiche terrestri da quelle marine. Il carbonio “marino” deriva infatti dal carbonio inorganico disciolto (bicarbonato disciolto) caratterizzato da un valore isotopico pari a circa 0‰, e quindi più alto di quello dell’anidride carbonica atmosferica, che è pari a circa -7‰. Tale differenza è mantenuta ad ogni livello trofico sia in ambiente marino che terrestre.

Passiamo adesso all’azoto. Esso è il principale elemento dell’atmosfera terrestre (circa il 78%), ma nonostante ciò la maggior parte degli organismi non sono in grado di usarlo in forma gassosa. L’azoto atmosferico deve essere quindi convertito in forme utilizzabili e ciò avviene naturalmente attraverso una serie di reazioni chimiche mediate dai microorganismi azotofissatori che vivono sia nel suolo che nelle acque dolci e salate (Clostridium, Azotobacter, Rhizobium leguminosarum, attinomiceti), e che producono un frazionamento isotopico dell’azoto differenziando i valori di δ15N. La prima fase della fissazione è quella che vede la formazione di azoto ammonico (δ15N pari a circa 1‰) da quello molecolare atmosferico.


Le reazioni che generano ulteriori trasformazioni dell’azoto nel suolo e nelle acque sono la mineralizzazione, la volatilizzazione, la nitrificazione e la denitrificazione, anch’esse per la maggior parte mediate da microrganismi. Ad esempio, la volatilizzazione è legata alla perdita di ammoniaca dal suolo all’atmosfera. Tale processo è caratterizzato da un elevato frazionamento che produce ammoniaca impoverita in δ15N, lasciando lo ione ammonio residuo presente nel suolo arricchito dell’isotopo pesante. L’intero processo di trasformazione dell’azoto in nitrato coinvolge diversi passaggi di reazione, ognuno dei quali produce un arricchimento nel substrato azotato residuo che può arrivare al 30‰. La denitrificazione (cioè l’utilizzo del nitrato al posto dell’ossigeno quale substrato per l’ossidazione della materia organica) è un processo mediato da batteri, in grado di provocare un elevato frazionamento isotopico del nitrato residuo. In relazione anche alle condizioni ambientali, l’attività batterica discrimina le forme isotopiche più “leggere”, determinando nel nitrato residuo un arricchimento dell’isotopo pesante rispetto all’azoto molecolare prodotto. Nella figura sottostante vengono riportati i principali processi che interessano il ciclo dell’azoto in ambiente marino.


La conoscenza dell’intero ciclo dell’azoto in termini di distribuzione isotopica, in associazione a quello del carbonio, è utile sia per comprendere le caratteristiche della catena trofica che le relazioni tra alcuni tipi di pressioni antropiche sull’ambiente e gli impatti sugli ecosistemi.

Ora passiamo all’ossigeno. Durante i fenomeni di evaporazione e precipitazione si verifica il frazionamento isotopico dell’ossigeno (e dell’idrogeno). L’entità del frazionamento dipende dalla temperatura e da altri fattori climatici e geografici come ad esempio latitudine, altitudine, stagionalità e continentalità. L’acqua tende ad evaporare con una reazione di equilibrio regolata dalla temperatura. La fase vapore è caratterizzata da un arricchimento in 16O (δ18O < 0) e in 1H, essendo le molecole di 1H16O più leggere e quindi più favorite all’evaporazione. La fase liquida sarà viceversa più ricca in 2H16O, 1H18O e 2H18O, molecole più pesanti (δ18O e δ2H >0).

Le precipitazioni si impoveriscono dell’isotopo pesante all’aumentare della latitudine e dell’altitudine. In una stessa regione, le precipitazioni relative ai mesi freddi sono caratterizzate da composizioni isotopiche negative mentre durante i mesi caldi risultano arricchite in isotopi pesanti (δ più positivi per l’effetto stagionalità). Le precipitazioni risultano infine più arricchite dell’isotopo pesante spostandosi da regioni costiere all’entroterra. Anche nel ciclo vegetativo delle piante, i processi di adsorbimento dell’acqua e di evapotraspirazione determinano un arricchimento degli isotopi pesanti (2H e 18O), in dipendenza dalla specie vegetale e dall’acqua “isotopicamente” diversa che le varie specie vegetali hanno a disposizione per la fotosintesi.

Il quarto elemento anch’esso molto utilizzato dagli archeologi e’ lo stronzio. Esso ha quattro isotopi naturali: 84Sr, 86Sr, 87Sr e 88Sr tutti quanti stabili. Le abbondanze isotopiche sono variabili a causa della formazione del radiogenico 87Sr attraverso il decadimento del 87Rb. Per questo motivo l’esatta composizione isotopica delle rocce dipende dal rapporto iniziale di Rb/Sr durante la fase di cristallizzazione, dall’età della roccia e da quanto 87Sr si è formato, dall’interazione con fluidi o attività metamorfica. Questa quantità rimane inalterata per tutta la vita della roccia, anche quando questa si sgretola andando a formare il terreno dove crescono piante che a loro volta entrano nel ciclo di vita degli erbivori.

Ma come facciamo a misurare quantità cosi esigue di elementi isotopici? Quali apparecchiature possono aiutare l’archeologo o il geologo ad effettuare studi sui rapporti isotopici? La risposta è la spettrometria di massa. Si tratta di una tecnica che viene utilizzata per separare molecole cariche, in base alla loro massa, o più correttamente in base al loro rapporto massa/carica grazie ad un campo magnetico. Quindi è una tecnica in grado di distinguere isotopi dello stesso elemento e di calcolarne i rapporti isotopici. Questi ultimi nel caso degli elementi leggeri (H, O, C, N, S) vengono misurati con la tecnica IRMS (isotope ratio mass spectrometry) mediante trasformazione in gas puri. Per altri elementi, quali ad esempio B, Fe e Sr viene utilizzata invece la tecnica ICP-MS (Inductively coupled plasma-mass spectrometry) che prevede la trasformazione del campione da analizzare in atomi e ioni.


Passiamo adesso a qualche esempio di applicazione nell’ambito archeologico. Una delle attività di ricerca in cui si fa largo uso degli isotopi stabili è quella dello studio dei resti fossili umani ed animali per risalire all’ecologia e alla dieta degli individui investigati e ricostruire le condizioni paleoclimatiche e paleoambientali dei loro habitat. A questo scopo vengono misurati i rapporti isotopici di C, N, S, O e Sr nelle frazioni minerali (carbonato idrossiapatite) e organiche (collagene) di ossa e denti. Questi risultano per vari motivi indicatori delle condizioni ambientali in cui gli organismi sono vissuti e possono fornire informazioni di tipo paleoidrologico, paleoclimatico, sulle diete, sulla disponibilità e qualità del cibo, ecc. Alcuni di questi parametri possono essere utilizzati inoltre per identificare la mobilità e le migrazioni degli individui studiati. Lo scheletro dei mammiferi è costituito da una frazione minerale di carbonato idrossiapatite (bioapatite) e da una struttura di sostegno costituita da fibre di collagene nelle quali sono dispersi i minerali. È stato osservato che entrambi i tessuti sono in grado di “registrare” nella loro composizione isotopica le condizioni dell’habitat. Il δ13C del collagene e della bioapatite deriva principalmente dalla dieta, a sua volta dipendente dalla distribuzione e dalla tipologia delle piante. Il δ15N risulta anche correlato direttamente con l’aridità dell’habitat. Il δ34S risulta invece dipendente anche dal contesto geologico, dato che subisce scarsi frazionamenti durante i vari passaggi della catena trofica, per cui può essere utilizzato per ricostruire la mobilità degli individui, esattamente come il 87Sr/86Sr. Il δ18O della porzione fosfatica (δ18Op) e carbonatica (δ18Oc) della bioapatite è dipendente da tutte le fonti di ossigeno (cibo, acqua, O2 respiratorio, ecc.) ma per il bilancio di massa l’acqua ingerita, essendo in maggiore quantità rispetto alle altre componenti, risulta generalmente responsabile della composizione isotopica del minerale dello scheletro. Prima di eseguire le misure con lo spettrometro di massa vengono utilizzate specifiche tecniche di estrazione per separare i tessuti di interesse. Un secondo esempio è quello dell’utilizzo degli isotopi stabili per determinare la provenienza dei marmi utilizzati nel passato come per esempio quelli mostrati qui di seguito.

La tecnica è stata utilizzata per la prima volta nel 1972, quando Harmon e Valerie Craig proposero l’uso degli isotopi stabili del C e del O per distinguere 170 diversi marmi da quattro cave greche.

Grazie a questi studi gli archeologi sono riusciti ad individuare le principali cave presenti nel mediterraneo al tempo dei romani (vedi mappa a di seguito) e nell’antica Grecia (vedi mappa b di seguito). Chiaramente la tecnica degli isotopi stabili è una di alcune tecniche a disposizione dell’archeologo per stabilire la provenienza del marmo. Solo grazie ad un’analisi comparata del campione è possibile stabilire con certezza la sua cava di provenienza.



Tra le tecniche maggiormente usate che affiancano l’analisi degli isotopi stabili c’è l’analisi macro e micro-petrografica e la diffrazione a raggi X. Queste tecniche permettono di stabilire le caratteristiche strutturali del campione, i minerali presenti e la tessitura. Nelle immagini seguenti vengono riportate le caratteristiche petrografiche di alcuni marmi osservati con un microscopio polarizzato.

Un’altra possibile applicazione degli isotopi stabili è quella della determinazione della provenienza di oggetti in metallo. In questo caso si usa l’isotopo del piombo i cui studi risalgono al 1970. La provenienza geologica del metallo può essere determinata grazie al frazionamento a cui viene sottoposto il piombo durante i processi metallurgici, come la fusione, raffinazione, lavorazione, colata o corrosione. Il piombo è presente in piccole quantità nel bronzo, nelle monete e in tanti altri oggetti metallici. In natura esiste come 4 isotopi stabili con numero di massa 204, 206, 207 e 208. Per stabilire la provenienza di un oggetto metallica si usa un diagramma del rapporto 206Pb/204Pb verso il rapporto 207Pb/206Pb o 208Pb/206Pb. In questo modo vengono riportati molti dati di miniere di diverse regioni del mediterraneo. Sul diagramma si formano diversi clusters a seconda delle diverse miniere. Qui di seguito un esempio per oggetti in rame da 170 posti del Sud-Ovest dell’Asia. Ci sono alcune punti che si sovrappongono e rendono difficile individuare in maniera univoca la regione della miniera se non si ricorre ad altri metodi.


Una migliore situazione e’ stata ottenuta per esempio per le miniere e oggetti di rame del periodo carcolitico (2900-1800 aC) della Serbia e della Bulgaria dove la sovrapposizione tra i diversi clusters è molto piccola il che si traduce in una migliore precisione del metodo degli isotopi stabili di Pb.



sabato 4 aprile 2020

Alla scoperta delle leggi di fisica con il machine learning



Nel lontano 1999, pubblicai sul TI Technical Journal (March 1999) un articolo scritto con alcuni miei colleghi dal titolo “Semiconductor Device Modeling using a Neural Network” (Link) in cui si proponeva l’uso di neural network (NN) per cercare di predire alcuni parametri critici dei dispositivi a semiconduttore partendo dai dati prodotti da un simulatore TCAD. Ma perche’ usare una NN per fare quello che puo’ fare un simulatore TCAD molto sofisticato e disponibile in tutti i centri di ricerca di grosse compagnie? Questione di memoria. Il CAD per produrre un risultato ha bisogno di tanta memoria e di tempi lunghi, addirittura di intere giornate. Una rete invece occupa poca memoria e una volta che ha appreso e’ velocissima. L’unica cosa importante e’ che con il simulatore bisogna esplorare uno “spazio delle fasi” del modello quanto piu’ largo e’ possibile,  all’interno del quale poi la rete fara’ le sue predizioni. Dopo 21 anni ritrovo in giro la stessa idea, allargata ovviamente a tutti gli algoritmi di machine learning e non solo le reti NN,  e applicata non solo nel campo dei semiconduttori ma ai concetti di base della fisica stessa. Procediamo con ordine.
Nelle ultime due decadi il Machine Learning e’ diventato uno dei pilastri dell’Intelligenza Artificiale, e parte della nostra vita anche se non tutti se ne accorgono. Il termine e’ stato usato per la prima volta da A. Samuel nel 1959. Ma quale e’ la differenza tra un algoritmo tradizionale e uno di machine elarning? Nel primo caso il programmatore conosce il modello e definisce i parametri e i dati necessari alla risoluzione del problema, mentre nel secondo caso non c’e’ un modello a priori e né una strategia. Si fa in modo che il computer impari da solo eseguendo l’attivita’ e migliorandola iterativamente. In questo caso si parla di apprendimento automatico. Con la continua crescita dei dati a disposizione e’ ragionevole pensare che la Data Analytics diventera’ sempre piu’ pervasiva e un ingrediente necessario del progresso tecnologico. Il Machine Learning puo’ apparire in diverse forme come:

1.     Ranking delle pagine web
2.     Identificazione di spamming
3.     Traduzione automatica
4.     Riconoscimento vocale
5.     Riconoscimento immagini
6.     Identificare gli interessi delle persone (come per es. su amazon, netflix etc)

per fare solo pochi esempi.  Il machine learning comunque e’ una tecnica molto potente non solo nel predire quale film ci piace vedere ma anche nella ricerca scientifica. Non a caso viene utilizzato su diversi siti web per la classificazione della forma delle galassie, l’individuazione delle tracce delle particelle elementari dopo una collisione, il riconoscimento di una forma tumorale dall’analisi di immagini TAC o della presenza di infezione da Covid-19 dall’analisi delle immagini a raggi X. E non solo. E’ notizia dell’ultima ora che il machine learning sta aiutando gli scienziati a predire i risultati di alcuni fenomeni naturali e in alcuni casi a rivelare le leggi alla base del fenomeno. Si avete capito bene, le leggi che sono alla base della natura. Pensate un attimo al racconto della mela di Netwon e supponete di trasferirlo nel mondo di oggi. Il nostro Netwon avra’ un cellulare, con cui potra’ organizzare l’esperimento della mela raccogliendo in una tabella l’altezza da cui cade la mela e il tempo necessario per arrivare a terra.  E supponiamo anche che Netwon abbia tanta pazienza da raccogliere ben 5000 misure. Questa potrebbe essere la tabella compilata. Il sistema di misure utilizzato e’ l’MKS cioe’ metri, Kg e secondi. Con y0 Newton ha indicato la posizione iniziale della mela cioe’ l’altezza da cui essa cade. La prima riga per esempio dice che la mela cadendo da un’altezza di 82 metri impiega 4 sec per arrivare a terra e cosi via.
  

 Ipotizziamo che Newton, uomo dei nostri tempi abbia a disposizione un software molto sofisticato per l’analisi dei dati e machine learning,  JMP Pro sviluppato dalla prestigiosa SAS. Tra i diversi modelli a disposizione e’ stata scelta una semplice rete neurale il cui schema e’ riportato qui di seguito:



Una rete neurale e’ un algoritmo che cerca di simulare il comportamento del nostro cervello. E’ costituita da uno strato di ingresso (il nostro y0) da alcuni strati nascosti (nel nostro caso 2, quelli con cerchi verdi) e uno strato di uscita (il nostro tempo t). Ogni strato e’ costituito da diversi nodi che rappresentano i nostri neuroni e delle linee che entrano ed escono da essi che simulano le nostre sinapsi. Ad ogni connessione e’ associato un peso che va a moltiplicare il valore che viene presentato su quel link. La stessa cosa accade per le altre connessioni. Dopo di che si sommano tutti i valori che insistono sullo stesso nodo. Questo costituisce il valore x di una funzione che viene specificata all’interno di ogni nodo (vedi immagine sopra dei cerchi verdi). Questa funzione viene chiamata funzione di attivazione. Tra quelle piu’ note abbiamo la funzione a gradino o sigmoide, quella lineare e quella Gaussiana. A seconda della funzione quindi avremo in uscita un certo valore dato da f(x)=y. 
  

 La potenza di questi oggetti, come per il nostro cervello sta nel mettere insieme piu’ nodi con connessioni non lineari tra loro. Questo costituisce una rete neurale. Dando alla rete sia i valori in input che in output, essa aggiusta in modo opportuno tutti i pesi delle connessioni per far si che l’errore in uscita sia il piu’ piccolo possibile. In partica si stabilisce una funzione di costo (anche detto Mean squared error - MSE) e si cerca dei metodi matematici per minimizzarla. Con la rete mostrata nell’immagine Diagram e con 5000 esempi raccolti dal nostro buon Netwon questo il risultato ottenuto.



Il campione a disposizione come vedete e’ stato diviso in 3 parti: un campione per il training, uno per la validation e un altro per il test.  Questo e’ vero per qualsiasi modello di machine learning utilizzato e non solo per le reti neurali. Questa divisione viene fatta per far si che la rete non impari a memoria perdendo cosi’ in generalizzazione della predizione. In effetti oltre al campione di training e di validazione dove la rete viene ottimizzata, si usa quello di test dove ci sono esempi che l’algoritmo non ha mai visto. Solo se le prestazioni della rete sono buone per entrambi il training e il test allora i suoi risultati si giudicano soddisfacenti. Nel nostro esempio i risultati sono eccezionali come ci dice il coefficiente R2 che e’ praticamente pari a 1 per tutti e 3 i gruppi. Ricordiamo che il parametro R2 esprime la bonta’ del nostro modello. Piu’ questo valore e’ vicino a 1 e piu’ il modello e’ buono. (per approfondimenti si puo’ consultare su questo blog il post ). Qui di seguito l’andamento previsto dalla rete della y0 e del tempo t. Come e’ possibile vedere la rete ha catturato correttamente l’andamento secondo la radice quadrata tra il tempo e l’altezza y0. 



Addirittura se andiamo a verificare il valore del tempo predetto in funzione dell’altezza y0 vediamo che la curva (in rosso) che meglio approssima i dati (in nero) e’ una radice quadrata e il coefficiente della radice quadrata di y0 e’ pari a 0.447 che e’ il valore previsto dalla teoria per l’accelerazione di gravita’.




Infatti la legge oraria della caduta dei gravi e’ data da:

Y=Yo+Vo*t-1/2*g*t^2

dove Yo e’ l’altezza, Vo la velocita’ iniziale e g l’accelerazione di gravita’, una costante pari a 9.8 m/sec2 in prossimita’ della superficie Terrestre. Per l’esperimento realizzato dal nostro Newton del 2020, Y=0 e Vo=0, cioe’

0=Yo-1/2*(9.8)*t^2

da cui

t=sqrt(Yo/5)=1/sqrt(5)*sqrt(Yo)=0.447*sqrt(Yo)  

avendo eliminata la soluzione negativa (in quanto il tempo e’ una quantita’ sempre positiva).
Adesso complichiamo leggermente le cose utilizzando l’intera legge oraria scritta precedentemente assumendo che la velocita’ iniziale non sia zero e che non ci sia un suolo su cui l’oggetto che cade possa fermarsi. Questo spiega i valori della y negativi nella tabella sottostante. I dati sono stati creati utilizzando la funzione RAND() di excel. Abbiamo introdotto anche una piccola variazione percentuale della costante g per cercare di confondere l’algoritmo. Come potete vedere e’ stata introdotta anche la massa m anche se dalla teoria sappiamo che essa non ha alcun impatto sulla caduta dei gravi come invece pensava Aristotele. Vediamo cosa succede utilizzando i predictive models di JMP Pro. 




Oltre alla rete neurale gia’ introdotta prima, sono stati utilizzati i seguenti modelli di machine learning: Random Forest (RF), Ensemble Boosted e KNN.
Prima di poter parlare di Random Forest dobbiamo introdurre il modello di decison tree da cui l’RF deriva. L’idea alla base dei modelli Tree e’ molto semplice. E’ quella dei mitici romani: dividi e governa. Dato un dataset con almeno una risposta Y (sia numerica che categorica) e tante features X (o anche predictors), l’algoritmo cerca di dividere il dataset in due gruppi per una certa X che massimizza la differenza tra i valori medi della Y nei due gruppi e che all’interno di ognuno di essi minimizza la standard deviation. Fatto cio’ ripete piu’ e piu’ volte questa divisione costruendo un vero e proprio albero con tanti rami e tante foglie (in giallo nell’immagine qui sotto).




 Purtroppo i decision tree sono molto sensibili ai dati in ingresso. Se i dati di training vengono cambiati (per esempio addestrando l’albero su un sottoinsieme dei dati di addestramento) l’albero decisionale risultante puo’ essere abbastanza diverso e con diverse previsioni.  Per questo motivo si e’ arrivati alle Random Forest (RF). Come il nome indica si tratta di un grande numero di alberi decisionali che operano come un tutt’uno nel senso che lavorano in parallelo. Le RF si basano sulla regola che un gruppo di persone e’ piu’ intelligente di un singolo individuo (saggio di J. Surowiecki, La saggezza della folla). Ogni singolo tree fa la sua previsione e la classe che ottiene piu’ voti diventa la previsione globale dell’algoritmo nel caso di classificazione oppure il valore medio nel caso di regressione. Individualmente, le previsioni fatte dai singoli alberi potrebbero anche essere non accurate, ma combinate insieme, si avvicineranno alla risposta corretta.
Passiamo adesso all’algoritmo di boosting. Come per le random forest anche questi modelli fanno parte dei cosiddetti modelli previsionali ottenuti tramite la composizione di vari modelli piu’ semplici (ensemble models). Il tutto nasce nel 1988 quando si capisce che un insieme di modelli di apprendimento deboli se messi insieme possono creare un modello robusto (l’unione fa la forza). Il concetto non e’ molto diverso dall’ottenere una funzione complessa a partire dalla somma di funzioni elementari semplici. Vediamo un esempio. Supponiamo di avere a disposizione dei dati (punti in blu nel grafico qui sotto) e vogliamo trovare la migliore curva che approssima questi punti. Poiche’ si vede un andamento globale a crescere la prima cosa che possiamo fare e’ provare con una funzione radice quadrata (in arancione). Non male. Osservando bene pero’ possiamo notare che ci sono delle oscillazioni intorno al valore della radice quadrata. Viene quindi spontaneo aggiungere alla radice quadrata una funzione seno (linea nera) che migliora l’approssimazione come si puo’ vedere dalla tabella excel riportata.



 


La somma dei residui al quadrato 

SUM(Y_actual-Y_predicted)^2

nel caso della radice quadrata e’ di circa 3800 mentre quella della funzione somma  

f(x)=a*sqrt(x)+b*(sin(c*x))                                                                a,b,c costanti

ha un valore di circa 200 cioe’ un fattore x19 di meno.
 E’ chiaro quindi che questo algoritmo lavora in serie cercando di migliorare un learner debole applicandone uno nuovo e cosi via. Allo stesso modo per un problema di classificazione o di regressione se M(x) e’ il primo modello avremo:

            Y=M(x)+e1                             con un’accuratezza per esempio pari a 84%

dove con e1 abbiamo indicato l’errore. Invece di procedere col costruire nuovi modelli, quello che si puo’ fare e’ cercare di modellizzare l’errore ottenendo qualche cosa come

e1=H(x)+e2

e a sua volta

e2=G(x)+e3

Se ci fermiamo qui otteniamo

Y=M(x)+H(x)+G(x)+e3

Arrivati qui possiamo assegnare degli opportuni pesi ai 3 learners M, G e H cercando di ottenere un’accuratezza migliore del primo learner M, cioe’

Y=aM(x)+bH(x)+cG(x)+e4                                                con Accuratezza>84%

L’ultimo modello utilizzato e’ il KNN cioe’ K-Nearest Neighbors model. Si tratta di uno dei modelli piu’ semplici del machine learning che produce dei buoni risultati in diverse situazioni. Si tratta di un algoritmo che cerca di predire una nuova istanza conoscendo i punti che sono separati in diverse classi. L’idea di base e’ quella di usare i k punti “piu’ vicini” (nearest neighbors) per effettuare la classificazione.



Questi algoritmi per poter lavorare correttamente hanno bisogno di una classe di training e di una metrica per calcolare la distanza (per esempio la distanza Euclidea) tra i vari records e il numero di vicini da utilizzare. Il processo di classificazione prima di tutto calcola la distanza rispetto ai record nel training set, identifica i k records piu’ vicini e utilizza le labels delle classe dei primi vicini per determinare la classe del record sconosciuto (per es. scegliendo quella che compare con maggiore frequenza) nel caso di classificazione oppure prendendo il valore medio dei primi vicini nel caso di regressione.  Qui di seguito un’immagine che mostra la definizione di primi vicini per diversi valori di k (da 1 a 3).


E’ chiaro che utilizzando valori diversi di k si otterranno risultati diversi come nel caso della regressione dove il valore medio dipende da quanti vicini vengono considerati. Bisogna quindi stabilire il valore ottimale di k. In genere l’errore per il gruppo di training e di validation ha un andamento come quello mostrato nell’immagine sottostante. Per piccoli valori di k l’errore sull’insieme di training e’ piccolo mentre quello sul validation set e’ grande. Chiaro segnale di overfitting nel senso che l’algoritmo ha una bassa generalizzazione della predizione avendo imparato a memoria. Al crescere di k vediamo che sia l’errore per il training set che per il validation set aumentano indicando che il modello e' in una condizione di underfitting. Comunque osservando la curva dell’errore del validation set (chiamata elbow a causa della sua forma a gomito) vediamo che essa mostra un minimo intorno a k=9. Questo valore di k e’ il valore ottimale del modello KNN.    



Dopo questo excursus sui modelli utilizzati, possiamo ritornare al nostro esperimento sintetico della caduta dei gravi. Utilizzando i 4 modelli descritti e presenti nel modulo di “Predictive modeling” di JMP Pro abbiamo ottenuto i seguenti risultati per la legge oraria del moto completa dei gravi. Il modello con le migliori prestazioni e’ stata la rete neurale seguita dal Boosted, dal KNN e per utlimo dal Random Forest. 



Osserviamo come la rete neurale (la stessa struttura di quella mostrata all’inizio di questo post) sia stata capace di prevedere l’andamento parabolico della y rispetto al tempo t.




Analizzando la sensibilta’ di ognuno dei parametri in X possiamo vedere come la rete abbia capito che la massa non ha alcun influenza sulla caduta di un corpo dando ragione a Newton e non ad Aristotele.



Interessante anche l’andamento della y con la g, l’accelerazione di gravita’ terrestre. Se questa diminuisce, come per esempio nel caso della luna, allora a parita’ di tempo t la y sara’ maggiore cioe’ il corpo avra’ percorso un tratto minore.   Altra osservazione.   Se guardiamo  al  grafico  che  mette in relazione  la y predetta  con quella  reale  si  puo’ vedere che se anche in termini di R2 i  3 modelli di NN,  Boosting e  KNN  sono molto simili essi  mostrano  una  varianza  intorno  alla  linea  centrale  diversa  con  la KNN essendo quella peggiore. 



Il modello della Random Forest e’ quella che merita un discorso a parte. Infatti e’ quello che mostra le performance peggiori in termini di R2. Proviamo a fare un grafico con i valori della y reali e quelli invece predetti dall’algoritmo. Notiamo subito che la retta e’ spezzata in due parti. Una pendenza piu’ o meno simile per i valori positivi e negativi con diversa intercetta. Anche considerando il valore assoluto della y questa anomalia rimane. Questo e’ un andamento che gli altri modelli non mostrano.




Passiamo adesso ad un altro fenomeno fisico molto noto, quello dell'interazione gravitazionale. Anche qui grazie a Newton sappiamo che la legge con cui due masse M1 e M2 si attirano dipende dalla loro distanza al quadrato secondo la legge:

F=-G (M1M2/r^2)

dove G e’ la costante gravitazione un numero uguale in tutti i punti dell’universo ed r la distanza tra le due masse. Come nei casi precedenti abbiamo generato 5000 records aggiungendo al data set una colonna V (velocita’) per cercare di “ingannare” il modello.
Esattamente come prima il miglior learner risulta essere la rete neurale, seguita dal modello Boosted, dal KNN e in ultimo dalle Random forest con una performance veramente poco convincente. 





Essendo l’R2 veramente molto basso si e’ provato a cambiare i parametri di default di JMP per migliorare l’R2. Uno di essi si e’ dimostrato essere quello giusto. Passando dal valore di default 2 a 5 siamo passati da un 47.7% al 91% come mostrato nella figura sottostante.  Il parametro in discussione e’ il numero di predittori che la foresta utilizza  ad ogni split. Avremo quindi una foresta che produce i suoi alberi considerando un solo predittore ad ogni split. Poi una seconda foresta che ne considerera’ 2 e cosi via fino a 5.  

Ancora una volta i learner sono riusciti a catturare il corretto andamento della legge fisica come possiamo vedere dall’andamento come 1/r^2 della forza F. E di nuovo i modelli hanno interpretato correttamente l’impatto nullo della velocita’ sulla forza F come stabilito dalla legge gravitazionale di Newton.



Questa la rete neurale utilizzata con due strati nascosti e tre diverse funzioni di attivazione.


  
Analizzare e comprendere i risultati dei diversi modelli di machine learning ancora richiede l’intervento umano essendo i risultati non chiari al pari di una legge matematica con tutta la sua eleganza e bellezza. Ma di sicuro essi possono aiutarci a capire come si comportano le leggi del nostro Universo semplicemente guardando nei dati, la ricchezza del futuro, il petrolio che fara’ muovere la tecnologia e la scienza dei prossimi anni. Pensate per esempio alla possibilita’ di applicare il machine learning a problemi complessi come quelli dei 3 corpi o all’ipotesi di Riemann sugli zeri della funzione zeta per citarne solo alcuni. 


  
http://www.wikio.it