Navigation

Tag 6: Zufallszahlen

A) Zufallszahlen, Wahrscheinlichkeitsverteilungen, Mittelwert und Standardabweichung

Mit Hilfe von experimentellen Messungen versucht man, allgemeingültige Aussagen zu treffen und Regeln für untersuchte Zusammenhänge aufzustellen. Man variiert einen Parameter (z.B. Menge an Dünger) und beobachtet den dadurch hervorgerufenen Effekt auf eine Messgröße. Dies wäre sehr einfach, wenn grundsätzlich jede Beobachtung immer gleich ausfallen würde, wenn man sie mehrfach wiederholt. In der Realität ist dies allerdings nicht der Fall: Messdaten hängen grundsätzlich zumindest in einem bestimmten Rahmen vom Zufall ab, denn in einem Experimente können niemals alle Zufallsfaktoren ausgeschlossen werden. (z.B. könnte es bei einer Studie über die Wirksamkeit eines Medikaments einen Einfluss haben, wie viel die Patienten geraucht haben oder ob sie gerade Stress hatten.) Auch wenn es eine eindeutige Abhängigkeit zwischen dem variierten Parameter und der gemessenen Größe gibt, werden die Messwerte unterschiedlich ausfallen, sie streuen um den erwarteten Wert.

Wiederholt man ein Experiment in genau gleicher Weise sehr häufig und schaut sich die dabei gewonnen Messwerte an, ergibt sich eine Häufigkeitsverteilung. Diese gibt an, wie häufig ein bestimmter Messwert beobachtet wurde und dient dazu, die Wahrscheinlichkeit dieses Messwertes abzuschätzen.

Für Messdaten (auch künstlich vom Computer erzeugte Zufallszahlen) werden Häufigkeitsverteilungen empirisch ermittelt, um dadurch auf die zugrundeliegende Wahrscheinlichkeitsfunktion zu schließen. Dafür benutzt man Histogramme. Diese teilen den gesamten Wertebereich der Variablen in mehrere Bereiche auf. Das Histogramm gibt für jeden der Bereiche an, wie häufig der Wert der Variable in einer Messung innerhalb des jeweiligen Bereichs lag. Dazu wird für jeden Teilbereich (sogenannte Klassen) ein Rechteck dargestellt, dessen Fläche die gemessene Häufigkeit repräsentiert. In Matlab werden Histogramme mit dem Befehl hist erzeugt. Dieser kann vielfältig eingesetzt werden:

 hist(v)   Teilt den Wertebereich des Vektors v in 10 gleich grosse Klassen ein. Wenn hist ohne Ausgabeargument aufgerufen wird, stellt es die Häufigkeit des Auftretens der Klassen als Balkengrafik dar.
 h=hist(v) Wenn hist mit einem Ausgabeargument aufgerufen wird, produziert es keine grafische Ausgabe, sondern gibt den Vektor der Häufigkeiten zurück. (Kann mit mehreren Eingabeargumenten kombiniert werden.)
 hist(v,nbins) Teilt den Wertebereich des Vektors v in nbins gleich grosse Klassen ein. (Mit oder ohne Ausgabeargument verwendbar)
