%==============================================================
%                    LINEARNI REGRESNI MODEL
%==============================================================

% (C) 1996-98, V. Vesely, Masaryk University of Brno (Czech Republic)

  echo on
  format compact
  format long
%==============================================================
% ANALYZA SOUBORU airpass.dat:
% Cestujici v mezinarodni letecke doprave
% x=mesicni souhrny v tisicich cestujicich
% t=1,2,...,144: leden 1949 az prosinec 1960
%==============================================================
  air_pass=[crpath,'/brockwel.dat/airpass.dat'];
  [x,xinf]=getdata(air_pass);   % nacteni a vykresleni
pause;
%---------------------------------------------------------------
% DEMONSTRACE UZITI LINEARNIHO REGRESNIHO MODELU:
% Prolozime regresni parabolu:
  n  = length(x)
  p  = 3;
  t  = (1:n).';                  % Sloupcovy vektor nezavisle promenne
  more on
  D  = [t.^2,t,ones(n,1)];       % Matice planu
  D(1:10,:)                      % vypis jejich prvych 10 radku
  more off
pause;
%---------------------------------------------------------------
  DDinv=inv(D'*D)
  b = DDinv*D'*x                 % Koeficienty regresni paraboly
pause;
%---------------------------------------------------------------
  r = x - D*b;                   % rezidua
  s = norm(r)/sqrt(n-p)          % odhad sigma
  alpha = 0.05                   % riziko
  tp = tinv(1-alpha/2,n-p)       % kvantil Studentova rozlozeni
  d  = sqrt(diag(DDinv))
  be = tp*s*d                    % polosirka bodoveho konfidencniho pasu pro b
pause;
%---------------------------------------------------------------
  tt = linspace(1,n+20).';       % sit bodu pro vykresleni
  xx = polyval(b,tt);            % prubeh paraboly na teto siti
  P  = [tt.^2,tt,ones(size(tt))];% P(k,:) = parametricka funkce pro vypocet
                                 % hodnoty paraboly v bode tt(k)
  max(abs(xx-P*b))               % P*b je totez jako polyval(b,tt)
  gamma = sqrt(diag(P*DDinv*P'));
  xxE1  = tp*s*gamma;            % xxE1(k) polosirka bodoveho konfidencniho
                                 % intervalu pro hodnotu regresni paraboly v
                                 % bode tt(k).
pause;
%---------------------------------------------------------------
  gamma = sqrt(1+diag(P*DDinv*P'));
  xxE2   = tp*s*gamma;           % xxE2(k) polosirka bodoveho konfidencniho
                                 % intervalu predikovane hodnoty v tt(k).
  plot(t,x,'y',tt,xx,'w',tt,[xx-xxE1,xx+xxE1],'w--',tt,[xx-xxE2,xx+xxE2],'w:');
pause;
%---------------------------------------------------------------
% Vyse uvedeny postup je zimplementovan v ponekud dokonalejsi podobe
% v procedure 'linregr'.
% Nastudujte nyni pomoci 'help linregr' pouzivani teto procedury,
% pripadne provedte jeji uplny vypis pomoci 'type linregr' a prostudujte
% podrobneji cely algoritmus.
keyboard;
%---------------------------------------------------------------
% Zopakujeme nyni celou parabolickou regresi s vyuzitim procedury
% 'linregr'.
% Pro ilustraci uzijeme oba mozne zpusoby:
% a) primy vypocet:
%    -------------
  [bE,rE,sE,PbE,CbE] = linregr(x,D,P); % = linregr(x,D,P,P,alpha);
  bE        % nyni spoctene koeficienty paraboly a jejich konf. polosirky
  [b,be]    % porovname s drive spoctenymi hodnotami
pause;
%---------------------------------------------------------------
  max(abs(rE(:,1)-r))  % maximalni odchylka mezi nyni a drive spoctenymi
                       % rezidui
pause;
%---------------------------------------------------------------
  sE(1)         % nyni spocteny odhad sigma porovname
  s             % s drive spoctenou hodnotou.
  % navic dostavame pro s levou a pravou polosirku konfidencniho intervalu:
  sE(2),sE(3)
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs(CbE-[xx,xxE1]))
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho polosirkou predikcniho konfidencniho intervalu.
  max(abs(PbE(1:length(tt),:)-[xx,xxE2]))
pause;
%---------------------------------------------------------------
% b) neprimy vypocet (pomoci operatoru):
%    ---------------------------------
  [R,E0,sE,PE0,CE0] = linregr([],D,P); % = linregr([],D,P,P,alpha);
  sE
  % R=DDinv*D'=inv(D'*D)*D', ostatni vystup odpovida sigma=sE=1;
  % Pozadovane hodnoty pro konkretni x pak tedy dopocitame takto:
  bb=R*x;
  bE=[bb,s*E0]    % koeficienty paraboly a jejich konf. polosirky
  [b,be]          % porovname s drive spoctenymi hodnotami
pause;
%---------------------------------------------------------------
  rr=x-D*bb;
  max(abs(rr-r))  % maximalni odchylka mezi nyni a drive spoctenymi rezidui
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs([P*bb,s*CE0]-[xx,xxE1]))
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho polosirkou predikcniho konfidencniho intervalu.
  max(abs([P*bb,s*PE0(1:length(tt))]-[xx,xxE2]))

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