%% Interaktyvus pratimas nr. 11. Statistinė duomenų analizė
% Šiame interaktyviame pratime pasistengsime suprasti kai kuriuos aspektus iš statistinės duomenų analizės. Pamatysime kaip galime nustatyti priklausomybę tarp dviejų kintamųjų, kaip sukurti likimybinį modelį prognozei atlikti bei kaip grafiškai įvertinti gauto modelio kokybę.
%% I dalis. Duomenų generavimas. Pirsono koreliacijos koeficientas
% Nagrinėsime žaislinį pavyzdį. Pradžioje sugeneruosime 1000 stebėjimų, pasiskirsčiusių pagal standartinį normalųjį skirstinį N\left(0;1\right), t.y. tokį, kurio vidurkis m=0, o dispersija \sigma^2 =1. Tai bus kintamasis x. Kintamajam y sukurti kurso klausytojai gali parinkti tris koeficientus a,b,c. MATLAB apskaičiuos y=a\;x^2 +b\;x+c\left(1+\varepsilon \right), t.y. vektorių sudaro 1000 stebėjimų, o \varepsilon yra paklaida, kurios dydį taip pat galite pasirinkti intervale \left\lbrack 0;10\right\rbrack. MATLAB taip pat pavaizduos visus gautus taškus \left(x_i ,y_i \right), čia i=1,2,\ldotp \ldotp \ldotp ,1000.
% Galite pasirinkti koeficientus a,b,c\in \left\lbrack -5;5\right\rbrack ir paklaidos koeficientą \varepsilon \in \left\lbrack 0;10\right\rbrack. Žemiau vaizduojama sugeneruotų duomenų sklaidos diagrama.

x = randn([1 1000]);

a = 2.2;
b = -2;
c = 3;

epsilon = 0.39;

err = 1 + epsilon * randn([1 1000]);

y = a*x.^2 + b*x + c * err;

figure
scatter(x, y)
hold on

% Vienas iš svarbiausių statistinių rodiklių yra Pirsono koreliacijos koeficientas, kuris parodo ar kintamieji x ir y yra susieti tiesiniu ryšiu. Pagal nusistovėjusias tradicijas tiesinio ryšio stiprumas remiantis Pirsono koeficientu yra
% labai stiprus, jeigu 0\ldotp 9\le |\rho \left(X,Y\right)|\le 1;
% stiprus, jegu 0\ldotp 7\le \left|\rho \left(X,Y\right)\right|<0\ldotp 9;
% vidutinis, jeigu 0\ldotp 5\le \left|\rho \left(X,Y\right)\right|<0\ldotp 7;
% silpnas, jeigu 0\ldotp 3\le \left|\rho \left(X,Y\right)\right|<0\ldotp 5;
% labai silpnas (arba neegzistuojantis), jeigu \left|\rho \left(X,Y\right)\right|<0\ldotp 3.
% Reikia pastebėti, kad koreliacijos koeficientas gali nurodyti tik tiesinio ryšio stiprumą, tačiau negali vertinti ryšio priežastingumo.
% Nustatykite mažą paklaidą bei fiksuote koeficientus b ir c. Keiskite koeficientą a ir stebėkite, kaip keičiasi Pirsono koreliacijos koeficientas. Taip pat pabandykite fiksuoti koeficientus ir stebėkite, kaip koreliacijos koeficientas priklauso nuo paklaidos.

rho = corrcoef(x, y);
disp(['Pirsono koreliacijos koeficientas tarp kintamųjų x ir y yra: ' num2str(rho(1,2), 3)])