hist(v,centers) Benutzt den Vektor centers als Mittelpunkte der Klassen, in die die Elemente von v aufgeteilt werden. Wenn hist ohne Ausgabeargument aufgerufen wird, stellt es die Häufigkeit des Auftretens der Klassen als Balkengrafik dar. (Mit oder ohne Ausgabeargument verwendbar)
 [h,xout]=hist(v) Wenn hist mit zwei Ausgabeargumenten aufgerufen wird, produziert es keine grafische Ausgabe, sondern gibt zwei Vektoren zurück: den Vektor h der Häufigkeiten und den Vektor xout der Klassenmittelpunkte. (Kann mit mehreren Eingabeargumenten kombiniert werden

Die meisten biologischen Daten lassen sich durch eine Normalverteilung (auch Gauß-Verteilung genannt)  beschreiben, bei der Messwerte umso häufiger auftreten, je näher sie am Erwartungswert, dem Mittelwert der Verteilung liegen.

Für eine normalverteilte Zufallsvariable x entspricht die Wahrscheinlichkeitsdichte folgender Formel:

Bild

 

wobei μ den Mittelwert und σ die Standardabweichung der Wahrscheinlichkeitsverteilung angibt.

Für normalverteilte Messwerte (oder auch mit einem Zufallsgenerator erzeugte Zufallszahlen) kennt man diese Kennwerte der den Daten zugrunde liegenden Normalverteilung nicht. Man kann sie jedoch näherungsweise aus den Messwerten x1 bis xn ermitteln:

Empirischer Mittelwert:

Bild

Empirische Standardabweichung:

Bild

 

Auch wenn es keine schlechte Übung ist, diese Formeln einmal in Matlab umzusetzen, kann man stattdessen auch einfach die Befehle mean und std benutzen.

Achtung: Die Berechnung von empirischem Mittelwert und empirischer Standardabweichung macht ausschließlich für NORMALVERTEILTE Werte Sinn!

Eine andere wichtige Verteilung, die in diesem Kurs auch betrachtet wird, ist die Gleichverteilung, bei der alle Werte in einem bestimmten Bereich mit gleicher Wahrscheinlichkeit auftreten. (Für gleichverteilte Werte ist es absolut sinnlos, Mittelwert und Standardabweichung zu berechnen.)

Bevor wir uns mit der statistischen Auswertung von echten Messdaten beschäftigen, erzeugen wir zunächst einmal selber "Messdaten" mit Matlab, nämlich Zufallszahlen. Diese werden beispielsweise gebraucht, wenn man Experimente plant, in denen bestimmte Reize in zufälliger Reihenfolge präsentieren werden sollen. Eine weitere wichtige Anwendung von Zufallszahlen sind Simulationen biologischer Prozesse. Wenn man Zufallszahlen künstlich erzeugt, ist (im Gegensatz zur Auswertung von Messdaten) die Wahrscheinlichkeitsverteilung bekannt (also die Wahrscheinlichkeit dafür, dass eine Zufallsvariable einen bestimmten Wert annimmt).

Im Rahmen des Kurses erzeugen wir normalverteilte und gleichverteilte Zufallszahlen, wofür es in Matlab die Funktionen randn und rand gibt. Eine weitere Funktion, um Zufallszahlen zu erzeugen, ist die Funktion randperm. v=randperm(n) liefert einen Vektor der ganzen Zahlen von 1 bis n in zufälliger Reihenfolge.

Aufgaben:

T6A1) Zufallszahlen werden in Matlab mit zwei Funktionen erzeugt: rand und randn. Beide Funktionen bekommen als Eingaben die Dimensionen der Matrix mit Zufallszahlen, die sie erzeugen und zurückliefern sollen (M=randn(7,10)).
Erzeugen Sie sich einige Beispiele der Zufallszahlen: Was passiert, wenn man die Funktionen mehrfach hintereinander in gleicher Weise aufruft?
In welchen Bereich liegen die Werte für die beiden Funktionen?
Was könnte der Unterschied zwischen den beiden Funktionen sein?

T6A2)  Erzeugen Sie sich jeweils einen sehr langen Vektor (z.B. 10000 Punkte) mit jeder der beiden Funktionen rand und randn.

  • Schauen Sie sich die jeweilige Verteilung der Zufallszahlen mit dem Befehl hist an.
  • Was sind die Unterschiede zwischen den beiden Verteilungen? 
  • Mit hist(v,n) teilt hist den Vektor v in n gleich große Bereiche ein. Sehen Sie sich die Verteilungen für verschiedene Werte von n an. 
  • Schätzen Sie aus der Abbildung ab: Was ist der Mittelwert? Was die Standardabweichung?
  • Berechnen Sie Mittelwerte, Standardabweichungen, Varianzen, Minima und Maxima Ihrer beiden Vektoren mit mean, std, var, min und max.

