% Simulace binarniho procesu a nahodne prochazky (Grossmanova)
echo on
du
p=0.5;
n=100;
t=1:n;
y=rand(1,n);
figure
x=-(y<=0.5)+(y>0.5);
plot(x,'o')
axis([0,100,-2,2])
p=0.6;
x=-(y<=p)+(y>p);
figure
plot(x,'o')
axis([0,100,-2,2])

%nahodna prochazka
x=x'*ones(1,n);
xt=triu(x);
S=sum(xt);
figure
plot(S)

p=0.5;
x=-(y<=p)+(y>p);
S=sum(triu(x'*ones(1,n)));
figure
plot(S)
echo off
% Linearni regresni model
help linregr
                  MULTIPLE LINEAR REGRESSION USING LEAST SQUARES
  for the linear model x = D*b + e, with random disturbances e: N(0,sigma^2*I)
  ----------------------------------------------------------------------------
  CALLING:
  1) REGRESSION RESULT:   [bE,rE,sE,PbE,CbE] = linregr(x ,D,P,C,alpha)
  2) REGRESSION OPERATOR: [R ,E0,1 ,PE0,CE0] = linregr([],D,P,C,alpha)
  INPUT:
  x  ... is an n x 1 vector of observations
         x==[] => R,E0,PE0 and CE0 are matrices defining the regression operator
  D  ... is an n x p full-column rank matrix of regressors
  P  ... is an m x p matrix of regressors for prediction (optional:default=[])
  C  ... is an q x p matrix of parametric functions (optional:default=P)
  alpha ... is a percentual confidence risk (optional: default=[]=5%)
  OUTPUT:
  1) REGRESSION RESULT:
  bE ... =[b,E] is an p x 2 matrix, where
          b are parameter estimates
          [b-E,b+E] are their confidence intervals of (100-alpha) %
  rE ... =[r,E] is an n x 2 matrix, where
          r are residuals x-D*b
          [r-E,r+E] are their confidence intervals of (100-alpha) %
          (!!! E is not implemented in this version => rE is of size p x 1)
  sE ... =[s,E1,E2] is an 1 x 3 vector, where
          s is the estimate of the standard deviation sigma
          [s-E1,s+E2] is its confidence interval of (100-alpha) %
  PbE ...=[Pb,E] is an (m+1) x 2 matrix, where
          Pb(1:m,:) are predicted values of x = P*b + e
          Pb(m+1,:) is  predicted value  of y = sum(P)*b + sum(e)
          [Pb-E,Pb+E] are their confidence intervals of (100-alpha) %
          PbE = [] if P was not supplied on input (P=[])
  CbE ...=[Cb,E] is an q x 2 matrix, where
          Cb are estimates of the parametric functions C*b
          [Cb-E,Cb+E] are their confidence intervals of (100-alpha) %
          CbE = [] if neither C nor P was supplied on input (C=P=[])
  2) REGRESSION OPERATOR:
  R  ... p x n matrix defining the regression operator b = R*x
  E0 ... is an p x 1 matrix defining the confidence vector E from bE
         by E=sigma.*E0 where sigma is an estimate of the standard deviation
         of random disturbances, usually sigma=norm(D*b-x)/sqrt(n-p).
  sE ... sE = s = 1 ... unit standard deviation (scalar)
  PE0... is an (m+1) x 1 matrix determining the confidence matrix E from PbE
         by E=sigma.*PE0 where sigma=norm(D*b-x)/sqrt(n-p) - see E0
  CE0... is an q x 1 matrix determining the confidence matrix E from CbE
         by E=sigma.*CE0 where sigma=norm(D*b-x)/sqrt(n-p) - see E0

% Demonstrace polynomialni regrese pro odhad trendu v 'airpass.dat'
crcv2
  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)
n =
   144
  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
ans =
     1     1     1
     4     2     1
     9     3     1
    16     4     1
    25     5     1
    36     6     1
    49     7     1
    64     8     1
    81     9     1
   100    10     1
  more off
pause;
%---------------------------------------------------------------
  DDinv=inv(D'*D)
DDinv =
   0.00000000290780  -0.00000042163123   0.00001025969336
  -0.00000042163123   0.00006515549845  -0.00177903082833
   0.00001025969336  -0.00177903082833   0.06426877113497
  b = DDinv*D'*x                 % Koeficienty regresni paraboly
b =
  1.0e+002 *
   0.00007008198318
   0.01640995151950
   1.12380037509439
pause;
%---------------------------------------------------------------
  r = x - D*b;                   % rezidua
  s = norm(r)/sqrt(n-p)          % odhad sigma
s =
  44.90560667970558
  alpha = 0.05                   % riziko
alpha =
   0.05000000000000
  tp = tinv(1-alpha/2,n-p)       % kvantil Studentova rozlozeni
tp =
   1.97693148863424
  d  = sqrt(diag(DDinv))
d =
   0.00005392403558
   0.00807189559224
   0.25351286187286
  be = tp*s*d                    % polosirka bodoveho konfidencniho pasu pro b
be =
   0.00478712285965
   0.71658501622562
  22.50568235957112
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)
ans =
    5.684341886080802e-014
  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;
return
%---------------------------------------------------------------
% 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
bE =
  1.0e+002 *
   0.00007008198318   0.00004787122860
   0.01640995151950   0.00716585016226
   1.12380037509439   0.22505682359571
  [b,be]    % porovname s drive spoctenymi hodnotami
ans =
  1.0e+002 *
   0.00007008198318   0.00004787122860
   0.01640995151950   0.00716585016226
   1.12380037509439   0.22505682359571
pause;
%---------------------------------------------------------------
  max(abs(rE(:,1)-r))  % maximalni odchylka mezi nyni a drive spoctenymi
