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

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

  echo on
  format compact
  format long
%==============================================================
% SIMULACE casove rady:
% parabola x(t) = b1*t^2 + b2*t +  b3 + wn(t), wn ~ WN(0,sigma^2)
% b1=1, b2=0, b3=2, sigma=0.2
% t=1:0.05:1
%==============================================================
  b1=1; b2=0; b3=2; sigma=0.2;
  t = (0:0.025:1).';
  n = length(t);
  y = b1*t.^2 + b2*t + b3;
  x = y + sigma*randn(n,1);
  tt = linspace(-0.2,1.2).';      % sit bodu pro vykresleni
  yy = b1*tt.^2 + b2*tt + b3;
  plot(tt,yy,'y',t,x,'co'); hold on
pause;
%---------------------------------------------------------------
% DEMONSTRACE UZITI LINEARNIHO REGRESNIHO MODELU:
% Prolozime regresni parabolu:
  p  = 3;
  D  = [t.^2,t,ones(n,1)];       % Matice planu
  DDinv=inv(D'*D);
  b = DDinv*D'*x                 % odhady koeficientu regresni paraboly
  r = x - D*b;                   % rezidua
  s = norm(r)/sqrt(n-p)          % odhad sigma
pause;
%---------------------------------------------------------------
  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
  Ib=[b-be,b+be]                 % intervaly spolehlivosti pro b
  (Ib(2,1)<=0) & (0<=Ib(2,2))    % true => b(2) lze povazovat za nulovy
pause;
%---------------------------------------------------------------
  P  = [tt.^2,tt,ones(size(tt))];% P(k,:) = parametricka funkce pro vypocet
                                 % hodnoty paraboly v bode tt(k)
  xx = P*b;                      % prubeh odhadnute paraboly na [-0.8,1.2]
  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(tt,xx,'r',tt,[xx-xxE1,xx+xxE1],'w--',tt,[xx-xxE2,xx+xxE2],'g:');
  hold off;
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);
  % nyni spoctene koeficienty paraboly a jejich konf. intervaly
  [bE(:,1)-bE(:,2),bE(:,1)+bE(:,2)]-Ib % 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(1)-sE(2),sE(1)+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]-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