T6A3) Modifizieren Sie Ihre Zufallsvektoren, indem Sie diese

  • mit verschiedenen Faktoren multiplizieren
  • verschiedene Zahlen hinzuaddieren
  • Wie wirken sich diese Änderungen auf die Verteilungen aus?
  • Wie wirken sie sich auf Mittelwert, Standardabweichung, Minimum und Maximum aus?

T6A4) Laden Sie den Vektor mit den Anzahlen an Sonnenblumenkernen von 100 Blumen [sbkerne.mat]. Berechnen Sie Mittelwert, Varianz und Standardabweichung. Wie habe ich diesen Vektor erzeugt? (Nein, ich habe mich nicht in den Garten gesetzt und gezählt...)

*T6A5) Setzen Sie die oben angegebenen Formel für die Wahrscheinlichkeitsdichte einer Normalverteilung in Matlab um. Schreiben Sie eine Funktion, die für die Eingabeargumente x (ein Vektor, der den Definitionsbereich angibt, z.B. x=-4:0.01:4), Mittelwert und Standardabweichung die entsprechende Wahrscheinlichkeitsdichteverteilung als Vektor zurück gibt und als Kurve grafisch darstellt (bitte mit Titel und beschrifteten Achsen). Variieren Sie die Parameter Mittelwert und Standardabweichung. Wie verändern diese die Kurve?

T6A5) Schreiben Sie eine Funktion wuerfel, die Ihnen jeweils eine ganze Zahl zwischen 1 und 6 zurückgibt.

T6A5) Benutzen Sie diese Funktion in einer weiteren Funktion wuerfel_verteilung, die als Eingabewert bekommt, wie oft gewürfelt wird, und als Ausgabe die Verteilung (als in einem Vektor gespeichertes Histogramm) der erzielten Würfelergebnisse zurückliefert.

B) Datenaufnahme

Die meisten biologischen Daten ändern sich nicht schrittweise, sondern kontinuierlich (z.B. Ph-Wert, Körpertemperatur, Membranpotential einer Nervenzelle etc.). Wenn man solche Daten aufzeichnet, steht man vor dem Problem, dass eine kontinuierliche Aufzeichnung grundsätzlich nicht möglich ist. Spätestens wenn ein Computer benutzt wird, um die Daten zu speichern, müssen sie diskretisiert werden -  es ist nur in bestimmten Zeitschritten eine Aufzeichnung möglich und durch die begrenzte (wenn auch sehr hohe) Genauigkeit der Darstellung von Nachkommazahlen sind auch die Messwerte nicht zu 100% kontinuierlich messbar. Grundsätzlich gilt: je größer die Abtastrate ist (also je mehr Messungen pro Zeiteinheit vorgenommen werden), desto genauer ist die Messung - aber auch desto "kostspieliger", denn es fallen mehr Daten an, die gespeichert und verarbeitet werden müssen (und häufig wird die Messhardware deutlich teurer, wenn man hohe Abtastraten erreichen will.)

T6C1) Lassen Sie eine Parabel zeichnen:  x=-5:1:5, y=x.^2, plot(x,y)
Was fällt Ihnen an dieser Parabel auf? Machen Sie die einzelnen Punkte sichtbar.
Machen Sie die Schritte kleiner und plotten Sie die Parabel erneut.
Erstellen Sie einen gemeinsamen plot, bei dem die oben erzeugten Werte mit einzelnen Punkten und die Werte für x_fein=-5:0.01:5 mit einer durchgezogenen Linie in eine Abbildung gezeichnet werden.

