Erste Datenanalysen

Fitten: Funktionen an Daten anpassen

Im letzten Abschnitt haben wir die Daten aus der Datei "scripts/examples/amplitude.dat" in die Tabelle data() geladen (falls das bei Dir nicht so sein sollte, wiederholst Du am besten den letzten Abschnitt). Wir haben außerdem gesehen, dass die Daten Amplituden vs. Frequenz-Kurven enthalten, die jeweils ein Maximum besitzen. Es liegt also nahe, dass es sich hier um Resonanz-Messungen handelt. Als erste Analyse kannst Du also versuchen, eine Resonanzkurve (auch Lorentz-Funktion genannt) an die Daten anzupassen, um die Parameter der Messung zu extrahieren.

Dazu definieren wir uns zunächst die Lorentz-Funktion. Das ist zwar nicht zwingend nötig, macht den Code im Folgenden aber insgesamt lesbarer:

<- define lorentz(x, x0, amplitude, damping, offset) := amplitude/sqrt((x^2-x0^2)^2+4*damping^2*x^2) + offset

-> Die Funktionsdefinitionen wurden erfolgreich aktualisiert.

Die Dir derzeit noch unbekannten Parameter der Lorentz-Funktion x0, amplitude, damping und offset kannst Du Dir nun mittels des Kommandos fit aus den Daten extrahieren:

<- fit data(:,1:2) -with=lorentz(x, x0, amp, damp, off)

-> Fitte "amp/sqrt((x^2-x0^2)^2+4*damp^2*x^2) + off" ... Erfolg.

===========================================================================================

-> NUMERE: FITERGEBNIS

===========================================================================================

-> Funktion: 0.43003/sqrt((x^2-0.00083679^2)^2+4*3.0264^2*x^2) + 0.0080817

[...]

-> Parameter Initialwert Anpassung Asymptotischer Standardfehler

----------------------------------------------------------------------------------------

amp 0 0.4300348 ± 369.36 (8.589e+004%)

damp 0 3.026374 ± 3.0706e+006 (1.015e+008%)

off 0 0.008081738 ± 4.5169 (5.589e+004%)

x0 0 0.0008367877 ± 2.2077e+010 (2.638e+015%)

[...]

===========================================================================================

Bei Dir wird im Kommandofenster hier noch viel mehr Text stehen. Wir haben die sehr lange Ausgabe, die fit erzeugt, an dieser Stelle etwas abgekürzt und wollen uns erst einmal auf das Ergebnis fokussieren. In der Tabelle findest Du die alphabetisch sortierte Liste der Parameter mit deren Initial- und angepassten Werten sowie einer Schätzung für deren Genauigkeit. Wenn Du die prozentuale Abweichung ansiehst, dann wirst Du sicher vermuten, dass hier etwas nicht stimmt. Falls Du es eher gewohnt bist, mit dem Bestimmtheitsmaß R² zu arbeiten, kannst Du Dir dieses mit der folgenden Zeile auch berechnen lassen:

<- 1 - chi^2 / sum((data(:, 2) - avg(data(:, 2)))^2)

-> ans = 0.08690886

Auch hier zeigt sich, dass offenbar kein besonders guter Match gegeben ist. Die Variable chi wurde dabei automatisch von fit erzeugt. Zur verbesserten Analyse stellen wir einfach mal die Daten mit dem Fit zusammen graphisch dar.

Leider ist das ein fundamentales Problem von Fitting-Algorithmen. Die Ergebnisse sind von den Startwerten abhängig und können durch leichte Änderungen dieser Parameter vollkommen anders ausfallen. Daher ist es sinnvoll, bei nicht-polynomialen Funktionen die ungefähren Parameter als Schätzungen vorzugeben.

Bis jetzt haben wir für alle Parameter beim Fitten den (automatischen) Startwert 0 verwendet (Korrekterweise ist das der automatische Initialwert der neu erzeugten Variablen). Das kannst Du aber auch explizit anders vorgeben, indem Du die params=PARAMS Option verwendest:

<- fit data(:,1:2) -with=lorentz(x, x0, amp, damp, off) params=[x0=2.5,amp=2,damp=0,off=0]

-> Fitte "amp/sqrt((x^2-x0^2)^2+4*damp^2*x^2) + off" ... Erfolg.

===========================================================================================

-> NUMERE: FITERGEBNIS

===========================================================================================

-> Funktion: 0.031963/sqrt((x^2-2.4575^2)^2+4*0.11405^2*x^2) + 0.007953

[...]

-> Parameter Initialwert Anpassung Asymptotischer Standardfehler

----------------------------------------------------------------------------------------

amp 2 0.0319629 ± 0.0087816 (27.47%)

damp 0 0.114046 ± 0.030319 (26.58%)

off 0 0.00795303 ± 0.0045607 (57.35%)

x0 2.5 2.457466 ± 0.015739 (0.6404%)

[...]

===========================================================================================

Wenn Du wieder das Bestimmtheitsmaß R² berechnest, so erhältst Du das folgende:

<- 1 - chi^2 / sum((data(:, 2) - avg(data(:, 2)))^2)

-> ans = 0.9443302

Fouriertransformationen mittels FFT

Die Frequenz- bzw. Fourieranalyse gehört mit zu dem wichtigsten Handwerkszeug in der Datenanalyse. Daher wollen wir Dir hier auch ein Beispiel zeigen, wie Du eine Fouriertransformation in NumeRe rechnen kannst. Dazu verwenden wir im folgenden das Kommando fft. Zuvor brauchen wir aber noch einen passenden Datensatz, den wir auch transformieren können. Gebe dazu die folgenden vier Zeilen in das Kommandofenster ein:

<- new signal();

<- signal(:, 1) = {0:0.001:0.5};

<- signal(:, 2) = sin(_2pi*50*signal(:, 1)) + sin(_2pi*120*signal(:, 1)) + gauss(0, 2);

<- signal(#, :) = "Time t [s]", "Signal x(t)";

Mit den ersten drei Zeilen haben wir einen etwas verrauschten Datensatz erzeugt, wie er Dir bei der Arbeit mit Audioaufnahmen begegnen kann. In der zweiten Zeile siehst Du außerdem eine Möglichkeit, einen Vektor mit iterativ befüllten Werten zu generieren: {A : B : C} = {A, A+B, A+2*B, ..., C}, wobei C nur dann Teil des Vektors wird, wenn es durch A+n*B ausdrückbar ist. Dabei ist außerdem {A : B} = {A : 1 : B}, bzw. {A : -1 : B}, falls B < A.

Für das Rauschen ist die Funktion gauss(x0,x1) verantwortlich, die Gauss-verteiltes Rauschen um die Position x0 mit der Halbwertsbreite x1 erzeugt. Wir haben außerdem die automatische Vektorisierung durch die Verwendung von Tabellenspalten (signal(:, 1)) genutzt. Die letzte Zeile ergänzt dann, wie Dir bereits bekannt ist, die Spaltenüberschriften.