Regressione NON lineare

boba74
Salve sono nuovo del forum
Ho un problema di tipo statistico. Ho un modello del tipo:

$Y=A+B*e^(C*X)$

dispongo di n dati (x,y) dove n>3 e dovrei stimare i parametri A, B, C mediante una regressione che in questo caso non è lineare (e non posso passare ai logaritmi, a causa del termine "A").
Esiste una formula ai minimi quadrati per la stima di A,B,Cc oppure si devono usare metodi numerici? In tal caso, esiste un programino semplice che mi permetta di risolvere il sistema?
Inoltre, già che ci siamo, dovrei anche fare i test di significatività sui parametri....

Grazie

Risposte
olaxgabry
Il fatto che tu non possa usare il metodo dei minimi quadrati ordinari, dipende sì dal termine $A$ ma non necessariamente: mi spiego meglio. Considera il modello di regressione

$Y=Be^{CX}*epsilon$

dove $epsilon sim N(0,sigma^{2})$: in maniera pratica $epsilon$ rappresenta l'errore che commetti in fase di stima o previsione. In questo caso tu puoi passare ai logaritmi ed ottenere il modello di regressione lineare

$ln(Y)=ln(B)+CX+ln(epsilon)$

Qui usi il metodo deiminimi quadrati per stimare $ln(B)$ e $C$. Le cose cambiano se il tuo modello è del tipo additivo, ovvero

$Y=Be^{CX}+epsilon$

in quanto l'uso dei logaritmi non ti porta a dei miglioramenti: in questa circostanza siamo nella stessa situazione del modello che proponi tu in quanto non possiamo stimare i parametri con il metodo dei minimi quadrati ordinari.
Cmq, tornando al tuo modello, per stimare i parametri si usano delle iterazioni numeriche come il metodo dei minimi quadrati non lineari: il primo problema è risolto.
Dato che devi fare inferenza sui parametri, sappi che se gli errori $epsilon$ sono distribuiti normalmente, allora gli stimatori dei parametri $,hat{A},hat{B},hat{C}$ sono consistenti, efficienti e hanno una distribuzione limite normale, per cui puoi usare tranquillamente i "classici test": ovvero valori della $t-Student$ o, se $n$ è grande, quelli della Normale.
Per quanto riguarda i software, puoi provare ad usare gretl: è gratuito ed ha un'opzione (sotto modelli) propio sui modelli non lineari; è facile da usare e poi, per ogni cosa tu voglia fare, appare sempre una finestra aiuto che ti dice come procedere.
Il link è
http://gretl.sourceforge.net/gretl_italiano.html

Al limite puoi anche trovarlo qui

www.freestatistics.altavista.org (qui ci sono tutti i software econometrici gratuiti: un altro carino è jmulti. Il più potente è R, ma non è semplicissimo).

Spero di esserti stato d'aiuto. Piccola curiosità: il modello lo devi usare per serie storiche o dati cross section?
Ciao

boba74
Grazie 10000!!!!
Ora provo un po' a smanettare.

Il modello devo usarlo per scopi "botanici", ovvero, per determinare il legame tra esposizione al freddo (in questo caso X) e la successiva esposizione al caldo(in questo caso Y) di una data specie vegetale, necessarie per la rimozione della dormienza invernale delle gemme. Spiegazione botanica:
Le specie cosiddette "dormienti" hanno un periodo di arresto delle gemme che inizia in autunno e termina in primavera, quando la dormienza è completamente rimossa e le gemme si possono aprire. Tuttavia le gemme non percepiscono il passare del tempo, sono solo in grado di "leggere" le temperature. In certe specie (quelle che seguono questo particolare modello), c'è una prima fase in cui è il freddo a determinare la rimozione della dormienza, in partricolare la variabile X (o CD, da Chill days), che convenzionalmente è il numero di giorni con temperatura media < 5°C contati a partire dal 1° novembre fino alla data di germogliamento. Successivamente nella seconda parte dell'inverno il freddo da solo non è più sufficiente a rimuovere la dormienza, la quale viene rimossa dall'esposizione al caldo, in particolare dalla quantità Y (o TT, Thermal Time), che convenzionalmente è la somma delle temperature giornaliere >5°C a partire dal 1° febbraio e fino alla data di germogliamento.
Si è visto che in queste specie il legame è del tipo esponenziale inverso (come da mia equazione), ossia che maggiore è il "freddo" percepito nella prima parte dell'inverno, e minore il fabbisogno residuo di calore necessario a rimuovere la dormienza, e viceversa. (alcuni alberi tipo il Faggio, se non ricevono abbastanza freddo rischiano di svegliarsi in ritardo, alcuni abeti anche mai, perchè la curva diventa talmente ripida che a bassi valori di X si ottengono valori eccessivamente alti di Y, non compatibili con le condizioni presenti in natura.....
Io parto da dati noti di X e Y di una data specie in luoghi e anni diversi, ricavati dalle T giornaliere misurate, e devo quindi determinare la curva caratteristica di quella specie.... L'esperimento è a scopo bonsaistico, per me solo un hobby, perchè io faccio tutt'altro.... 8-)

