Navigation

ACHTUNG! NEUE, VERÄNDERTE VERSION! 

Tag 8: Modellbildung

A) Kurvenanpassung

Kurvenanpassung dient dazu, Messdaten durch eine Kurve - also einen mathematisch beschreibbare Funktion - zu beschreiben. Die Idee dabei ist, dass die Messdaten systematisch von einem vorgegebenen Parameter abhängen, aber statistisch um die zu erwartenden Werte streuen.

Graphisch trägt man die Werte des Parameters (z.B. Zeit, Düngemittelmenge, Sonnenscheindauer etc) auf der x-Achse auf, die Messwerte und die theoretischen Werte in Abhängigkeit von diesem Parameter auf der y-Achse. Häufig werden dabei die Messwerte als Punkte (oder andere Symbole) dargestellt, die theoretischen Werte als durchgezogene Kurve.

Diese Darstellung hilft dabei, die Idee der Interpolation zu verstehen:  Durch sie lässt sich abschätzen, welche Messwerte zu erwarten gewesen wären, wenn man zwischen den gewählten Parameterwerten weitere Messpunkte ermittelt hätte.

Außerdem ermöglicht Kurvenanpassung auch Extrapolation, die Vorhersage, wie sich die Messwerte außerhalb des getesteten Parameterbereichs verhalten.

Um eine Funktion zu finden, die die Abhängigkeit der vorhandenen Messdaten vom variierten Parameter möglichst gut beschreibt, möchte man eine Kurve finden, die möglichst nah an den Messdaten liegt. Hierfür (oder generell wenn man versucht, Messdaten durch eine Bildungsregel (also ein Modell) zu beschreiben), braucht man ein Maß dafür, wie gut das Modell die Daten widerspiegelt. Dadurch kann man verschiedene Modelle miteinander vergleichen und sich für das am besten passende entscheiden. Ein Standardverfahren, um die optimale Funktion (bzw. das optimale Modell) zu bestimmen, ist die mittlere quadratische Abweichung zwischen den Messwerten und den geschätzten Werten zu minimieren. (Prinzipiell gibt es aber noch viele andere Möglichkeiten, die Güte eines Modells zu messen.)

Hierzu wird zunächst für jeden Datenpunkt das Residuum bestimmt, also die Differenz zwischen dem empirisch ermittelten Wert yi und der Schätzung für diesen Wert (Schätzwerte macht man gerne mit einem Dach kenntlich):

Bild

Die mittlere quadratische Abweichung zwischen allen Datenpunkten und der Schätzung wird ermittelt als Summe der quadrierten Residuen, geteilt durch die Anzahl der Datenpunkte n:

 Bild

Was bewirkt das Quadrat?


 In diesem Kurs benutzen wir Polynome der Form

p(x)=p1xn+p2xn-1+...pnx+pn+1

 um den empirischen Daten eine Funktion anzupassen.

Wie sieht in dieser Darstellung eine Gerade mit Steigung 2.5 und y-Achsenabschnitt 7 aus?

Natürlich lassen sich nicht alle Abhängigkeiten durch Polynome adäquat darstellen. Welche Gegenbeispiele fallen Ihnen ein?

In Matlab stehen für die Kurvenanpassung eines Polynoms an Messdaten die Befehle polyfit und polyval zur Verfügung:
[p,S]=polyfit(x,y,n) gibt im Vektor p=[p1 p2... pn pn+1] diejenigen Koeffizienten für die Gleichung

p(x)=p1xn+p2xn-1+...pnx+pn+1

