Codice in R

[olaxgabry]

Ciao a tutti,
ho un problema con un codice che ho scritto in R. La compilazione diventa lunghissima: in parte il problema è di R, ma forse si potrebbe rendere il codice migliore. Mi hanno consigliato di usare il C, ma ne so veramente poco.
Intanto posto il codice fatto in R: si tratta del metodo di Eulero per la simulazione delle equazioni differenziali stocastiche.
In breve faccio questo: ho l'equazione differenziale stocastica lineare

$dX_{t}=a*X_{t}*dt + b*X_{t}dW_{t}$

La soluzione è

$X_{t}=X_{0}*e^{(a-\frac{b^{2}}{2})*t+b*W_{t}}$

Per utilizzare Eulero prima discretizzo facendo una partizione dell'intervallo

$0=t_{0}
e suppongo per comodità che l'ampiezza degli intervalli sia pari a $delta$. L'approssimazione con Eulero è

$Y_{n+1}=Y_{n}+a*Y_{n}*delta+b*Y_{n}*\Delta W_{n}$

Quello che faccio nel codice è calcolarmi l'errore all'ultimo istante, facendo ripetere la procedura N volte: in sostanza ottengo un vettore con N componenti (gli errori all'ultimo istante)

## Metodo di Eulero con un N ripetizioni: delta potenza di 0.5 e t0<-0.

Eulero<-function(a,b,x0,delta,N)
{ 
   error<-c(0)                 ## vettore errori
   for(k in 1:N)
   {
     y0<-x0
     n<-1/delta
     MB<-rnorm(n,0,sqrt(delta))            ## genero n numeri casuali
     somma<-0
     y<-c(y0)
     x<-c(x0)
     for(j in 2:(n+1))
     {
       y_approx<-y[j-1]+a*y[j-1]*delta+b*y[j-1]*MB[j-1]                ## soluzione approssimata
       y<-c(y,y_approx)
       t_i<-(j-1)*delta
       somma<-somma+MB[j-1]
       esponente<-(a-(b^2)/2)*t_i+b*somma_MB                        
       x_reale<-x0*exp(esponente)                                               ## soluzione reale
       x<-c(x,x_reale)
     }
     e<-abs(x[n+1]-y[n+1])                ## errore all'ultimo istante
     error<-c(error,e)                         ## aggiorno il vettore degli errori
    }
    error<-error[-1]                           ## elimino la prima componente che ho messo pari a 0
    mean_error<-mean(error)             ## calcolo la media degli errori
    result<-list(errors<-error,mean<-mean_error)
    result
}

Grazie a chiunque voglia aiutarmi.
Ciao

[olaxgabry]

Sembra interessante, la guarderò con calma domani e ti farò sapere.
Come sempre grazie Sergio.

[olaxgabry]

"Sergio":
In ogni caso.... potresti fornire un set di dati in input ed un outout corretto?

In che senso set di input ed output corretti? Scusa l'ignoranza in materia.

[olaxgabry]

"Sergio":Be', la tua funzione ammette come parametri a, b, x0, delta e N e ritorna un vettore.
Per eventualmente giocarci un po' (soprattutto con una eventuale versione in C), comincerei con qualcosa del tipo:

a <- ?
b <- ?
x0 <- ?
delta <- ?
N <- ?

Cosa dovrei mettere al posto dei punti interrogativi?
E se l'eventuale versione in C funzionasse, che vettore dovrei avere come risultato usando quegli input?

Hai perfettamente ragione. Come input ho iniziato con questi:

a<- 1.5
b<- 1
x0<- 1
N<-1000 e 10000 
delta<- 2^(-7) (fino a 2^(-12))

L'output che si ottiene è uno scalare che nel codice è la media. Poi ci ho messo anche il vettore degli errori, ma questo non mi interessa. Con errore intendo la differenza in modulo tra il valore della soluzione reale e quella approssimata all'ultimo istante: calcolo questo per ogni ripetizione N ottenendo un vettore di errori, su questo calcolo la media aritmetica.

[olaxgabry]

Perché sono un "moccolone" e invece di scrivere somma ho messo somma_MB: ho fatto una modifica all'istante, ma quella riga non l'ho corretta.
Scusa la svista, ora ti posto il codice corretto.

## Metodo di Eulero con un N ripetizioni: delta potenza di 0.5 e t0<-0.

Eulero<-function(a,b,x0,delta,N)
{ 
   error<-c(0)
   for(k in 1:N)
   {
     y0<-x0
     n<-1/delta
     MB<-rnorm(n,0,sqrt(delta))
     somma<-0
     y<-c(y0)
     x<-c(x0)
     incr<-c(0)
     for(j in 2:(n+1))
     {
       y_approx<-y[j-1]+a*y[j-1]*delta+b*y[j-1]*MB[j-1]
       y<-c(y,y_approx)
       t_i<-(j-1)*delta
       somma<-somma+MB[j-1]
       esponente<-(a-(b^2)/2)*t_i+b*somma
       x_reale<-x0*exp(esponente)
       x<-c(x,x_reale)
     }
     e<-abs(x[n+1]-y[n+1])
     error<-c(error,e)
    }
    error<-error[-1] 
    mean_error<-mean(error)
    result<-list(errors<-error,mean<-mean_error)
    result
}

L'ho provato e non dovrebbero esserci problemi.

[olaxgabry]

"Sergio":Perdona l'ignoranza.
Nella tua funzione si esegue N volte un ciclo all'inizio del quale ci sono alcune assegnazioni.
A occhio, n <- 1/delta vale per tutti i cicli (n non viene riassegnato) quindi potrebbe essere spostato all'inizio (cioè eseguito una sola volta invece che N).
Vale la stessa cosa per MB<-rnorm(n,0,sqrt(delta))? In altri termini: MB deve essere creato N volte, oppure può bastare crearlo una volta ed usarlo poi N volte lasciandolo sempre uguale?

MB devo essere creato N volte. Hai ragione su

n<- 1/delta

Il discorso è che alla fine ho dovuto creare un ciclo for per fare variare delta, dove quest'ultimo è sempre una potenza di $1/2$: tipo delta tra $[2^(-12),2^(-7)]$. Ti posto quest'ultima cosa in modo da essere più chiaro.
L'output sono i vari delta con l'errore medio associato. I valori h e k sono dove faccio variare delta, cioè $delta\in[2^(-k),2^(-h)]$

## Eulero: ordine di convergenza.


Eulero<-function(a,b,x0,h,k,N)
{ 
   DELTA<-c(0)
   ERRORE<-c(0)
   for(s in h:k)
 {
   delta<-2^(-s)
   DELTA<-c(DELTA,delta)
   error<-c(0)
   for(k in 1:N)
   {
     y0<-x0
     n<-1/delta
     MB<-rnorm(n,0,sqrt(delta))
     somma<-0
     y<-c(y0)                                
     x<-c(x0)
     for(j in 2:(n+1))
     {
       y_approx<-y[j-1]+a*y[j-1]*delta+b*y[j-1]*MB[j-1]
       y<-c(y,y_approx)
       t_i<-(j-1)*delta
       somma<-somma+MB[j-1]
       esponente<-(a-(b^2)/2)*t_i+b*somma
       x_reale<-x0*exp(esponente)
       x<-c(x,x_reale)
     }
     e<-abs(x[n+1]-y[n+1])
     error<-c(error,e)
    }
    error<-error[-1] 
    mean_error<-mean(error)
    ERRORE<-c(ERRORE,mean_error)
 } 
    DELTA<-DELTA[-1]
    ERRORE<-ERRORE[-1]
    result<-list(delta<-DELTA,errori_assoluti<-ERRORE)
    result
}

Grazie per il tempo che mi stai dedicando e scusa se non ti ho parlato subito della variazione di delta.

[olaxgabry]

"Sergio":[quote="olaxgabry"]Grazie per il tempo che mi stai dedicando e scusa se non ti ho parlato subito della variazione di delta.

Figurati! Tra l'altro sono davvero solo pochi minuti perché..... devo studiare analisi [/quote]
Allora ti ringrazio a maggior ragione !
Io nel frattempo sto buttando giù una "traduzione" in java, sperando che le cose migliorino. Ho notato una cosa: quando metto in funzione Task Manager Windows, nella finestra prestazioni noto che l'utilizzo del PC è del 50%. Non è che aumentandone l'utilizzo (non so in che modo) l'esecuzione possa velocizzarsi?
Ciao e buona analisi.

[olaxgabry]

Effettivamente hai ragione (anche se fortunatamente pare non dia problemi nell'esecuzione). Possiamo cambiare il k nella function con l, ovvero

Eulero<-function(a,b,x0,h,l,N) {
...
  for(s in h:l) {
  ...
     for(k in 1:N) {
     ...
     }
  ...
  } 

Sul C immagino che sia come per il java, quindi ci penso subito per attrezzarmi a riguardo e far girare il codice sulla generazione dei numeri pseudocasuali.
A proposito, ormai non ho più parole per ringraziarti soprattutto perché mi aiuti nonostante tu sia sei molto impegnato. Grazie!

[olaxgabry]

Pensi che eclipse possa andare bene anche per il C? Lo uso per java. Dato che ha molti plugin, forse si potranno utilizzare anche per il C. Intanto con il java dei miglioramenti ci sono stati nell'esecuzione, però mi interessa molto il discorso tra C ed R.
Un'opinione: mi è stato consigliato di non utilizzare il java perché non è il massimo per scopi "scientifici", che ne pensi? A me tanto male non sembra, però io sono molto ignorante in materia.

p.s: comunque ti devo fare i complimenti per la preparazione: R, java, C..........da paura!

[olaxgabry]

Domani provo anche io seguendo i tuoi consigli.
Effettivamente con il java un pò di lentezza di esecuzione quando i numeri sono molti c'è, però in 30 minuti il risultato esce mentre con R non se ne vedeva la traccia.
Ci credi che in java avevo dei problemi per una semplice istruzione che permetteva di tener conto solo delle prime 4 cifre decimali: certe cose non le capirò mai !
Ti farò sapere al più presto.
Ciao

[olaxgabry]

"Sergio":Funziona.
Ho scaricato l'installer, ma ho dovuto scegliere "Previous", che scarica e installa la penultima versione, perché un file di quella "Current" sembra corrotto.
Poi ho aggiunto C:\MinGW\bin al PATH, ho aperto una finestra DOS su una directory in cui avevo messo rnorm.c e ho dato il comando:
gcc -lm -o rnorm.exe rnorm.c
Ho poi eseguito rnorm.exe.
Forse è un po' vecchio stile per voi giovani, ma a me l'interfaccia a carattere non dispiace...

Ho problemi di connessione (non sto a casa) e non riesco a scaricare il tutto (queste chiavette internet fanno pensa, se non becchi un posto dove hai campo sei finito). In settimana torno a casa e faccio tutto.
Ciao e scusami per l'attesa.

[olaxgabry]

Mitico!!! Sicuramente mi sarà utile soprattutto perché non ho intenzione di continuare le cose in java soprattutto per i grafici. Sembra veramente un ottimo libro, seguirò il tuo consiglio.
Grazie ancora.