Ah, aggiungo. Dato il modello, il parametro A rappresenta il minimo fabbisogno di calore a prescindere dall'esposizione al freddo, perciò non può essere nullo. Al contrario, esistono specie che non hanno necessità di esposizioneal freddo, in tal caso la curva si riduce a una retta orizzontale costante pari ad A. La significatività sui parametri mi serve quindi anche per stabilire se il modello è applicabile a una data specie, o se invece si possa ridurre a una relazione lineare tra X e Y, o a una costante X=A.
Poi appena avrò dati disponibili sarò certamente a richiedere chiarimenti, perchè di statistica ne ho vista pochissima, essendo ingegnere. 8-)

olaxgabry
Interessante. A proposito: dato che hai detto che di statistica ne hai vista poca, per caso sai come effettuare i test di significatività sui paremetri?
Cmq, se prendi come livello di significativa $alpha=0,05$, affinché i tuoi parametri siano significativamente diversi da zero ti basta che il modulo del rapporto tra la stima del parametro e la sua deviazione standard sia maggiore di 2. In gretl comprare sotto il nome rapporto t o statistica t.
Se hai problemi posta e vediamo se ti posso aiutare.
Ciao e buona regressione!

boba74
sì, ho fatto alcune regressioni lineari su dati pluviometrici durante la mia tesi, e ho applicato la procedura della t di Student per la significatività, anche se devo dire che in alcuni punti la mia memoria è un po' arrugginita...
A questo proposito avrei una domanda: ho un articolo in cui si parla di un esperimento fatto proprio con questo modello su una serie di dati, e alla fine c'è una tabella dove si vedono le varianze e i test di significatività su questi parametri. Bene, in pratica risulta che solo il parametro B è significativo con p<0,05, mentre A e C non sono significativi (cioè hanno un intervallo di confidenza molto grande che comprende anche lo zero....), la varianza spiegata è circa 89%.
Non riesco a capire cosa implica questo risultato: vuol dire che posso scartare i parametri A e C, oppure il fatto che non siano significativi è solo dovuto ai pochi dati usati per stimare i parametri? (n=10) Perchè così "a occhio" mi sembra che i dati si adattino bene al modello esponenziale (R^2=0,7), mentre se utilizzassi un modello lineare avrei una dispersione molto maggiore...

olaxgabry
Non avendo i dati e la loro visualizzazione sul piano, ti parlo per esperienza. Da quello che mi hai detto, i parametri $A$ e $C$ del modello

$Y=A+BE^{CX}$

sono risultati non significativi (quindi con un p-value diciamo sopra a 0,1). Hai due parametri non significativi: quello che si fa solitamente è prendere quello meno significativo e eliminarlo dal modello, ristimando nuovamente i parametri del "nuovo" modello (in questo caso, però, data la tua funzione ha senso fare ciò se A fosse il meno significativo: magari mandami i valori per questi due parametri.). Poi un consiglio: attento a non fare troppo affidamento all' $R^{2}$; infatti, quando elimini la costante, ad esempio dal modello lineare classico, l'$R^{2}$ può portarti a considerazioni errate.
Poi può darsi che effettivamente l'esponenziale non sia la scelta migliore: solitamente si testano anche altre funzioni, tipo la logistica o l'esponenziale modificata.
Ciao

boba74
Ecco il grafico. Il modello è stato applicato su 3 specie, ottenendo 3 stime di parametri (sul grafico sono riportati tutti, ma usando simboli diversi, anche se alla fine non è molto leggibile.....):



E questa è la tabella riassuntiva:



Come si vede, per tutti e 3 le specie, i parametri A e C (in questo caso chiamati "a" ed "r") risultano non significativi, anche se trattandosi di unità di misura diverse non so dire qual'è il "meno signficativo" dei due....
Comunque come dicevi, l'unico che posso togliere dal modello è "a", perchè altrimenti la formula si ridurrebbe a una costante....
Secondo te potrebbero esserci altri modelli possibili?

