// Regresní analýza + validace
// ----------------------------------------------------------
clc, clear, close, mode(0)
getd C:\_Skola\Statistika\Funkce_03-2017\statistika  // natažení funkcí

// Ve dnech x jsme měřili výšku hladiny [m] v řece. 
// Získali jsme data
x=[1 2 3 4 5 6];
y=[2 2.3 2.3 2.2 2.5 2.4];
// Chceme určit:
// a) zda změřené údaje jsou vhodné pro regresní analýzu,
// b) jaký je trend pohybu hladiny (nárůst nebo pokles),
// c) jaká bude hladina v den 10,
// d) kdy hladine dosáhne výšky 2.5 m

// ŘEŠENÍ - regresní přímka
n=length(x);
// průměry
mx=mean(x);
my=mean(y);
// součty čtverců
Sx=sum((x-mx)**2);
Sy=sum((y-my)**2);
Sxy=sum((x-mx).*(y-my));

// koeficienty regresní přímky
disp 'Odhady koeficientů regresní přímky'
b1=Sxy/Sx
b0=my-b1*mx

// predikce
yp=b0+b1*x;
plot(x,y,'.',x,yp,'markersize',12)
legend('data','reg. přímka',2);

// VALIDACE REGRESE
// 1. korelační koeficient a Pearsonův test
disp 'Korelační koeficient'
r=Sxy/sqrt(Sx*Sy)
t=r/sqrt((1-r**2)/(n-2));
q=cdft('PQ',t,n-2);
disp 'Pearsonův test'
pv_t1=2*min(q,1-q)            // oboustranný t-test

// test pomocí funkce
pv_t2=pearson_test(x,y);      // stejné jako pv_t1

// 2. tabulka a Fisherův test
disp 'Fisherův test'
// součty čtverců
Sr=sum((yp-my)**2);
Se=sum((y-yp)**2);
Sc=sum((y-my)**2);            // platí: Sc=Sr+Se
// stupně volnosti
Dr=1;
De=4;                         // počet dat - 2
// průměry čtverců
MSr=Sr/Dr;
MSe=Se/De;
// statistika pro Fisherův test
F=MSr/MSe;
// p-hodnota
pv_F1=1-cdff('PQ',F,Dr,De)    // pravostranný F-test

// Fisherův test (pomocí funkce)
pv_F2=f_test_reg(x,y);        // stejné jako pv_F1

// 3. koeficient detreminace
disp 'Koeficient determinace'
R=r**2

// 4. standardní chyba reziduí   
disp 'Standardní chyba reziduí'
SE=sqrt(Se/(n-2))

// test na nezávislost reziduí     H0: jsou nezávislé
disp 'Test nezávislosti reziduí'
pv_wz=wz_test(y,yp)

// test na autokorelace            výsledek je v  res
disp 'Test na přítomnost autokorelace'
res=autoreg_test(x,y)

// DRUHÝ ZPŮSOB VÝPOČTU ... b=[b0, b1]'
disp '--------------------------------'
disp 'Druhý způsob výpočtu koeficientů'
Y=y(:);
X=[ones(x(:)) x(:)];
b=inv(X'*X)*X'*Y