zurück, die den quadratischen Fehler zwischen Messwerten und geschätzten Werten minimiert.
S ist eine Struktur, die polyval als Eingabe braucht (Strukturen behandeln wir morgen ausfuehrlicher).
y_schaetz = polyval(p,x,S) bekommt als Eingabe den von polyfit zurückgegebenen Vektor p und die Struktur S, sowie einen Vektor von x-Werten. Für diese x -Werte berechnet polyval die gemäß der gefundenen Gleichung vorhergesagten y-Werte (die x Werte werden in die Gleichung eingesetzt).
Versieht man die Funktion polyval mit einem zweiten Ausgabeparameter [y_schaetz,delta] = polyval(p,x,S) wird im Vektor delta ein Maß dafür angegeben, wie sicher diese Schätzung ist. Im Bereich y_schaetz-delta bis y_schaetz+delta liegen 50% der Messpunkte.

Trotz dieser Befehle möchte ich Sie bitten, zunächst einmal für einen einfachen Fall Kurvenanpassung “von Hand” vorzunehmen, damit Sie die Grundidee besser verstehen:

T8A1) Als Vorbereitung beschaeftigen wir uns zuerst mal mit dem Fehlermaß, der mittleren quadratischen Abweichung: Schreiben Sie eine Funktion, die einen Vektor gemessener und einen gleichlangen Vektor geschaetzter y-Werte bekommt und den mittleren quadratische Fehler zurueckgibt. 

*T8A2) Erweitern Sie Ihre Funktion so, dass Sie auch fuer eine Matrix von Messwerten (bei der jedem x-Wert mehrere y-Werte zugeordnet sind) und einen Vektor der geschaetzten y-Werte (fuer jeden x-Wert  einer) den mittleren quadratischen Fehler berechnet. 

T8A3) Jetzt zur eigentlichen Kurvenanpassung:
Wenn man Nervenzellen mit Strom reizt, reagieren sie darauf mit einer Folge von Aktionspotentialen.
Stellen Sie sich ein Neuron vor, das auf einen Stromimpuls von 1 nA mit einer Frequenz von 32 Sp/s reagiert, bei 2 nA produziert es 60 Sp/s.
a) Welche Antwortfrequenz erwarten Sie bei einer Strominjektion von 1.5 nA?
b) Wie lautet die lineare Gleichung fuer die Gerade, die beide Punkte verbindet?
c) Plotten Sie die Gerade und die beiden Datenpunkte

T8A4) Sie fuehren jetzt eine weitere Messung durch und stellen fest, dass bei 1.5 nA Strominjektion 44 Sp/s erzeugt werden.
a) Erweitern Sie Ihren plot um den neuen Datenpunkt.
b) Berechnen Sie den quadratischen vertikalen Abstand des Datenpunkts von der Geraden.
c) Schaetzen Sie "per Auge" eine Gerade, die alle drei Punkte gut annaehert und plotten sie zusammen mit den Datenpunkten.
d) Berechnen Sie Ihren mittleren quadratischen Fehler.
f) Benutzen Sie jetzt den Befehl polyfit, um die beste lineare Naeherung zu finden und polyval um den Fehler zu berechnen und die Gerade zu plotten.
[p,S]=polyfit(x,y,n) gibt im Vektor p=[p1 p2... pn pn+1] diejenigen Koeffizienten fuer die Gleichung

p(x)=p1xn+p2xn-1+...pnx+pn+1

zurueck, die den quadratischen Fehler zwischen Messwerten und geschaetzten Werten minimiert.
S ist eine Struktur, die polyval als Eingabe braucht (Strukturen behandeln wir morgen ausfuehrlicher).
y_schaetz = polyval(p,x,S) bekommt als Eingabe den von polyfit zurueckgegebenen Vektor p und die Struktur S, sowie einen Vektor von x-Werten. Fuer diese x-Werte berechnet polyval die gemaess der gefundenen Gleichung vorhergesagten y-Werte.
Versieht man die Funktion polyval mit einem zweiten Ausgabeparameter [y_schaetz,delta] = polyval(p,x,S) wird im Vektor delta mit angegeben, wie sicher diese Schaetzung ist. Im Bereich y_schaetz-delta bis y_schaetz+delta liegen 50% der Messpunkte.
Plotten Sie die von Matlab berechnete Kurve und deren delta-Werte als Fehlerbalken gemeinsam mit den Datenpunkten in eine Abbildung.
Wie gut war Ihre Schaetzung?
*g) Eigentlich braucht man die Funktion polyval nicht unbedingt fuer unsere Anwendung. Schreiben Sie eine Funktion, die Vektor p der Koeffizienten und den x-Vektor als Eingangsargumente bekommt und den resultierenden Vektor der Schaetzwerte zurueckgibt. Testen Sie, ob ihre Schaetzungen denjenigen der Funktion polyval entsprechen.

