%reg_cwicz
close all; clear all; clc;
%Dane
wspZ=0.05;
Ldemp=30; xmin=0.1; xmax=1.1;
dx=(xmax-xmin)/(Ldemp-1); x=[xmin:dx:xmax]';%rand(Ldemp,1);
[Yemp,Yteor,wspZ,Af]=obiekt(x,-wspZ);
%mamy dane teraz regresja
Kd=10; %wielomian kdego-1 stopnia
%obliczamy macierz wejsc uogulnionych FI
FI(:,1)=ones(Ldemp,1);
for (k=2:Kd) FI(:,k)=x.^(k-1); end
%teraz obliczenia formalne
G=inv(FI'*FI);  %Macierz Gaussa
Aob=G*FI'*Yemp; % Wspolczynniki modelu
Yob=FI*Aob;     %Wartosci Y w obserwacjach
%----liczymy reszty (błąd modelu)
E=Yemp-Yob;
%liczymy wariancje skaldowej losowej i macierz kowariancji wspolczynnikoa, KA
varZ=(E'*E)/(Ldemp-Kd); sigZ=sqrt(varZ)
KA=G*varZ;
%liczymy odchylenie standardowe bledu estymacji dla obserwacji x i licznosc bledow e w 1.sigmowym przedziale
lmalych=0;
for(n=1:Ldemp)
    varYm=FI(n,:)*KA*FI(n,:)'; varYe=varZ+varYm;
    sgYm(n)=sqrt(varYm); sgYe(n)=sqrt(varYe);
    if(abs(E(n))<=sgYe) lmalych = lmalych+1; end;
end
udzialMalych=lmalych/Ldemp*100;
%liczymy odchylenie standardowe bledu estymacji dla zadanych dowolnych xv, np innych niz x
Ldv=500; Xvmin=0; Xvmax=1.55;
dx=(Xvmax-Xvmin)/(Ldv-1); xv=[Xvmin:dx:Xvmax]';%rand(Ldemp,1);
for(n=1:Ldv)
    fi(1)=1; for(k=2:Kd) fi(k)=xv(n)^(k-1); end %obliczamy wejscia uogolnione dla zadanych xv(n)
    Yv(n)=fi*Aob; varYm=fi*KA*fi'; varYe=varZ+varYm;
    sigYm(n)=sqrt(varYm); sigYe(n)=sqrt(varYe);
end
%-----Koniec------
%Generujemy duzo zaklucen, aby zobaczyc jak moga wygladac dane
Ldem=10000;
dx=(Xvmax-Xvmin)/(Ldem-1); xd=[Xvmin:dx:Xvmax]';%rand(Ldemp,1);
[Ye,Yt]=obiekt(xd,wspZ);
figure(1);
subplot(1,2,1);
plot(xd,Ye,'c.', x,Yteor','k',x,Yemp','k.');
xlabel(sprintf('Dane Y_{emp}k i regresja Y_{ob}r: Kd=%d Ldemp=%d sigZ=%.2f (wspZ=%.2f) udzMalych=%.2f %%', Kd,Ldemp,sigZ,wspZ,udzialMalych)); axis('tight');
subplot(1,2,2);
plot(Yob,E.^2,'k.'); xlabel(sprintf('e^2 v.s. Y_{ob}')); axis('tight');
pause();
subplot(1,1,1);
plot(xd,Ye,'c.', x,Yteor','k',x,Yemp','k.',x,Yob,'r'); hold on
plot(xv,Yv,'r--',xv,Yv-sigYm,'b--',xv,Yv+sigYm,'b--',xv,Yv-sigYe,'k--',xv,Yv+sigYe,'k--'); hold off;
xlabel(sprintf('Dane Y_{emp}k i regresja Y_{ob}r: Kd=%d Ldemp=%d sigZ=%.2f (wspZ=%.2f) udzMalych=%.2f %%', Kd,Ldemp,sigZ,wspZ,udzialMalych)); axis('tight');