ans =
    3.410605131648481e-013
                       % rezidui
pause;
%---------------------------------------------------------------
  sE(1)         % nyni spocteny odhad sigma porovname
ans =
  44.90560667970557
  s             % s drive spoctenou hodnotou.
s =
  44.90560667970558
  % navic dostavame pro s levou a pravou polosirku konfidencniho intervalu:
  sE(2),sE(3)
ans =
   4.68542453707688
ans =
   5.93077926860596
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs(CbE-[xx,xxE1]))
ans =
  1.0e-012 *
   0.45474735088646   0.07105427357601
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho polosirkou predikcniho konfidencniho intervalu.
  max(abs(PbE(1:length(tt),:)-[xx,xxE2]))
ans =
  1.0e-012 *
   0.45474735088646   0.04263256414561
pause;
%---------------------------------------------------------------
% b) neprimy vypocet (pomoci operatoru):
%    ---------------------------------
  [R,E0,sE,PE0,CE0] = linregr([],D,P); % = linregr([],D,P,P,alpha);
  sE
sE =
     1
  % 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
bE =
  1.0e+002 *
   0.00007008198318   0.00004787122860
   0.01640995151950   0.00716585016226
   1.12380037509439   0.22505682359571
  [b,be]          % porovname s drive spoctenymi hodnotami
ans =
  1.0e+002 *
   0.00007008198318   0.00004787122860
   0.01640995151950   0.00716585016226
   1.12380037509439   0.22505682359571
pause;
%---------------------------------------------------------------
  rr=x-D*bb;
  max(abs(rr-r))  % maximalni odchylka mezi nyni a drive spoctenymi rezidui
ans =
    2.842170943040401e-013
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs([P*bb,s*CE0]-[xx,xxE1]))
ans =
  1.0e-012 *
   0.28421709430404   0.07105427357601
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]))
ans =
  1.0e-012 *
   0.28421709430404   0.02842170943040

%================================== KONEC ==================================
echo off
% Demonstrace polynomialni regrese na simulovanych datech
crcv2d
  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
b =
   0.49599060591642
   0.50902445135145
   1.93783180997219
  r = x - D*b;                   % rezidua
  s = norm(r)/sqrt(n-p)          % odhad sigma
s =
   0.19173273888458
pause;
%---------------------------------------------------------------
  alpha = 0.05                   % riziko
alpha =
   0.05000000000000
  tp = tinv(1-alpha/2,n-p)       % kvantil Studentova rozlozeni
tp =
   2.02439416391196
  d  = sqrt(diag(DDinv))
d =
   1.99729923048355
   2.06590170644431
   0.44656083447096
  be = tp*s*d                    % polosirka bodoveho konfidencniho pasu pro b
be =
   0.77523699145384
   0.80186453742111
   0.17332930015327
  Ib=[b-be,b+be]                 % intervaly spolehlivosti pro b
Ib =
  -0.27924638553742   1.27122759737026
  -0.29284008606966   1.31088898877256
   1.76450250981891   2.11116111012546
  (Ib(2,1)<=0) & (0<=Ib(2,2))    % true => b(2) lze povazovat za nulovy
ans =
     1
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;
return
%---------------------------------------------------------------
% 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
ans =
  1.0e-013 *
   0.08548717289614   0.11990408665952
  -0.21982415887578  -0.19317880628478
  -0.05551115123126  -0.05329070518201
pause;
%---------------------------------------------------------------
  max(abs(rE(:,1)-r))  % maximalni odchylka mezi nyni a drive spoctenymi
ans =
    1.598721155460225e-014
                       % rezidui
pause;
%---------------------------------------------------------------
  sE(1)         % nyni spocteny odhad sigma porovname
ans =
   0.19173273888458
  s             % s drive spoctenou hodnotou.
s =
   0.19173273888458
  % navic dostavame pro s levou a pravou polosirku konfidencniho intervalu:
  [sE(1)-sE(2),sE(1)+sE(3)]
ans =
   0.15669279845827   0.24710098048152
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs(CbE-[xx,xxE1]))
ans =
  1.0e-013 *
   0.15987211554602   0.01221245327088
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho polosirkou predikcniho konfidencniho intervalu.
  max(abs(PbE(1:length(tt),:)-[xx,xxE2]))
ans =
  1.0e-013 *
   0.15987211554602   0.00777156117238
pause;
%---------------------------------------------------------------
% b) neprimy vypocet (pomoci operatoru):
%    ---------------------------------
  [R,E0,sE,PE0,CE0] = linregr([],D,P); % = linregr([],D,P,P,alpha);
  sE
sE =
     1
  % 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
bE =
   0.49599060591643   0.77523699145384
   0.50902445135143   0.80186453742111
   1.93783180997218   0.17332930015327
  [b,be]-bE       % porovname s drive spoctenymi hodnotami
ans =
  1.0e-013 *
  -0.09992007221626  -0.01998401444325
   0.22648549702353  -0.01443289932013
   0.05329070518201   0.00055511151231
pause;
%---------------------------------------------------------------
  rr=x-D*bb;
  max(abs(rr-r))  % maximalni odchylka mezi nyni a drive spoctenymi rezidui
ans =
    1.776356839400251e-014
pause;
%---------------------------------------------------------------
  % maximalni odchylky mezi nyni a drive spoctenym regresnim prubehem a
  % jeho konfidencni polosirkou.
  max(abs([P*bb,s*CE0]-[xx,xxE1]))
ans =
  1.0e-013 *
   0.18207657603853   0.01276756478319
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]))
ans =
  1.0e-013 *
   0.18207657603853   0.00888178419700

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