T8A5) 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.
a) Laden Sie Reize und Raten ein und plotten diese als einzelne Datenpunkte.
b) Finden Sie die beste lineare Nährung für die Daten und plotten sie dazu.
c) Welche Rate erwarten Sie in Antwort auf einen 5nA starken Reiz? Welche
bei -2nA?

T8A6) Um unsere Extrapolation zu prüfen wurde die Messreihe erweitert, unter
[allereizeT8.mat] und [alleratenT8.mat] sind die entsprechenden Ergebnisse
gespeichert.
a) Plotten Sie alleraten gegen allereize als Einzelpunkte zusammen mit der bereits erstellten linearen Nährung für den in T8A3) bearbeiteten Teil der Daten.
b) Erstellen Sie eine neue lineare Nährung für alle Datenpunkte.
c) Wird die Nährung besser, wenn man ein Polynom zweiten Gerades verwendet?
Wie sieht sie für Polynome höheren Grades aus? Verwenden Sie zur Abschätzung der Güte der Nährung  die mittlere quadratische Abweichung der Messpunkte von der geschätzten Kurve.
d) 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).
e) 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? Untersuchen Sie systematisch, wie sich die Anpassung mit wachsender Koeffizientenzahl verbessert.

T8A7) Berechnen Sie für jeden der Stromwerte, mit dem Messungen durchgeführt wurden, Mittelwert und Standardabweichung und plotten diesen mit errorbar in eine Abbildung. (Hinweis: für diese Aufgabe ist es günstiger, aus dem langen Vektor der Messwerte eine Matrix zu machen. Das geht mit reshape.)

*T8A8) Mathematisch Interessierte könnten sich unter dem Stichwort "least squares fitting" ansehen, wie die Datenanpassung in Matlab vorgenommen wird und welche Möglichkeiten die Curve Fitting Toolbox noch bietet.



B) Variation von Modellparametern

Das Ziel vieler Experimente ist es, eine Abhängigkeit zwischen Parameterwerten und einer Messgröße zu finden. Man versucht, diese Abhängigkeit mathematisch durch eine Formel zu beschreiben. Dadurch kann man für gegebene Parameterwerte vorhersagen, welche Messwerte man erwartet. Dieses Vorgehen bezeichnet man auch als Simulation. Im Experiment kann dann überprüft werden ob die Vorhersagen den tatsächlichen Messungen entsprechen. Wenn dies nicht der Fall ist, muss das Modell (also die Formel zur Beschreibung der Abhängigkeit) entsprechend angepasst werden.

T8B1) Tatsächlich handelt es sich bei den gerade verwendeten Daten nicht um echte Messungen, sondern um mit Rauschen überlagerte Punkte auf einer Sigmoidfunktion.
Sigmodidfunktionen spielen in der Biologie eine wichtige Rolle, da es kaum Beispiele für unbegrenzt lineare Abhängigkeiten gibt.
Sigmoidfunktionen kann man mit der Formel
y=a./(1+(exp(b*(c-x))))
darstellen.
a) Erzeugen Sie einen Vektor x, der von -10 bis 10 in kleinen Schritten läuft.
Setzen Sie die Parameter a=120; b=1; c=2. Berechnen Sie y entsprechend der Sigmoidformel.
b) Plotten Sie die Sigmoidfunktion zusammen mit alleraten und Ihrem besten Fit.
Wie ähnlich sind sich die beiden Linien?
*c) Mit welcher Standardabweichung habe ich die Datenpunkte erzeugt?