T6C2) Wenn eine kontinuierliche Größe gemessen wird, ist die Wahl der Abtastrate wichtig - zu geringe Abtastraten können zu grundlegend falschen Ergebnissen führen. Dieser Effekt soll mit folgendem Skript verdeutlicht werden:
[aliasing_effekt.m]

T6C3) Schauen Sie sich den Effekt einmal für echte Daten an: Im Rahmen eines Elektrophysiologie-Praktikums werden Intrazellulärmessungen des Membranpotentials einer Nervenzelle des Blutegels vorgenommen. Die Zelle wird mit einer zeitlichen Auflösung von 10000 Punkten/Sekunde (10kHz) mit elektrischem Strom gereizt. Der Zeitverlauf dieses Stroms (in nA) ist in [stimulus.mat] gespeichert. Die Spannungsantwort (in mV) wird ebenfalls mit 10kHz aufgezeichnet. Sie ist in [spikes.mat] gespeichert. (Die schnellen, stereotypen Änderungen der Membranspannung, mit denen Nervenzellen kommunizieren, nennt man Aktionspotentiale oder auf Englisch spikes.)

  • Laden Sie die Dateien ein und plotten Sie sie übereinander in zwei subplots einer Abbildung.
  • Passen Sie die Zeitachsen in den Plots so an, dass Sekunden dargestellt werden und beschriften Sie die Achsen.
  • Erzeugen Sie einen Vektor, in den nur jeder 10. Punkt des Vektors spikes übernommen wird (dies entspricht einer Aufnahme mit 1kHz).
  • Erzeugen Sie einen weiteren Vektor, in den nur jeder 100. Punkt des Vektors spikes übernommen wird (dies entspricht einer Aufnahme mit 100Hz).
  • Öffnen Sie ein neues Graphikfenster, und legen Sie darin 3 untereinander liegende subplot-Fenster an.
  • Zeichnen Sie jeweils einen der drei Vektoren mit den Aufnahmen bei 10kHz, 1kHz und 100Hz in die drei Fenster ein.
  • Wie unterscheiden sich die Zeitverläufe?
  • Benutzen Sie die Lupen-Funktion und die Darstellung einzelner Punkte auf dem Graphen, um sich den Zeitverlauf eines einzelnen Aktionspotentials bei den drei Vektoren anzusehen.
  • Schauen Sie sich auch einen Zeitausschnitt an, bei dem "nichts" passiert, also der Messwert um einen festen Wert fluktuiert.
  • Berechnen Sie für alle drei Vektoren Mittelwert und Standardabweichung der ersten 200ms der Messung.

C) Sparse Matrizen:

Für große zweidimensionale Matrizen, bei denen ein hoher Prozentsatz der Elemente Nullen sind, bietet Matlab einen speziellen Datentyp an, die Sparse Matrizen. Diese benötigen wesentlich weniger Speicherplatz als normale Matrizen bestehend aus dem Datentyp double, bei denen jede 0 mit 32bit Speicherplatz kodiert wird. Bei sparse Matrizen werden nur diejenigen Elemente, die ungleich 0 sind, gemeinsam mit ihren Zeilen- und Spalten-Indizes in der Matrix abgespeichert. 

Eine Sparse Matrix wird erzeugt mit dem Befehl sparse, z.B.
A=[0 0 0; 0 9.5 0; 1 0 0];
S1=sparse(A)
ergibt:
S1 =
(3,1) 1.0000
(2,2) 9.5000

Mit dem Befehl full wird eine sparse Matrix in eine normale zweidimensionale Matrix des Typs double übersetzt: F=full(S1) ergibt die gleiche Matrix, die oben als A definiert wurde. 

Mit sparse Matrizen lassen sich alle Matrix-Manipulationen (z.B. Indizierung, Anfügen von Zeilen oder Spalten etc) und die meisten Operationen mit der gleichen Syntax ausführen, wie bei den normalen Matrizen (z.B. e=S1(2,2); S1=[S1; 0 0 7]; S2=2*S1).

