Salta al contenuto

Ricostruiamo l’Hockey Stick con R

Ultimamente ci siamo occupati di analizzare le temperature e la CO2 tramite il programma open source R. Presto pubblicheremo altre analisi e metodi statistici più raffinati. In attesa di questi nuovi tutorial, diversamente dal solito, oggi analizzeremo il programma scritto da Jeff Id che cerca di replicare il sistema utilizzato per creare l’hockey stick.

Prima di procedere, ricordo che R è disponibile sia per Microsoft Windows che per Linux.

La domanda che molti si staranno ponendo è: perchè tornare di nuovo sull’Hockey stick? La risposta non è univoca e possiamo dire che principalmente si torna sull’argomento in quanto l’ultima versione dell’HS è stata edita nel 2008 e poi perchè è sempre una notevole fonte di riflessione sulla ricerca scientifica. In questo articolo, però, affronteremo il “semplice” aspetto scientifico.

Secondo l’autore del programma1 (scritto in R) che andremo ad analizzare, quello che è successo con l’HS è molto semplice: dal set iniziale di proxy sono stati eliminati quelli che mostravano una correlazione debole con le temperature. La correlazione viene calcolata tra gli anni più recenti dei proxy e quelli più recenti delle temperature. Questa operazione preliminare ha fatto sì che Mann scartasse, dal set iniziale di 1209, ben 725 proxy che non rispondevano al criterio imposto di correlazione. Su questo fatto torneremo più avanti, perchè abbiamo visto di recente come i proxy siano spesso impossibili da legare in alcun modo alle temperature.

La fase successiva prevede di calibrare la deviazione standard dei proxy selezionati al punto precedente, in modo tale che questa coincida con la deviazione standard del segnale di controllo. In questo modo i proxy selezionati avranno delle oscillazioni che coincideranno molto bene con le oscillazioni delle temperature.

Infine calcoliamo la media del set di proxy selezionati e fin qui elaborati.

Stiamo per verificare, tramite questo programma in R, come il succitato metodo consenta di calibrare i set di proxy in modo tale che il segnale di controllo emerga sempre e comunque. Qualsiasi segnale di controllo si stia cercando. Questo è il punto focale di tutta la ricerca, ovvero il metodo CPS (Composite Plus Scale) che è stato inventato ad hoc per gestire un elevato numero di dati proxy e per ricreare un qualsiasi segnale di controllo. L’autore del programma intravvede malafede in questo metodo, noi semplicemente sottolineiamo come spesso in climatologia, non paghi dei metodi statistici correnti, ci si avventuri nell’invenzione di qualche nuovo filtro o metodo di regressione.

Un avvertimento per i gentili lettori: stiamo per iniziare una trattazione piuttosto lunga e corposa, il consiglio è di scaricare lo script e di tenere la vostra versione di R aperta e a portata di mouse.

Introduciamo l’argomento con un semplice esempio, anche per chi non mastica troppo di statistica. Per fare ciò utilizzeremo uno strumento che oggi sempre più persone hanno sul proprio computer di casa: un elaboratore di immagini. La tecnologia oggi è molto raffinata e ci consente di mettere sempre più spesso una “toppa” a fotografie scattate male o distrattamente. Ipotizziamo di fotografare una nuvola e di sovraesporla (un classico): vogliamo il cielo blu intenso e purtroppo la nuvola si brucia completamente. Ebbene, gli editor digitali di fotografie ci consentono di rimediare al nostro errore, andando a correggere il contrasto complessivo della foto. Se siamo fortunati, la sovraesposizione ha lasciato qualche dettaglio (in termini digitali, abbiamo ancora delle informazioni da poter elaborare). Niente di meglio per il programma che va a caccia di quelle poche informazioni presenti e riesce a restituirci i dettagli della nuvola. Se siamo stati sfortunati (o di mano troppo pesante), all’interno della nuvola avremo pochissime o nessuna informazione e, a quel punto, nessuno programma potrà mettere in luce alcun dettaglio.