%
%% II dalis. Tiesinė regresija
% Šioje dalyje atliekama tiesinė regresija, t.y. ieškoma tiesės lygtis, kuri sieja kintamuosius x ir y. MATLAB  apskaičiuoja tiesės y=b\;x+c koefcientų b ir c taškinius įverčius taip, kad mažiausių kvadratų metodo prasme modelio paklaida būtų mažiausia. Iš Tikimybių Teorijos ir Statistikos kurso yra žinomos šių įverčių išraiškos:
% \begin{array}{l}
% \hat{b} =\rho \left(X,Y\right)\cdot \frac{{\bar{\;s} }_y }{{\bar{\;s} }_x };\\
% \hat{\;c} =\bar{\;y} -\hat{\;b} \cdot \bar{\;x} ,
% \end{array}
% čia  \rho \left(X,Y\right) yra Pirsono koreliacijos koeficientas, o {\bar{\;s} }_x ir {\bar{\;s} }_y yra empiriniai pataisytieji standartiniai nuokrypiai.
% Tikimybinio modelio tikslumas gali būti vertinamas grafiškai. Kuo arčiau nubrėžtos tiesės yra pavaizduoti duomenys, tuo modelis yra geresnis. Siūlome išbandyti kelis variantus (pagal tiesinio ryšio stiprumą) ir įvertinti, kaip keičiasi modelio adekvatumas.
% Žemiau nubraižytas gautos tiesės grafikas, kad galetumėte grafiškai įvertinti, ar apskaičiuotas tikimybinis modelis tinka sugeneruotiems duomenims. Taip pat palyginimui pateiksime jūsų pasirinktus koeficientus ir MATLAB paketo apskaičiuotus įverčius.

coeffs = polyfit(x,y,1);
lr = polyval(coeffs, x);
plot(x, lr);
disp(['Jūsų pasirinkti koeficientai buvo: a = ' num2str(a) ', b = ' num2str(b) ', c = ' num2str(c)]);
disp(['Koeficientų taškiniai įverčiai buvo: b = ' num2str(coeffs(1)) ', c = ' num2str(coeffs(2))]);

%
%% III dalis. Interpoliacijos uždavinys
% Šioje dalyje nagrinėsime interpoliacijos uždavinį, t.y. bandysime priklausomybę tarp kintamųjų x ir y aprašyti daugianariu. Paprastumo dėlei apsiribosime kvadratiniu daugianariu. Taigi, šioje interaktyvaus pratimo dalyje MATLAB  apskaičiuoja parabolės y=a\;x^2 \;+\;b\;x+c koefcientų a,\;b ir c taškinius įverčius taip, kad mažiausių kvadratų metodo prasme modelio paklaida būtų mažiausia. Tam MATLAB naudoja Vandermondo matricas, apie kurias kalbėjome paskaitoje apie matricų taikymus. Praktikoje tuo atveju, kai sukurto modelio kokybė yra prasta, galima spręsti kelis interpoliacijos uždavinius keičiant daugianario eilę. Tokiu būdu analitikas gali pasirinkti optimaliausią variantą.
% Šio kurso klausytojus siūlome įvertinti abiejų modelių adekvatumą ir pasirinkti geresnį modelį. tam tikslui žemiau bus nubraižytas gautų tiesės ir parabolės grafikai, kad galetumėte grafiškai palyginti modelių kokybę. Taip pat palyginimui pateiksime jūsų pasirinktus koeficientus ir MATLAB paketo apskaičiuotus įverčius parabolės koeficientams.

sqfit = polyfit(x, y, 2);
t = min(x):0.001:max(x);
sq = polyval(sqfit, t);

plot(t, sq)
hold off
disp(['Jūsų pasirinkti koeficientai buvo: a = ' num2str(a) ', b = ' num2str(b) ', c = ' num2str(c)]);
disp(['Koeficientų taškiniai įverčiai buvo: a = ' num2str(sqfit(1)) ', b = ' num2str(sqfit(2)) ', c = ' num2str(sqfit(3))]);

%% Pabandykite patys
% Pabandykite paeksperimentuoti su didesnės eilės daugianariais, t.y. pakeiskite trečios dalies kodą taip, kad galima būtų pasirinkti interpoliuojamo daugianario eilę.
% Įvertinkite grafiškai gauto modelio adekvatumą.