Accelerare simulazione dell'algoritmo movimento Levy

Question

Ecco il mio piccolo script per la simulazione di Levy movimento:

clear all; 
clc; close all; 
t = 0; T = 1000; I = T-t; 
dT = T/I; t = 0:dT:T; tau = T/I; 
alpha = 1.5; 
sigma = dT^(1/alpha); 
mu = 0; beta = 0; 
N = 1000; 
X = zeros(N, length(I)); 
for k=1:N 
    L = zeros(1,I); 
    for i = 1:I-1 
     L((i + 1) * tau) = L(i*tau) + stable2(alpha, beta, sigma, mu, 1); 
    end 
    X(k,1:length(L)) = L; 
end 

q = 0.1:0.1:0.9; 
quant = qlines2(X, q, t(1:length(X)), tau); 
hold all 
for i = 1:length(quant) 
    plot(t, quant(i) * t.^(1/alpha), ':k'); 
end 

Dove stable2 restituisce un stable random variable con determinati parametri (è possibile sostituirlo con normrnd(mu, sigma) per questo caso, non è fondamentale); qlines2 restituisce i quantili necessari per il tracciamento.

Ma non voglio parlare di matematica qui. Il mio problema è che questa implementazione è piuttosto lenta e vorrei accelerarla. Sfortunatamente, l'informatica non è il mio campo principale - ho sentito parlare di metodi come la memoizzazione, la vettorizzazione e che ci sono molte altre tecniche, ma non so come usarle.
Ad esempio, sono abbastanza sicuro che dovrei sostituire questo sporco doppio per-loop in qualche modo, ma non sono sicuro di cosa fare invece.
EDIT: Forse dovrei usare (e imparare ...) un altro linguaggio (Python, C, qualsiasi altro funzionale)? Ho sempre pensato che Matlab/OCTAVE sia progettato per il calcolo numerico, ma se il cambiamento, allora per quale?

Answer 1

Il bit cruciale è, come hai detto, il ciclo for, Matlab non gli piace, quindi la vettorizzazione è davvero la parola chiave. (Insieme con preallocare lo spazio.

ho solo cambiato per la sezione ciclo un po 'in modo che non c'è bisogno di ripristinare L più e più volte, invece salviamo tutti L s in una matrice più grande (anche io elimiated il comando length(L)).

L = zeros(N,I); 
for k=1:N 
    for i = 1:I-1 
     L(k,(i + 1) * tau) = L(k,i*tau) + normrnd(mu, sigma); 
    end 
    X(k,1:I) = L(k,1:I); 
end 

Ora si può già vedere che X(k,1:I) = L(k,1:I); nel loop è obsoleto e questo significa anche che siamo in grado di cambiare l'ordine dei cicli. Questo è fondamentale, perché i i -steps sono ricorsivi (dipende dalla passaggio precedente) che significa che non possiamo vettorializzare questo ciclo, possiamo solo vettorizzare lo k -loop.

Ora il vostro codice originale aveva bisogno di 9,3 secondi sulla mia macchina, il nuovo codice deve ancora circa lo stesso tempo)

L = zeros(N,I); 
for i = 1:I-1 
    for k=1:N 
     L(k,(i + 1) * tau) = L(k,i*tau) + normrnd(mu, sigma); 
    end 
end 
X = L; 

Ma ora siamo in grado di applicare la vettorializzazione, invece di looping cittadi tutte le righe (il loop over k) possiamo invece eliminare questo ciclo e fare tutte le righe in "una volta".

L = zeros(N,I); 
for i = 1:I-1 
    L(:,(i + 1) * tau) = L(:,i*tau) + normrnd(mu, sigma); %<- this is not yet what you want, see comment below 
end 
X = L; 

Questo codice richiede solo 0,045 secondi sulla mia macchina. Spero che tu abbia ancora lo stesso risultato, perché non ho idea di cosa stai calcolando, ma spero anche che tu possa vedere come si fa a vettorizzare il codice.

PS: Ho appena notato che ora usiamo lo stesso numero casuale nell'ultimo esempio per l'intera colonna, questo ovviamente non è quello che vuoi. Instad dovresti generare un intero vettore di numeri casuali, ad esempio:

L = zeros(N,I); 
for i = 1:I-1 
    L(:,(i + 1) * tau) = L(:,i*tau) + normrnd(mu, sigma,N,1); 
end 
X = L; 

PPS: Ottima domanda!