% Data la seguente serie storica di 23 elementi % y=[48.6 49.1 49.3 47.8 47 46.9 48.9 48.7 48.9 52.7 50.8 ... % 50.6 52.4 55 57.8 60.1 59.8 61.9 57.5 58.1 55.9 54.6 54.5]'; % Calcolare tramite implementazione manuale % 1) La serie delle medie mobili di 3 termini % % 2) La serie delle medie mobili semplici di 5 termini % % Fare il plot della serie storica originale, della serie delle medie % mobili di lunghezza 3 e 5. Aggiungere la legenda delle 3 linee al % grafico. Commentare il grafico % % 3) La serie delle medie semplici di lunghezza k (con k numero intero % positivo dispari) % % 4) La serie storica delle medie ponderate di lunghezza k (k numero % positivo e dispari) i pesi decrescono in maniera lineare e sommano ad % uno. % % 5) La serie storiche delle medie semplici di lunghezza k % in modo da ridurre la dimensione della finestra vicino ai punti finali % dell'input per includere solo gli elementi esistenti. % % 6) La serie storica delle medie mobili con il metodo del livellamento % esponenziale % % 7) Utilizzando la funzione movmean verificare che è possibile ottenere % le medie ottenute ai punti 3) e 5) % % 8) Calcolare la media mobile semplice calcolata tramite la funzione % movavg ed analizzare la differenza con l'applicazione della funzione % movmean % % 9) Calcolare la serie delle medie mobili con il metodo del livellamento % esponenziale con la funzione movavg. % % 10) Controllare l'uguaglianza delle implementazioni manuali con quelle % automatiche. %% Serie storica iniziale y=[48.6 49.1 49.3 47.8 47 46.9 48.9 48.7 48.9 52.7 50.8 ... 50.6 52.4 55 57.8 60.1 59.8 61.9 57.5 58.1 55.9 54.6 54.5]'; T=length(y); %% 1) Implementazione manuale media mobile semplice di lunghezza 3 % Creo un vettore colonna di nan di lunghezza T mm3=nan(T,1); for i=2:T-1 mm3(i)=mean(y(i-1:i+1)); end %% 2) Implementazione manuale media mobile semplice di lunghezza 5 mm5=nan(T,1); for i=3:T-2 mm5(i)=mean(y(i-2:i+2)); end %% Grafico % Fare il plot della serie storica originale, della serie delle medie % mobili di lunghezza 3 e 5. Commentare il grafico close all plot(y,'LineWidth',2) hold("on") plot(mm3,'LineWidth',2,'LineStyle','--') plot(mm5,'LineWidth',2,'LineStyle',':') title('y MM 3 e 5 termini') legend({'y' 'MM 3 termini' 'MM 5 termini'},'Location','best') % Commento la media mobile di 5 termini è molto più smooth di quella di 3 % termini. La media mobile di breve "svolta" (cambia direzione) molto % prima della media di lungo periodo. Ad esempio il minimo di mm3 è al % tempo t=5 al contrario il minimo di mm5 è al tempo t=6; % Similmente, il max di mm5 è in corrispondenza di t=18 al contrario il max % di mm3 è in corrispodenza di t=17 %% Implementazione manuale media mobile semplice di lunghezza k (k numero positivo e dispari) k=5; % supponiamo che k sia sempre un numero dispari z=(k-1)/2; mmkFill=nan(T,1); for i=z+1:T-z mmkFill(i)=mean(y(i-z:i+z)); end %% 4) Implementazione manuale media ponderata di lunghezza k (k numero positivo e dispari) % i pesi decrescono in maniera lineare k=5; % supponiamo che k sia sempre un numero dispari % Pesi a somma 1 w=((1:k)/(0.5*k*(k+1)))'; % Osservazione: 0.5*k*(k+1) = sum(1:k)= somma dei primi k numeri naturali z=(k-1)/2; mmkweighted=nan(T,1); for i=z+1:T-z mmkweighted(i)=sum(y(i-z:i+z).*w); end % Dato che una FAQ era di spiegare da dove vengono questi pesi, aggiungo di % seguito una breve spiegazione % I pesi decrescono in maniera lineare e sommano ad % uno. Di conseguenza: % se ci sono 3 termini l'ultimo viene pesato 3, il penultimo viene pesato % 2, il primo viene pesato 1. La somma dei pesi deve fare 1 di conseguenza % divide per la somma ossia per 1+2+3=3*(1+3)/2; % % se ci sono 5 termini l'ultimo viene pesato 5, il penultimo viene pesato % 4, ......., il primo viene pesato 1. La somma dei pesi deve fare 1 di % conseguenza divide per la somma ossia per 1+2+3+4+5=5*(1+5)/2 % % ..... % % se ci sono k termini l'ultimo viene pesato k, il penultimo viene pesato % k-1, ......., il primo viene pesato 1. La somma dei pesi deve fare 1 di % conseguenza divide per la somma ossia per 1+2+3+4+5+ ...+k =k(1+k)/2 % % Spero che ora la linea di codice % % w=((1:k)/(0.5*k*(k+1)))'; % % sia chiara %% MEDIE MOBILI SHRINK % 5) Implementazione manuale media mobile semplice di lunghezza k (k numero positivo e dispari) % in modo da ridurre la dimensione della finestra vicino ai punti finali % dell'input per includere solo gli elementi esistenti. k=5; % supponiamo che k sia sempre un numero dispari z=(k-1)/2; mmkShrink=nan(T,1); for i=1:T iminusz=max([1, i-z]); iplusz=min([i+z,T]); mmkShrink(i)=mean(y(iminusz:iplusz)); end %% Implementazione manuale media mobile esponenziale % 6) La serie storica delle medie mobili con il metodo del livellamento % esponenziale mme=zeros(length(y),1); mme(1)=y(1); alpha=2/(k+1); for i=2:T mme(i)=(1-alpha)*mme(i-1)+alpha*y(i); end %% 7) Implementazione automatica k=5; % Con l'opzione discard le prime due e le ultime due osservazioni vengono % rimosse. La serie delle medie mobili presenta lunghezza T-s+1 auto1=movmean(y,k,'Endpoints','discard'); % Con l'opzione fill le prime (k-1)/2 e le ultime (k-1)/2 osservazioni % vengono riempite con NaN mmkFillAUTO=movmean(y,k,'Endpoints','fill'); % Con l'opzione shrink si riduce la dimensione della finestra vicino ai punti finali % dell'input per includere solo gli elementi esistenti. mmkShrinkAUTO=movmean(y,k,'Endpoints','shrink'); %% Media mobile semplice calcolata tramite la funzione movavg mmkShrinkAUTO1=movavg(y,'simple',k); % Se si utilizza movavg l'ultima media mobile viene centrata al tempo T, % la penultima al tempo T-1, .... disp('Confronto tra mm calcolata tramite movmean e movavg') aa=[y mmkShrinkAUTO mmkShrinkAUTO1]; % Mostro le prime 8 righe head(table(y,mmkShrinkAUTO,mmkShrinkAUTO1)) % mmkShrinkAUTO1(1)= y(1) % mmkShrinkAUTO1(2)= mean(y(1:2)) % mmkShrinkAUTO1(3)= mean(y(1:3)) % .... % mmkShrinkAUTO1(k)= mean(y(1:k)) % mmkShrinkAUTO1(k+1)= mean(y(2:k+1)) ............. % mmkShrinkAUTO(1) è uguale mmkShrinkAUTO1(1+(k-1)/2); % Ad esempio, se k=5 si ha che mmkShrinkAUTO(1) è uguale a % mmkShrinkAUTO1(3) %% Media mobile esponenziale implementazione automatica % 9) Calcolare la serie delle medie mobili con il metodot del livellamento % esponenziale con la funzione movavg. mmeAUTO= movavg(y,'exponential',k); %% Controllo l'uguaglianza dell'implementazione manuale con quella automatica. % % 10) Controllare l'uguaglianza delle implementazioni manuali con quelle % automatiche. diff1=max(abs(mmkFill-mmkFillAUTO),[],"omitnan"); diff2=max(abs(mmkShrink-mmkShrinkAUTO)); diff3=max(abs(mme-mmeAUTO)); tol=1e-12; assert(diff1