T8B2) Es lohnt sich, die Sigmoidfunktion genauer kennenzulernen. Schreiben Sie ein Skript, in dem Sie systematisch die drei Parameter der Sigmoidfunktion variieren und die Auswirkung grafisch als Kurvenscharen darstellen. Was bewirken die drei Parameter? (Tipp: gehen Sie von den Werten a=1, b=1, c=1 aus, um die Effekte gut zu sehen.)

**T8B3) Schreiben Sie ein Programm, das für eine gegebene Menge an Messwerten eine Sigmoidfunktion bestimmt, die den mittleren quadratischen Abstand zu den Messwerten minimiert.
Testen Sie diese Funktion, indem Sie künstliche Daten herstellen, bei denen eine Sigmoidfunktion mit  Rauschen überlagert wird.


D) Hausaufgaben

T8H1) Sie impfen 20 Petrischalen jeweils mit 100 Einzellern und zählen 10 Stunden lang 10 mal pro Stunde,  wie viele Einzeller sich in der Petrischale befinden. Dadurch erhalten Sie folgende Matrix [populationen.mat]. Schauen Sie sich die Messpunkte grafisch an. Ermitteln Sie das Polynom, das das Wachstum der Population am besten beschreibt.

* T8H2) Fünf Ratten leben in einem Käfig. Bei unserer Beobachtung interessieren wir uns im Moment nur für die Frage wie viele der Ratten zu einem gegebenen Zeitpunkt schlafen und wie viele wach sind.

a) Schreiben Sie eine Simulation des Schlaf-Wach-Verhaltens der Ratten, in der sie jeweils einmal in der Stunde bestimmen, wie viele Ratten wach sind, basierend auf folgenden Regeln:

•    Ausgangszustand: 3 Ratten sind wach, 2 schlafen
•    Wenn zu einem gegebenen Zeitpunkt nur 0 oder 1 Ratte wach ist, ist die Wahrscheinlichkeit für jede der Ratten 80%, dass sie bei der nächsten Messung schläft.
•    Wenn 2 oder mehr Ratten wach sind, ist die Wahrscheinlichkeit für jede Ratte nur 20%, dass sie bei der nächsten Messung schläft (denn die Tiere stören sich gegenseitig.)

b) Lassen Sie Ihre Simulation über einen langen Zeitraum laufen und protokollieren Sie die Anzahl der wachen Ratten. Wie groß ist der Prozentsatz der Zeit, den die Tiere schlafend verbringen?

**c) Jetzt machen wir die Sache mal etwas komplizierter: Die fünf Ratten sind in Einzelkäfigen untergebracht, die nebeneinander stehen. Jetzt stören sich nur noch diejenigen Ratten gegenseitig, die direkt nebeneinander untergebracht sind. Nun sehen die Regeln so aus:

•    Ausgangszustand: Ratte 1, 3 und 5 sind wach, 2 und 4 schlafen
•    Wenn alle Nachbarn schlafen, ist die Wahrscheinlichkeit für jede Ratte 80%, dass sie bei der nächsten Messung schläft.
•    Wenn mindestens ein Nachbar wach ist, ist die Wahrscheinlichkeit für jede Ratte nur 20%, dass sie bei der nächsten Messung schläft

Lassen Sie Ihre Simulation über einen langen Zeitraum laufen und protokollieren Sie, welche der Ratten wann wach ist.
Wie groß ist der Prozentsatz der Zeit, den die einzelnen Tiere schlafend verbringen?
Gibt es Unterschiede zwischen den Tieren? Sind diese signifikant?
Stellen Sie den Verlauf des Schlaf-Wach-Musters der Tiere grafisch dar.


 Zum 9. Kurstag



Webmaster:l4 Jun0ieqtta Kre9atzbery1kg (jutta.kretnczber1dhg@uofcl.dez21) (Stand: 07.11.2019)