olaxgabry
Se prendiamo come esempio U.glabra, il parametro meno significativo è $a$ per un livello di significatività $alpha=0,05$ (basta che prendi il modulo del rapporto tra il coefficiente e la sua deviazione standard e poi prendi quello minore). Poi bisogna vedere se il parametro $a$ sia effettivamente necessario: però, se così fosse, non avrebbe senso fare il test di significatività dato che deve essere presente nel modello.
Io proverei con altri tipi di regressione: prova con quelle che ti ho proposto.
Al limite, anche se i tuoi paremetri non dovessero essere significativi (bada che puoi sempre commettere un'errore di prima specie) ma il tuo grafico ti sembra buono, prova a vedere come si comporta in fase di stima: al limite ti consiglio alcune misure per capire la bontà del tuo modello (il più usato è il MAPE che sarebbe la media aritmetica, in percentuale, degli errori relativi in valore assoluto).
Cmq, se hai voglia, postami il grafico dei tuoi dati.
Ciao

boba74
purtroppo io avrò dati disponibili solo a 2009 inoltrato... forse anche 2010, perchè le misure di temperatura fanno riferimento a quest'inverno e arrivano fino a primavera prossima, tuttavia già da maggio-giugno potrei fare delle "prove" su alcuni dati che già possiedo, anche se inizialmente dovrò accontentarmi di un n molto basso....
Comunque ti ringrazio della tua disponibilità, per ora mi basta informarmi sulla teoria, in modo da sapere poi come procedere una volta che i dati saranno disponibili.
Comunque, riflettendo sui dati, si vede che in effetti il parametro a ha un significato fisico di "minimo valore di Thermal Time" necessario al germogliamento, anche se per alti valori di Chill days i dati sembrano non avere un limite inferiore (vedi punti più a destra nel grafico), mentre il modello tende a sovrastimare questo fabbisogno. Può anche darsi che esposizioni al freddo molto elevate (fuori dai valori registrati dal campione), possano darmi valori ancora più piccoli di TT, al punto che potrebbe anche essere a=0. In tal caso diventerebbe plausibile eliminarlo dal modello....
Tuttavia io cercherò anche altri modelli, perchè secondo me per alcune specie potrebbe anche esserci un valore minimo di Chill days (a sinistra), per il quale la curva in oggetto tendere asintoticamente a TT=infinito (ossia un minimo fabbisogno di freddo sotto il quale non è possibile andare) allora in tal caso serve un'altro modello, chessò, un iperbole....

orbit1
ragazzi un aiuto,è la prima volta che ho tra le mani statistica quindi scusate l'ignoranza
cmq nell'adattamento dei modelli non lineari,dopo che riscrivo l'equazione con i logaritmi come calcolo i nuovi a e b?
esempio. y=(a)moltiplicato per (e^bx) -> (lny)=(lna)+bx -> w=A+bx ora quando devo calcolare A e b le calcolo cosi
A=M(W)-bM(X) b=cov(X,W)/varianzaX
ma dove vado a prendere i valori di W? sono quelli di y?

olaxgabry
"orbit":
ragazzi un aiuto,è la prima volta che ho tra le mani statistica quindi scusate l'ignoranza
cmq nell'adattamento dei modelli non lineari,dopo che riscrivo l'equazione con i logaritmi come calcolo i nuovi a e b?
esempio. y=(a)moltiplicato per (e^bx) -> (lny)=(lna)+bx -> w=A+bx ora quando devo calcolare A e b le calcolo cosi
A=M(W)-bM(X) b=cov(X,W)/varianzaX
ma dove vado a prendere i valori di W? sono quelli di y?

Allora hai il modello moltiplicativo

$Y_i=a*e^(b*X_{i})*\epsilon_{i}$ $i=1,...,n$

dove $a,b$ sono costanti ed $\epsilon$ è un white noise. Usando la trasformata logaritmica hai

$\ln(Y_{i})=\ln(a) + b*X_{i} + \ln(\epsilon_{i})$

Poni $W_{i}=\ln(Y_{i})$ e $A=\ln(a)$, $u_{i}=\ln(\epsilon_{i})$ per cui il modello diventa

$W_{i}=A + b*X_{i}+u_{i}$

Su questo modello puoi utilizzare gli OLS. Ora il valore di $W_{i}$ è semplicemente il logaritmo in base $e$ di $Y_{i}$.
Chiaro?

orbit1
si grazie 1000!me ne sono reso conto subito dopo averlo scritto XD avevo pure annotato lny=w mi sara sfuggito x distrazione :oops:

olaxgabry
Figurati, per così poco.
Ciao e buona regressione.

Rispondi
Per rispondere a questa discussione devi prima effettuare il login.