Nützliche Funktionen für den Umgang mit sparse Matrizen:

 issparse(s) Test, ob eine Matrix sparse ist
find
finde Elemente ungleich 0 (wie sonst auch)
 nonzeros Werte der Elemente ungleich 0
 speye generiert sparse Identitätsmatrix
 spfun wendet Funktion auf Elemente ungleich 0 an.
 sprand(S) erzeugt eine sparse Matrix mit der gleichen Struktur wie bei Matrix S aber gleichverteilten zufälligen Elementen
 sprandn(S)

erzeugt eine sparse Matrix mit der gleichen Struktur wie bei Matrix S aber normalverteilten zufälligen Elementen

 spones
Ersetzt die Elemente ungleich 0 mit einsen
 spy Visualisierung des sparse-Musters

T8C1) Erzeugen Sie eine übersichtliche Matrix, in der viele Nullen enthalten sind. 

  • Erzeugen Sie unter einem anderen Namen eine sparse Matrix aus ihrer ursprünglichen Matrix.
  • Probieren Sie aus, ob Sie mit dem Befehl full die ursprüngliche Matrix herausbekommen.
  • Führen Sie einige Matrix-Operationen sowohl auf der sparse- als auch auf der full-Version der Matrix aus.
  • Überprüfen Sie, ob am Ende beide Matrizen den gleichen Inhalt darstellen.

T8C2) Mit sparse Matrizen lässt sich im Arbeitsspeicher und beim Abspeichern in Dateien Platz sparen. Jedoch benutzt Matlab beim Speichern in Dateien ohnehin eine geschickte Art der Komprimierung, so dass sich hier der Unterschied zwischen sparse und normaler Matrix nicht so stark auswirkt:

  • Matlab Erzeugen Sie eine wirklich große Matrix (mindestens 1000x1000) bestehend aus Nullen. 
  • Ersetzen Sie an zufälligen Stellen 1% der Nullen mit Zufallszahlen.
  • Erzeugen Sie aus dieser Matrix eine Sparse-Matrix mit einem anderen Namen.
  • Speichern Sie beide Matrizen jeweils einzeln in eine Datei.
  • Vergleichen Sie die Dateigrößen der beiden Dateien.
  • Erzeugen Sie eine gleichgroße Matrix, die komplett mit Zufallszahlen gefüllt ist. Und erzeugen Sie auch für diese eine Sparse Matrix-Variante.
  • Speichern Sie beide Zufalls-Matrizen einzeln in Dateien.
  • Was fällt Ihnen bei den vier Dateigrößen auf?

D) Hausaufgaben:

T6H1) Im Programm [vogelfang.m] werden drei verschiedene Arten von Zufallszahlen benutzt. Vollziehen Sie dieses Programm nach.
Nehmen Sie schrittweise folgende Änderungen vor:
a) Bei Amseln gibt es 60% Weibchen.
b) Bei Spatzen streut das Gewicht von Weibchen 3 Mal mehr als das Gewicht von Männchen.
c) Es kommen 25% Meisen und 25% Spatzen in der Gegend vor.

T6H2) In einem psychophysikalischen Experiment sollen einem Versuchstier drei verschiedene Töne in zufälliger Reihenfolge vorgespielt werden, aber jeder Ton soll genau 5 Mal vorkommen.

  • Wir kümmern uns erstmal nicht um die Generierung der Töne, sondern nennen sie einfach Bedingung 1, 2 und 3.
  • Überlegen Sie sich einen Algorithmus, der die Reizbedingungen in die richtige Reihenfolge bringt und setzen Sie diesen in ein Programm um.
  • Testen Sie das Programm, indem Du Sie es mehrfach laufen lassen. Tut es immer, was es soll? Sind die Ergebnisse jedes mal gleich?
  • Erweitern Sie Ihr Programm so, dass es N (eine beliebige Anzahl) Reize, die M mal (also beliebig oft) vorgespielt werden sollen, in eine Reihenfolge bringt.