Tutto questo preambolo per dire cosa? Uno dei pilastri fondanti dell’indagine statistica è la presenza di dati. Se abbiamo un campione numeroso, ci possiamo lavorare tranquillamente. Se, diversamente, non abbiamo un campione numeroso, anzi addirittura abbiamo scarsità di dati, non c’è algoritmo che tenga: non possiamo ricostruire qualcosa da una serie rarefatta. Mutatis mutandis, da una serie rarefatta possiamo tirare fuori qualsiasi cosa.

Ecco cosa succede con il metodo statistico sopra descritto: dai proxy così selezionati e scalati, in realtà possiamo ottenere qualsiasi tipo di segnale, semplicemente andando a cambiare il segnale di controllo. Volete temperature in aumento a fine millennio? Nessun problema. Volete temperature (sempre con gli stessi proxy) in diminuzione a fine millennio? Neanche questo è un problema. Infine, volete temperature che oscillano a fine millennio? Nemmeno a dirlo, si può fare anche questo.

Il listato originale lo potete trovare a questo link, per una migliore comprensione, tuttavia, lo analizzeremo nei suoi punti più interessanti.

Inizializziamo lo script e carichiamo i dati. Se caricate lo script con una connessione internet attiva, il programma scaricherà automaticamente tutto quanto necessario per eseguire i calcoli (diversamente dovrete scaricarvi i dataset manualmente e andare a modificare il percorso sul vostro disco rigido). Ricordo che lo script è stato elaborato da Jeff Id.

