%==============================================================
% ANALYZA SOUBORU uspop.dat:
% Populace USA v desetiletych intervalech
% x=absolutni pocet obyvatel
% t=1,2,...,20: 1790 az 1980
%==============================================================
echo on
format compact
us_pop=[crpath,'/brockwel.dat/uspop.dat'];
[x,xinf]=getdata(us_pop);   % nacteni a vykresleni
n=length(x);
pause;
%---------------------------------------------------------------
xacf=acf(x);
plotdata(xacf); % vypocet a vykresleni autokorelacni funkce
pause;

%---------------------------------------------------------------
% Nalezeni trendu
%---------------------------------------------------------------
% Predpokladame polynomialni trend. Stupen odhadneme pomoci diferenci
% 1. diference:
figure;
dx1=x(2:n)-x(1:n-1); plot(dx1); n1=n-1;
% dx1 neosciluje kolem konstanty => trend neni linearni
pause;

% 2. diference:
dx2=dx1(2:n1)-dx1(1:n1-1); plot(dx2);
% dx2 osciluje kolem konstanty => trend lze povazovat za parabolicky
pause;

% Vypocet regresni paraboly
%----------------------------------------------------------------
% a) primo pomoci 'polyfilt'
%    ----------------------
t=(1:n)'; tt=linspace(1,n).'; N=length(tt);
[y,E,r,s,bE]=polyfilt(2,n,x,t,tt,[],figure); bE
pause;

% b) neprimo pomoci 'polyfilt'
%    ------------------------
[FIR,E0]=polyfilt(2,n,[],t,tt);
yfir=FIR*x; Efir=s*E0;
max(abs(y-yfir))
max(abs(E-Efir))
% Vysledky jsou az na numericke chyby totozne
pause;

% c) primym sestavenim regresniho modelu pro 'linregr'
%    ------------------------------------------------
D=[t.^2 ,t ,ones(n,1)];
C=[tt.^2,tt,ones(N,1)];
[bElm,rE,sE,PbE,CbE]=linregr(x,D,C);
ylm=CbE(:,1); Elm=[CbE(:,2),PbE(1:N,2)];
max(abs(y-ylm))
max(abs(E-Elm))
% Vysledky jsou az na numericke chyby totozne
pause;

% d) Vypocet pomoci funkci MATLABu (toolbox stats)
%    ---------------------------------------------
[b,S] = polyfit(t,x,2);      % b=koeficienty regresni paraboly
[yy,delta]=polyconf(b,tt,S); % yy=regresni prubeh,
                             % delta=polosirka 95% intervalu spolehlivosti
                             %       predikovanych hodnot
max(abs(y-yy))
max(abs(E(:,2)-delta))
% Vysledky jsou az na numericke chyby totozne
pause;

% Standardni funkce MATLABu 'polyval' pocita pouze prubeh zadaneho polynomu:
yp = polyval(b,tt);
max(abs(yy-yp))


%================================== KONEC ==================================
echo off
