% primo test di applicazione del test di Anderson % lo script fa uso del "System Identification Toolbox": prima viene % generata una sequenza di dati ottenuti applicando ad un processo % stocastico stazionario ARX(2,2). I dati raccolti vengono % utilizzati poi per identificare un modello di tipo ARX, i cui ordini sono % assegnati dall'utente. % La bont? del risultato viene valutata applicando il test di Anderson ai % residui di predizione del predittore ad 1 passo identificato in % precedenza. % % Provare ad identificare modelli di complessit? certamente sbagliata per % verificare se il test viene superato o meno. Provare anche con gli ordini % corretti, sempre per verificare se il test viene superato o meno. % Testato in Matlab R2012a % % corso "Sistemi Dinamici" 2019/2020 % G. Fenu % clc clear close all %------- Creazione dei dati ------ ndati=20000; y=zeros(ndati,1); y(1)=randn;y(2)=randn; e=randn(size(y)); u=2*randn(size(y)); %----- il processo VERO: ARX(2,2) for i=3:length(y) y(i)=1.2*y(i-1)-0.32*y(i-2)+u(i-1)+0.5*u(i-2)+e(i); end figure;subplot(3,1,1); plot(y);title('uscita misurata y(k)'); subplot(3,1,2); plot(u,'r');title('ingresso u(k)'); subplot(3,1,3); plot(e,'g');title('rumore e(k)'); %------ Ripartizione dati per validazione -------- while 1 traindata=round(input('Inserire percentuale di dati per l''addestramento (0 per terminare): ')*ndati/100); if traindata==0 break end ytrain=y(1:traindata); utrain=u(1:traindata); ytest=y(traindata+1:end); utest=u(traindata+1:end); datitrain=iddata(ytrain,utrain); datitest=iddata(ytest,utest); na=input('Ordine modello rispetto ad y? (n_a): '); nb=input('Ordine modello rispetto ad u? (n_b): '); % identifico un ARX model = arx(datitrain,[na nb 1]); %----- Validazione ------ prederr=pe(model,datitest); err=prederr.OutputData; figure;subplot(2,1,1); plot(ytest,'x-r');hold on;plot(ytest-err,'ob'); legend('uscita y(k)','uscita predetta y_p(k)'); subplot(2,1,2);plot(err,'k');title('errore di predizione'); % -= Verifica Gauss =- passointegrale=1e-6; beta=1.64; %1.96 % 0.05;%2.58; % 0.99 intervallo=-beta:passointegrale:beta; puntigaussiana=1/sqrt(2*pi)*exp(-intervallo.^2/2); area=trapz(intervallo,puntigaussiana); alfa=1-area; % --- Test di Anderson --- N=length(err); taumax=floor(N/5); gamma0=1/N*sum(err.^2); gamma=zeros(taumax,1); for j=1:taumax somma=0; for k=1:N-j somma=somma+err(k)*err(k+j); end gamma(j)=somma/N; end rho=gamma/gamma0; %>< nout=length(find(abs(rho)>beta/sqrt(N))); if ((nout/taumax)