###########################################################################
##################### TURNKEY DATA DOWNLOADER #############################
###########################################################################
method=”online”
if (method==”offline”) {load(”d:/climate/data/mann08/mann.tab”)
load(”d:/climate/data/mann08/mannorig.tab”)
load(”d:/climate/data/mann08/idorig.tab”)
names(mannorig)=idorig$idorig
load(”d:/climate/data/mann08/details.tab”)
source(”d:/climate/scripts/mann.2008/instrumental.txt”)} else {
url=”http://www.climateaudit.org/data/mann.2008″
download.file(file.path(url,”mann.tab”),”temp.dat”,mode=”wb”);load(”temp.dat”)
download.file(file.path(url,”mannorig.tab”),”temp.dat”,mode=”wb”);load(”temp.dat”)
download.file(file.path(url,”details.tab”),”temp.dat”,mode=”wb”);load(”temp.dat”)
download.file(file.path(url,”idorig.tab”),”temp.dat”,mode=”wb”);load(”temp.dat”)
names(mannorig)=idorig$idorig
source(”http://www.climateaudit.org/scripts/mann.2008/instrumental.txt”)}

Nulla da dire fin qui, carichiamo una serie di dataset ospitati sul sito Climate Audit, in buona sostanza è l’insieme dei proxy necessari per l’analisi.

############################################################################
################ FUNCTIONS ################################################
############################################################################

######## CREATE TIME SERIES ##########
gg = function(X)
{
ts(X[,2],start=X[1,1])
}

######## GAUSSIAN FILTER FUNCTION TESTED AT CA #########
######## set to 11 years ###############################
source(”http://www.climateaudit.org/scripts/utilities.txt”)

ff=function(x)
{
filter.combine.pad(x,truncated.gauss.weights(11) )[,2]
}
##########################
##### END FUNCTIONS ######
##########################

Con questo secondo blocco creiamo una serie di funzioni che ci torneranno utili nel corso dell’elaborazione. La prima funzione “gg” imposta lo schema base per creare le serie storiche (il comando ts è già stato analizzato a fondo su Climate Monitor, vi invito a dare un’occhiata ai tutorial pubblicati). La funzione “ff”, invece, si occupa di filtrare le serie storiche. Il filtro applicato è un gaussiano centrato su periodo di tempo di 11 anni2 .

mano=mannorig #copy non-infilled PROXY series to different memory
#mano=mann #copy REGEM infilled PROXY series to different memory

In questo blocco abbiamo la prima scelta possibile. La variabile “mano” è di default settata per caricare i proxy “mannorig” (serie non completata tramite il metodo RegEM), tuttavia è possibile scegliere la serie “mann” che al contrario è completata con il sistema RegEM. Per fare questo basta deselezionare la riga “mannorig” con il simbolo # ed eliminandolo dalla riga “mann”.

for(i in 1:lent)
{
mano[[i]][,2]=ff(mano[[i]][,2])
}

Il blocco precedente è molto interessante, in quanto ci consente di filtrare ricorsivamente tutti i proxy con la funzione “ff” che abbiamo incontrato precedentemente. E’ un normale ciclo “for” utilizzato in modo sapiente all’interno di R.

#m=c(1:leng)/-leng ###### 0 to -1 NEGATIVE linear slope
#m=c(1:leng)/leng ###### 0 to 1 POSITIVE linear slope
#m=sin( (1:leng)/(leng/cycles)*pi) #sin wave with # cycles
m=-sin( (1:leng)/(leng/cycles)*pi) #-sin wave with # cycles

Questa sezione è molto importante perchè troviamo le funzioni che riproducono il segnale di controllo. Ne abbiamo di diverso tipo: temperature crescenti, decrescenti e oscillanti. Anche in questa sezione possiamo operare una scelta (l’invito è provarle tutte per vedere come le serie proxy si adattano perfettamente a ogni tipo di segnale).

Jeff Id opera una serie di calcoli preparatori che qui tralasceremo, per passare direttamente al cuore di tutto il programma: il Composite Plus Scale (CPS).

for(j in 1:lent)
{
if(pass[j]==TRUE)
{
y=(gg(mano[[j]]))
val = window(y, start(m)[1],end=end(m)[1])
m0=mean(val,na.rm=T)
sd0=sd(val,na.rm=T)

################# CPS METHOD IS APPLIED HERE #######################

yscale=y-m0
yscale=yscale/sd0
yscale=yscale*ssig
yscale=yscale+msig

x=ts.union(x,ts (yscale,time(y)[1]))
}
}

In questo ciclo “for” molto importante, innanzitutto scremiamo i proxy, scartando quelli che non rientrano nella correlazione desiderata con le temperature, quindi le serie vengono troncate negli anni in cui esistono i dati del segnale di controllo. L’abbiamo detto sopra, vengono utilizzati gli ultimi anni disponibili nei proxy e gli ultimi anni disponibili nel segnale di controllo delle temperature. Le operazioni di centratura e scalatura vengono effettuate solo sulla base di questo periodo.

Calcoliamo quindi la media m0 e la deviazione standard sd0. Adesso siamo pronti per il CPS: sommiamo tutti i proxy, non prima di averli centrati su m0 e scalati con sd0. In particolare per la deviazione standard, essendo a denominatore:

yscale=yscale/sd0

al crescere della deviazione standard (quindi dell’oscillazione complessiva del proxy), cresce la scalatura del proxy derivato. Questo è importante: abbiamo eliminato i proxy che non coincidono con la correlazione da noi cercata, e ora scaliamo in modo più deciso quei proxy che oscillano di più (al di fuori del periodo di calibrazione, non dimentichiamolo!).

Infine i proxy vengono sommati ricorsivamente e assemblati per avere un unico proxy del quale viene poi calcolata la media. Questo è il valore che viene rappresentato nei grafici.

Segnale di controllo crescente

Ricostruzione delle temperature tramite segnale di controllo crescente.

Segnale di controllo decrescente

Ricostruzione delle temperature tramite segnale di controllo decrescente.

Ricostruzione tramite segnale di controllo oscillante

Ricostruzione delle temperature tramite segnale di controllo oscillante.

Reblog this post [with Zemanta]
Related Posts Plugin for WordPress, Blogger...Facebooktwitterlinkedinmail
  1. http://noconsensus.wordpress.com/2009/06/20/hockey-stick-cps-revisited-part-1/ []
  2. Un approfondimento molto tecnico http://www.ltrr.arizona.edu/~dmeko/notes_8.pdf []
Published inAttualitàClimatologiaR

12 Comments

  1. Carlo Manzo

    Secondo me, dovrebbe prevalere il principio di cautela, visto che i due opposti scenari (quello della prosecuzione della politica energetica attuale, e quello che vede invece il contenimento delle emissioni) sono associati a un rischio nettamente diverso. Nel primo caso si prospettano situazioni disastrose per non dire raccapriccianti qualora fosse più o meno corretta la teoria convenzionale sul riscalamento climatico. Cosa succederebbe se, giustamente, anche i paesi emergenti volessero intraprendere questo modello di sviluppo? Oggi il 20% della popolazione consuma l’80% dell’energia (in gran parte derivante dal petrolio). Senza spinte politiche su differenziazione delle fonti energetiche, uso di tecnologie a zero impatto, efficienza energetica di aziende, abitazioni e sistemi di trasporto, mi pare evidente che il modello di sviluppo attuale (senza limiti) non possa essere esteso a tutto il pianeta. E questo porterà inevitabilmente a situzazioni conflittuali ben più ampie di quelle che già stiamo vivendo attualmente.

  2. L’incertezza delle ricostruzioni vale in entrambi i sensi e deriva in larga misura dalla mancanza di completezza delle informazioni, nonchè dall’oggettiva difficoltà di raccordarle alle osservazioni moderne, le quali a loro volta sono incomplete e disomogenee, pur essendo molto più precise. Ad avvalorare la tesi dell’esistenza di periodi con importanti oscillazioni climatiche anche superiori a quella attuale, arrivano però anche le evidenze storiche ed i contributi interdisciplinari. La discussione è però parecchio più complessa, sia della diatriba sulla mazza da hockey di Mann che delle rasoiate di Occam e riguarda soprattutto il fatto che senza contributo antropico, non esisterebbe ragione di preoccuparsi di mitigare il clima, quanto piuttosto di adattarsi ad esso, qualunque direzione dovesse prendere. Assegnare un peso tanto rilevante alla forzante antropica, pur in presenza di enormi incertezze al riguardo, getta le basi per proiezioni la cui validità è tutta da verificare, ma costituisce già oggi strumento di giustificazione di decisioni politiche ed economiche. Queste riguardano tutti indistintamente, specie se prese a scala globale e, per dirla con un’autorevole voce di questo blog, Luigi Mariani, sbagliare la diagnosi mette a serio rischio di sbagliare la cura. Ciò detto, piuttosto che avallare faraonici piani di mitigazione (che hanno già dimostrato di essere scientificamente inefficaci e politicamente inattuabili), sarebbe forse meglio preoccuparsi di ciò che già è, piuttosto che di quello che forse sarà. Ambiente, energia, risorse alimentari ed idriche tanto per cominciare. Lavorando seriamente su queste, scenderanno anche le emissioni e saremo tutti più contenti. Forse lo saranno meno nelle varie cittadelle finanziarie dove si pasteggia a base di emission trading, ma, tant’è tutto non si può avere.
    gg

  3. Carlo Manzo

    Salve,
    premetto che da poco che frequento questo sito e sono più “appassionato” che esperto di clima. Da quanto ho capito, anche leggendo il blog, le tesi negazioniste si fondano soprattutto sulla constatazione dell'”optimum medioevale” secondo il ragionamento per cui, se vi sono state oscillazioni positive di lunga durata nel passato, in assenza di emissioni antropogeniche di gas serra, allora anche il periodo caldo attuale può essere spiegato da fattori naturali (es. variazioni di irradianza solare). Ammesso che siano attendibili le ricostruzioni passate tramite proxy, cosa discutibile come giustamente mostrato in questo post, a me questa pare un’argomentazione debole. Si critica molto l’uso della mazza da hokey di Mann, ma stando così le cose sembra che gli scettici facciano un uso improrpio del “Rasoio di Occam”

  4. @ Pellegrini

    I punti da lei messi in luce sono assolutamente interessanti e costituirebbero una analisi probante (o meno) del lavoro di Jeff Id. R è uno degli strumenti statistici più utilizzati, tuttavia questo tipo di analisi può essere condotta anche con (per esempio) Matlab, o magari scrivendo ex-novo una routine in C++.

    Come per lei, anche per noi il tempo è tiranno e purtroppo (o per fortuna) ogni settimana ricevo decine di mail con richieste di analisi di dati o elaborazioni statistiche. La lista di cose da fare è lunghissima, ahimè. Vedrò cosa si riuscirà a fare.

    CG

    • Luca Galati

      Mi scusi:
      lei solleva il problema dell’eliminazione eccessiva dei dati proxy da parte di Mann (40%) e poi fa statistiche con non più del 20% dei dati, ovvero peggio di quello che avrebbe fatto Mann?
      Non ho capito: dov’è il mio errore?

      Un’altro cosa non ho capito: i proxy di per se ricostruiscono la temperatura passata quando non si è in presenza di dati diretti sulla temperatura stessa: perchè c’è bisogno di confrontarli con la temperatura misurata nell’ultimo secolo?

  5. Giovanni Pellegrini

    Volevo ravvivare un attimo questo post. Penso che sarebbe interessante fare due ulteriori verifiche, che nemmeno jeff id, su the air vent, fa:

    1) Applicare il metodo CPS, ma introducendo poi un metodo a posteriori per controllare i trend su diverse finestre temporali, e tenere soltanto i dati che, calibrati su una finestra, sono coerenti con i dati dell’altra. In questa maniera mi aspetterei che un segnale sinusoidale come quello dell’ultimo esempio venga automaticamente scartato.

    2) Sarebbe inoltre interessante produrre uno script con il metodo EIV e il controllo al punto 1.

    In questa maniera si potrebbe anche stabilire se effettivamente i metodi CPS ed EIV, applicati in tutto e per tutte come li applicati Mann (quindi con il controllo a posteriori su finestre temporali diverse), sono in grado di ricostruire dei segnali di temperatura, al contrario di quello che mostra Jeff Id sul suo blog, in assenza di controllo a posteriori.

    Farei questi conti io, ma proprio non ho tempo e, lo ammetto, R non riscontra i miei favori.

    Cordiali Saluti

    Giovanni Pellegrini

  6. @ Pellegrini

    Teoricamente sarei in ferie, in ogni caso ho letto il suo dubbio.

    Chiaramente il numero di proxy ha un valore statistico, più che fisico.

    A lei viene un dubbio, giustamente, perchè vede tre grafici con una percentuale utilizzata di proxy decisamente più bassa rispetto al 40% di Mann. In realtà non ho volutamente inserito il grafico col 40% perchè, implicitamente e me ne scuso, l’avevo lasciato come “compito” a casa per chi avesse voluto cimentarsi col programma.

    Non ho sottomano il grafico, appena possibile ne metterò una copia, ma le posso anticipare che il comportamento è esattamente lo stesso.

    C.G.

  7. Giovanni Pellegrini

    @Claudio Gravina e Tutti

    Bump!!! Mi scuso di bussare e della mia impazienza! C’e’ qualcuno che sa rispondere alla mia domanda poco sopra?

    Cordiali Saluti

    Giovanni Pellegrini

  8. Giovanni Pellegrini

    @Claudio Gravina

    Ho una domanda riguardo questa ricostruzione.

    Mann nella sua calibrazione riesce a tenere circa il 40% di proxy. Nei tre casi qui esposti vengono tenuti rispettivamente il 20%, il 5% ed infine il 9% dei proxy. Questa differenza di percentuali potrebbe avere un significato fisico?
    Intendo dire che, avendo un numero altissimo di proxy (che ne so, 10000 time series), mi aspetto, prendendo un subset arbitrariamente piccolo, di essere in grado di riprodurre qualsiasi segnale prendendo solo quei pochi proxy che si correlano bene. Ma e’ possibile correlare bene un’alta percentuale di proxy su qualsiasi segnale?

    Cordiali Saluti

    Giovanni Pellegrini

Lascia un commento

Il tuo indirizzo email non sarà pubblicato.

Questo sito usa Akismet per ridurre lo spam. Scopri come i tuoi dati vengono elaborati.

Categorie

Termini di utilizzo

Licenza Creative Commons
Climatemonitor di Guido Guidi è distribuito con Licenza Creative Commons Attribuzione - Non commerciale 4.0 Internazionale.
Permessi ulteriori rispetto alle finalità della presente licenza possono essere disponibili presso info@climatemonitor.it.
scrivi a info@climatemonitor.it
Translate »