% T8A5u6 -  lineare Nährung und Extrapolation
   % Um das Neuron besser zu charakterisieren wurde eine Messreihe mit 
   % Stromreizen zwischen 0 und 3.5nA durchgeführt und jeweils 
   % 10 mal die Antwortrate gemessen. Die Reize sind unter [reizT8.mat]
   % die resultierenden Raten unter [ratenT8.mat] abgespeichert.
   % Um unsere Extrapolation zu prüfen wurde die Messreihe erweitert, unter
   % [allereizeT8.mat] und [alleratenT8.mat] sind die entsprechenden Ergebnisse
   % gespeichert.
         
%% Laden Sie Reize und Raten ein und plotten diese als einzelne Datenpunkte.

load reizT8
load ratenT8

plot(reizT8,ratenT8,'.')
xlabel('Rate [nA]')  % Beschriftung der x-Achse
ylabel('Antwortrate') % Beschriftung der y-Achse
xlim([0 3.5])           %Wertebereich x-Achse
ylim([0 110])          %Werteberecih y-Achse
hold on

%% Finden Sie die beste lineare Nährung für die Daten und plotten sie dazu.

data=ratenT8;
para=reizT8;
poly=1;    

[koeff,stru]=polyfit(para,data,poly); 

% angenäherte Gerade 
X=reizT8;     
P=koeff(1)*X + koeff(2);  % p(x)=x+c
plot(X,P)


%% Welche Rate erwarten Sie in Antwort auf einen 5nA starken Reiz? Welche bei -2nA?

x=-2:0.5:5;
p=koeff(1)*x + koeff(2);  % lineare Nährung für den Bereich -2 bis 5nA

sprintf('Bei einer Rate von -2nA erwarte ich einen Antwort von %g und bei einer Rate von 5nA einen Wert von %g.',p(1) ,p(end))

    
%% Plotten Sie alleraten gegen allereize als Einzelpunkte zusammen mit der bereits erstellten linearen Nährung.

load allereizeT8
load alleratenT8

figure
plot(allereizeT8,alleratenT8,'.')  % alle Daten als Einzelpunkte
hold on
plot(x,p)   % bereits erstellte lineare Nährung für den Bereich -2 bis 5nA


%% Erstellen Sie eine neue lineare Nährung für alle Datenpunkte.

data=alleratenT8;
para=allereizeT8;
poly=1;    

[allekoeff,allestru]=polyfit(para,data,poly); 

% angenäherte Gerade für alle Werte mit Polynom 1 (BLAU)
figure
alleX=allereizeT8;     
alleP=allekoeff(1)*alleX + allekoeff(2);  % p(x)=x+c
plot(alleX,alleP)

hold on
plot(allereizeT8,alleratenT8,'.')  % alle Datenpunkt mit anzeigen

%% Wird die Nährung besser, wenn man ein Polynom zweiten Gerades verwendet? 

poly=2;    % poly=2, dadurch wird die Linie gebogener

[allekoeff,allestru]=polyfit(para,data,poly); 

% angenäherte Gerade für alle Werte mit Polynom 2 (ROT)
alleX=allereizeT8;     
allePP=allekoeff(1)*(alleX.^2)+allekoeff(2)*alleX + allekoeff(3);  % p(x)=x^2+x+c
plot(alleX,allePP,'r')


%% Wie sieht sie für Polynome höheren Grades aus? 

% angenäherte Gerde für alle Werte mit Polynom 6 (GRÜN)
poly=4;     

[allekoeff,allestru]=polyfit(para,data,poly); 

alleX=allereizeT8;     
allePPPP=allekoeff(1)*(alleX.^4)+allekoeff(2)*(alleX.^3)+allekoeff(3)*(alleX.^2)+ allekoeff(4)*alleX + allekoeff(5);  
plot(alleX,allePPPP,'g')


%% Verwenden Sie zur Abschätzung der Güte der Nährung  die mittlere quadratische Abweichung der Messpunkte von der geschätzten Kurve.

% für Polynom 1
rP=alleratenT8-alleP ;     % Berechnung der Distanz zwischen den geschätzten und gemessenen Daten 
rsquP=(rP.^2);          % quadratischer vertikaler Abstand
MWrsquP=mean(rsquP)     % Mittelwert der mittleren quadratischen Abweichung aller Messpunkte


% für Polynom 2
rPP=alleratenT8-allePP ;     
rsquPP=(rPP.^2);  
MWrsquPP=mean(rsquPP)

% für Polynom 4
rPPPP=alleratenT8-allePPPP ;     
rsquPPPP=(rPPPP.^2);  
MWrsquPPPP=mean(rsquPPPP)

sprintf('Je höher das Polynom ist desto kleiner wird die mittlere quadratische Abweichung.')


% Eine alternative Möglichkeit zur Fehlerabschätzung ist der von [p,S]=polyfit(x,y,n) zurückgegebene Wert S.normr.
% Dieser Wert gibt den Betrag (euklidische Norm) der Residuen an. (Man kann ihn sich einfach z.B. mit residuen_betrag=S.normr
% anzeigen lassen. Was es mit dieser seltsamen Schreibweise auf sich hat, sehen wir morgen). 



%% Wie viele Koeffizienten werden gebraucht, um die Daten mit polyfit zufriedenstellend anzunähern?
% Wann wird die Schätzung nicht mehr besser, wenn man zusätzliche Koeffizienten hinzunimmt? 

% Polynom 6
poly=6;     

[allekoeff,allestru]=polyfit(para,data,poly); 

alleX=allereizeT8;     
allePPPPPP=allekoeff(1)*(alleX.^6)+allekoeff(2)*(alleX.^5)+allekoeff(3)*(alleX.^4)+allekoeff(4)*(alleX.^3)+allekoeff(5)*(alleX.^2)+ allekoeff(6)*alleX + allekoeff(7);  

rPPPPPP=alleratenT8-allePPPPPP ;     
rsquPPPPPP=(rPPPPPP.^2);  
MWrsquPPPPPP=mean(rsquPPPPPP)

% Polynom 8
poly=8;     

[allekoeff,allestru]=polyfit(para,data,poly); 

alleX=allereizeT8;     
allePPPPPPPP=allekoeff(1)*(alleX.^8)+allekoeff(2)*(alleX.^7)+allekoeff(3)*(alleX.^6)+allekoeff(4)*(alleX.^5)+allekoeff(5)*(alleX.^4)+allekoeff(6)*(alleX.^3)+allekoeff(7)*(alleX.^2)+ allekoeff(8)*alleX + allekoeff(9) ;  

rPPPPPPPP=alleratenT8-allePPPPPPPP ;     
rsquPPPPPPPP=(rPPPPPPPP.^2);  
MWrsquPPPPPPPP=mean(rsquPPPPPPPP)

% zwischen Polynom 6 und 8 verringert sihc die Abweichung nur noch um etwa 0,2