Home Forum Statistica con R Come funziona optim()?

This topic contains 6 replies and has 2 voices.

Viewing 7 posts - 1 through 7 (of 7 total)
  • Author
    Posts
  • #4839

    blackema
    Participant

    Salve ho da risolvere un problema di ottimizzazione di parametri teorici di una funzione per adattarla all’andamento dei dati sperimentali.
    Girando in rete ho trovato un paio di suggerimenti che indicano di usare una funzione del pacchetto di base stats che si chiama “optim()”.
    Per caso qualcuno ha esperienza con questa funzione?
    Potreste darmi qualche suggerimento sul suo funzionamento?
    Ringrazio anticipatamente chiunque possa darmi un suggerimento in merito a quanto esposto.

    #4840
    Davide Massidda
    Davide Massidda
    Moderator

    Io la conosco abbastanza. La cosa, per la verità, non è proprio semplicissima. Provo a farti un esempio.

    Immagina di voler individuare i parametri b0 e b1 di un modello lineare:

    y = b0 + b1 * x

    Hai prima di tutto bisogno di una funzione da ottimizzare. Diciamo che i valori migliori per b0 e b1 sono quelli che minimizzano la sommatoria dei quadrati residua, ovvero:

    RSS = sum( (y_estimated – y_observed)^2 )

    Dati due vettori di valori osservati x e y, possiamo utilizzare optim() per individuare i valori di b0 e b1 in grado di minimizzare il valore RSS.

    Dunque, dobbiamo convertire in linguaggio R le formule precedenti:

    RSS <- function(par, predictor, y_observed)
    {
        b0 <- par[1]
        b1 <- par[2]
        y_estimated <- b0 + b1 * predictor
        residuals <- sum( (y_estimated - y_observed)^2 )
        return(residuals)
    }

    L’argomento par rappresenta i parametri da individuare. La funzione optim richiede dei parametri di partenza per b0 e b1, che dobbiamo inventarci noi. Valori di partenza ragionevoli per un modello di regressione sono b0=0 e b1=1:
    start <- c(0,1)

    A questo punto, possiamo lanciare optim:
    optim(par=start, predictor=x, y_observed=y, fn=RSS, method="Nelder-Mead")
    dove x e y sono degli oggetti R che contengono i dati osservati, rispettivamente per il predittore e per la variabile risposta.

    Nelder-Mead è una delle possibili routine di calcolo (è un’implementazione del metodo del simplesso).

    #4842

    blackema
    Participant

    Innanzitutto grazie per la risposta, visto che ad oggi ho trovato abbastanza difficoltà nel capire optim() con degli esempi semplici.

    Le spiego nel dattaglio la mia problematica, praticamente sto scrivendo un codice, per la tesi in campo ambientale, che mi permetta di fare delle analisi statistiche sul campione di dati che ho a disposizione e in più restando sempre in ambiente R volevo stimare l’andamento dei parametri sperimentali (5 per la precisione) con il modello teorico della procedura suggerita dalla metodologia che sto applicando.
    La formulazione teorica citata con la quale devo ottimizzare i dati sperimentali è praticamente la seguente:

    la funzione teorica quindi è data dalla somma di 2 funzioni di densità di probabilità (normale + log-normale), tutto ciò per risalire ai 5 parametri(1.media normale,2.deviazione standard normale,3.media log-normale, 4.deviazione standard log-normale,5. A il parametro di peso) per ottenere alla fine:

    Seguendo un esempio simile che ho trovato in rete avevo provato a impostare il processo di ottimizzazione dei parametri in questo modo:

    fitnormlnorm<-function(par, val) {
       p <- 1/(1+exp(-par[5]))
       return(-sum(log(k*p*dlnorm(val, par[1], abs(par[2]), log = FALSE)+
                         k*tn*(1-p)*dnorm(val, par[3], abs(par[4]), log = FALSE))))
    }
    #valori di partenza
    # Mean 1
    m1=1.8; s1=0.9 
    # Mean 2
    m2=3; s2=0.5
    p=0
    
    par<-c(m1, s1, m2, s2,p)
    
    result2<-optim(par, fitnormlnorm, method="BFGS", val=X,
                    hessian=FALSE, control=list(trace=1))
    
    result2  

    dove par è il vettore dei parametri da ottimizzare, val rappresenta il mio campione in ingresso(X), il problema è che andando ad effettuare una procedura di validazione su un campione(X) di cui conosco i 5 parametri non riesco a ottenere il fitting.
    Il problema per come l’ho impostato potrebbe andare bene oppure sono fuori strada per come è stata impostata la function?
    Se riuscisse a darmi un ulterirore suggerimento mi sarebbe di grande aiuto.Grazie

    • This reply was modified 3 months, 4 weeks ago by  blackema.
    #4844
    Davide Massidda
    Davide Massidda
    Moderator

    Secondo me sei sulla giusta strada e stai lavorando molto bene.

    Questo progetto mi ricorda molto da vicino quello che diverso tempo fa avevo fatto per il pacchetto retimes (qui), la cui principale funzione (timefit) fa una cosa molto semplice: dato un vettore di dati, si assume che questi siano distribuiti secondo una ex-gaussiana, quindi ne si stimano i patametri tramite massima verosimiglianza utilizzando optim.

    Che è esattamente quello che vuoi fare tu, solo che la tua distribuzione di riferimento è una log-normale.

    I passaggi quindi sono questi:
    1) Crei una funzione che, preso un vettore di parametri per la log-normale e un vettore di dati osservati, restituisce il valore della funzione di verosimiglianza.
    2) Passi a optim dei parametri di start, un vettore di dati osservati, e la funzione creata al punto (1); optim si occuperà di individuare i parametri che minimizzano la funzione di verosimiglianza.

    Ora, se il tuo codice è corretto, perché optim non individua dei buoni parametri? L’unica ipotesi che mi viene in mente è che tu stia utilizzando una routine di minimizzazione non adeguata.

    Per esempio, mi ero reso conto che per la distribuzione ex-gaussiana non tutti i metodi di stima erano buoni: in condizioni un po’ critiche ce n’erano alcuni che prendevano degli svarioni assurdi.

    In ogni caso, vorrei segnalarti il pacchetto gamlss, in cui la distribuzione log-normale è già implementata: se non hai l’esigenza di costruire una funzione ex novo, potresti sfruttare quello (posso spiegarti come).

    #4847

    blackema
    Participant

    Inizialmente nel mio lavoro stavo cercando qualcosa nei package messi a disposizione dagli altri utenti ma purtroppo la funzione a cui devo fare riferimento per la procedura di ottimizzazione è la somma di una distr. normale e di una log-norm per cui data la particolarità del caso mi ha costretto ad optare(dopo tantissimi giorni di ricerca) per una procedura come quella esposta.
    Il problema di non essermi potuto confrontare con qualcuno che conoscesse bene lo strumento optim() mi ha un pò demoralizzato visto che quando ho fatto la procedura di validazione di un campione del quale sapevo gia i parametri(i 5 della mia procedura) non riuscivo, anche dando dei valori di start molto prossimi a quelli esatti, una soluzione prossima ai reali valori dei parametri.
    Facendo un ulterirore approfondimento mi sono reso conto che la vicinanza dei parametri alla soluzione reale una volta effettuata l’ottimizzazione veniva influenzata da:
    a)numero di campioni a disposizione(purtroppo nel mio caso sono < 10 dati )
    b)set di dati con valori abbastanza distanziati tra loro

    Queste considerazioni potrebbero essere la fonte della non perfetta ottimizzazione?

    #4848
    Davide Massidda
    Davide Massidda
    Moderator

    Ecco, quando parlavo di condizioni critiche, intendevo proprio questo. Dieci dati sono troppo pochi per individuare i parametri di una distribuzione così complessa… non ce la si può fare.

    Se vuoi capire se il tuo metodo funziona bene, simula campioni di 1000 osservazioni e vedi che cosa restituisce optim. Se le stime sono buone, prova a restringere la numerosità campionaria… in ogni caso, la vedo molto grigia con co così pochi dati.

    #4849

    blackema
    Participant

    Immaginavo allora che fosse per questo che le cose non andavano bene, purtroppo i dati a disposizione sono pochi e mi devo accontentare di questi perchè sono frutto di una campagna di indagini non ripetibile.
    Vedrò di effettuare altre prove anche se i tempi a disposizione non sono proprio lunghi il che, in mancanza di alternativa dovrò a malincuore(visto il tempo dedicato) virare su altro.

    Ti ringrazio tantissimo per il supporto e il tempo che mi hai dedicato.

Viewing 7 posts - 1 through 7 (of 7 total)

You must be logged in to reply to this topic.