*T6H3)  Die Messwerte einer Apparatur ist selbst ohne biologisches Präparat nicht perfekt rauschfrei. Um das Geräterauschen abzuschätzen, wurden im Elektrophysiologie-Praktikum für die Apparatur mit einer Modellzelle (einem elektronischen Schaltkreis, der die Membraneigenschaften einer Nervenzelle nachbaut) 100 Messungen mit dem gleichen Reiz [stimulus.mat] durchgeführt und die Antworten als Matrix unter [antworten1khz.mat] abgespeichert.

  • Schauen Sie sich ein beliebiges "trial" zusammen mit dem Reiz an (entsprechend Aufgabe C3).
  • Berechnen und plotten Sie den Zeitverlauf der über die 100 Messungen gemittelten Antwort.
  • Berechnen und plotten Sie den Mittelwert und die Standardabweichung der jeweils letzten 300ms für jede Messung (Mittelung über die Zeit). Gibt es eine Tendenz? Gibt es Ausreißer?
  • Sind Mittelwert und / oder Standardabweichung vor, während und nach der Reizung unterschiedlich?

*T6H4) In den Neurowissenschaften ist es eine verbreitete Herangehensweise, die Zeitpunkte des Auftretens einzelner Aktionspotentiale als unwichtig anzusehen und nur die Spikerate (also Anzahl aufgetretener Aktionspotentiale in einem gegebenen Zeitfenster) zu betrachten. Wenn das Auftreten einzelner Spikes simuliert werden soll, wird dafür häufig ein spezieller Zufallsprozess verwendet, der sogenannte Poissonprozess. Ein Poissonprozess hat die Eigenschaft, das zu jedem Zeitschritt unabhängig von allen anderen Zeitschritten entschieden wird, ob ein Ereignis stattfindet oder nicht. Wenn z.B. ein Neuron mit 10 Spikes/s feuert und eine Sampelrate von 1000Hz benutzt wird, ist die Wahrscheinlichkeit für das Auftreten eines Spikes in jedem Zeitschritt 1%.

a) Schreiben Sie eine Funktion, die als Eingabewerte eine Spikerate, die Samplerate  und die Länge der gewünschten simulierten Antwortspur bekommt und einen Vektor mit  der simulierten Antwort des Neurons zurückgibt.

b) Erweitern Sie Ihre Funktion so, dass sie nicht einen festen Wert für die Spikerate als Eingabe bekommt, sondern einen Vektor, der die Spikerate für jeden Zeitpunkt der zu simulierenden Antwort angibt (die Länge braucht man dann nicht mehr zu übergeben, sie ergibt sich aus der Länge des Vektors der Spikeraten).

*T6H5) Noch eine kleine Grafik-Übung: Nehmen wir an, Sie haben ein paar besonders empfindliche Fische in Ihrem Aquarium, in dem die Messreihe  [phWerte.mat] vorgenommen wurde. Diese Fische vertragen den Bereich von pH 6.5 bis pH 7.5 gut, sind aber darüber und darunter gefährdet. Definieren Sie außerdem in der Messreihe die Abweichungen um mehr als 2 Standardabweichungen als Messfehler. Plotten Sie die ph-Werte so, dass sie die Messwerte innerhalb und außerhalb Toleranzbereichs sowie die Messfehler mit Symbolen in drei unterschiedlichen Farben darstellen.

*T6H6) Die Funktion hist kann man auch auf Matrizen anwenden, dann  wird jede Spalte als Histogramm dargestellt. Erzeugen Sie sich eine Matrix mit zwei gegeneinander verschobenen Verteilungen und stellen diese grafisch dar.

Zum 7. Kurstag



Webmaster: Jutta Kretzberg (jutt19ffpa.kgoblqretzberg@uol.qujfde) (Stand: 21.08.2020)