Verfasst am: 06.07.2015, 13:24
Titel: Fitten von Modell an (Mess-)Daten ohne Toolbox-Abhängigkeit
Hallo zusammen,
es kommt immer häufiger vor, dass Leute gegebene (Mess-)Daten an ein vorgelegtes Modell fitten wollen und dafür die optimalen Parameter finden möchten. Oftmals führt als Lösung der Weg über die statistics- oder die curve-fitting-Toolbox.
Was aber, wenn man für die angegebenen Toolboxen keine Lizenzen verfügt, beispielsweise jenseits der Universität? Es gibt für genau so einen Fall die Funktion
fminsearch
. Da ich diesen Weg schon einige Male gegangen bin, habe ich einen Standardfall vorbereitet, der aus gegebenen Daten und einer Modellfunktion mit Angabe einer ersten Schätzung für die Parameter optimale Parameter zurückliefert.
Der damalige Anwendungsfall war das Finden der optimalen Lage der Kontrollpunkte für Bézierkurven. Die Variable
order
gab dafür die Anzahl der Kontrollpunkte an. In der Variable
p
standen die x- und y-Koordinaten dieser Kontrollpunkte. Diese Details waren jedoch alle in der
fit_function
(s.u.) platziert, sodass das anfängliche Rahmenwerk sehr allgemein ist. Aber ganz so kompliziert möchte ich es hier nicht machen.
In dem angegebenen Beispiel erzeuge ich zufällige Daten für x und für y verrauschte Daten eines Polynoms. Anschließend möchte ich die optimalen Parameter (Koeffizienten des Polynoms) finden. Natürlich macht man das eigentlich mit
polyfit
. Das mache ich aber absichtlich nicht, da ich gerne
fminsearch
benutzen möchte.
In der
fit_function
habe ich noch andere Fälle als Kommentar angefügt, die verdeutlichen sollen, dass grundsätzlich alles denkbare damit gelöst werden kann. Ich kann damit durch die Messdaten auch eine Summe eines Sinus und einer Normalparabel legen - sofern ich denn möchte. Es müssen nur die Parameter aus
p
entsprechend benutzt werden.
Viel Spaß und viel Erfolg mit dem Beispiel.
Code:
function main
% Beispiel für eine Optimierung der Parameter für einen Fit an Messdaten mit fmisearch - ohne Abhängigkeit zu irgendeiner Toolbox.
%% ----- BEISPIELDATEN ERZEUGEN----- % % Wir erzeugen als Beispieldaten einige zufällige x-Werte und als % zugehörige y-Werte ein Polynom verrauscht mit Zufallszahlen. % In einem echten Fall hat man für gewöhnlich bereits Messdaten für x und y % vorliegen. Beispielsweise könnte x die Zeit beschreiben und y eine % gemessene Größe (Temperatur, Druck, Windgeschwindigkeit, ...).
x = sort(rand(30, 1))*4;
p_original = [1.5, -1.5, -5, 1];
noise = rand(size(x)) - 0.5;
y = polyval(p_original,x) + noise;
%% ----- WEITERE OPTIMIERUNGSPARAMETER ----- % % Für bestimmte Fälle, insbesondere bei komplizierten fit-Funktionen will % man weitere Optimierungsparameter übergeben. Als Beispiel übergeben wir % noch den Grad des Polynoms, das wir an die Daten fitten wollen.
order = 3;
%% ----- FINDEN DER PARAMETER ----- % % In diesem Abschnitt werden die optimalen Parameter gesucht. Wenn wir % wirklich ein Polynom vom grad order fitten wollen, würde man das % selbstverständlich mit polyfit() erledigen: % % p_fit = polyfit(x, y, order); % % Wir wollen jedoch in diesem Beispiel mit fminsearch arbeiten. So haben % wir die Möglichkeit, ohne die Abhängigkeit von speziellen Toolboxen % beliebig komplexe fit-Funktionen zu fitten. Dafür rufen wir die eigens % geschriebene Funktion optimierung() mit unseren Daten x und y, sowie dem % zusätzlichen Argument "order" auf.
p_fit = optimierung(x, y, order);
%% ----- OPTIONAL: ERGEBNISSE DARSTELLEN ----- % % Anschließend bietet es sich an, das Ergebnis zu visualisieren. figure
x_fit = linspace(min(x), max(x), 100); % Mehr Stützstellen für weichen plot plot(x, y, 'ko')% Original gestörte Daten hold on
plot(x_fit, polyval(p_original, x_fit), 'k-')% Original Polynom plot(x_fit, fit_function(x_fit, p_fit), 'g')% Ergebnis von fminsearch legend('location', 'northwest', 'Gestörte Datenpunkte', ...
'Originales Polynom', 'Lösung mit Optimierer') end
%% Hilfsfunktionen function p_fit = optimierung(x, y, order) % In dieser Methode wird die Optimierung durchgeführt, indem fminsearch % mit der zu minimierenden Funktion objFun und zusätzlich einem Startwert % p0 aufgerufen wird. % Wir legen der Einfachheit halbe einer Polynom oder Ordnung order durch x % und y. Dafür brauchen wir order+1 Koeffizienten.
p0 = ones(1, order+1); % wir raten: alle Koeffizienten sind 1
p_fit = fminsearch(@objFun, p0);
% ----- BEGIN: NESTED objFUN -----% function e = objFun(p) % Die Funktion objFun bekommt nur p als Parameter übergeben und liefert den % Wert e zurück, welcher dabei minimiert wird. e ist dabei das Maß, das % beschreibt, wie gut die aktuell gefundenen Parameter sind, um mit dem % Modell (fit_function) die Daten (x und y) zu beschreiben. Dafür % vergleicht man das Ergebnis der fit_function mit den Meswerten. Meistens % ist dafür "least squares" also die Summe der Fehlerquadrate angebracht. % In dieser "nested function" leben auch x und y (und ebfennalls auch p0, % und order), sodass dieser Vergleich gemacht werden kann.
yfitted = fit_function(x, p);
e = sum((yfitted-y).^2);
end % ----- END: NESTED objFUN -----%
function yfitted = fit_function(x, p) % In der fit_function steht die gesamte Magie beschrieben, wie aus den % x-Werten und dem Vektor "p" zu jedem x-Wert ein y-Wert aus unserem Modell % gemacht wird. Hierbei können beliebig andere Funktionen aufgerufen % werden, je nach Komplexität des Problems. % % Wir wollen in diesem Beispiel ein Polynom fitten (Fall 1). Dabei ist p % das Polynom in der Matlab-Repräsentation als Vektor der Koeffizienten, das % an den Stellen x ausgewertet wird. % % Zusätzlich sind noch andere Beispiele angeführt (Fall 2 bis 4), deren % Sinnhaftigkeit aber lieber nicht näher untersucht werden sollte, sie % dienen bloß zur Veranschaulichung des Prinzips für komplexere % fit_functions.
% Fall 2-4: Hier kann auch irgendetwas seltsames stehen. % Das versteckt man alles in dem p, das übergeben wurde. Die fit-function % muss bloß die Information in dem p richtig deuten.
% Fall 2: Summe aus Sinus-Welle und skalierter Standardparabel % A = p(1); % w = p(2); % f = p(3); % sin_wave = A*sin(2*pi*f*x + w); % parabola = p(4)*x.^2; % yfitted = sin_wave + parabola;
% Fall 3: Summe aus einer Exponentialfunktion und einer beliebigen Gerade % exp_fun = p(1)*exp(p(2)*x); % straight_line = polyval(p(3:4), x); % yfitted = exp_fun + straight_line;
% Fall 4: Produkt aus Wurzel und Exponentialfunktion + Verschiebung % yfitted = p(1)*sqrt(x).*exp(p(2)*x + p(3)) + p(4);
hallo nras
sehr schönes tutorial. vielen dank dafür.
ich lass gleich nochmal den link zum kurzen beispiel von matlab da
Optimal Fit of a Nonlinear Function dort wird auch mit fminsearch gearbeitet.
grüße
_________________
ich sehe es so, dass
fminsearch
dann die Methode der Wahl ist, wenn weder Statistics Toolbox noch Optimization Toolbox zugänglich sind.
Falls eine der Toolboxen zur Verfügung steht, sind
lsqcurvefit
(Optimization Toolbox) oder
fitnlm
(Statistics Toolbox; in älteren Versionen NonLinearModel.fit) einfacher in der Anwendung.
Da die Algorithmen hinter
lsqcurvefit
und
fitnlm
gradientenbasiert sind, konvergieren sie zudem meist schneller als
fminsearch
(Direktsuche).
Die von fitnlm erzeugten NonLinearModel-Objekte machen viele Informationen, z.B. R², leicht zugänglich und haben ein paar Methoden, z.B. coefCI zur Berechnung von Konfidenzintervallen, die man darauf anwenden kann:
Falls die Modellfunktion linear in den Parametern ist, sind
lsqlin
bzw.
fitlm
(früher: LinearModel.fit) am besten geeignet und erfordern keine Startwerte.
Die Curve Fitting Toolbox beinhaltet schließlich noch ein interaktives Tool, das keinerlei Programmierkenntnisse erfordert.
Wie eingangs erwähnt ist dieses Tutorial für genau diesen Fall gedacht.
Das war bei mir angekommen. Ich denke auch, dass das Tutorial diesen Zweck erfüllt.
Es ist aber ja gut möglich, dass Leute auf diesen Thread stoßen, die die benötigten Toolboxen haben. Daher hielt ich diesen Thread für einen guten Platz, auf die Toolbox-Funktionalitäten und ihre Unterschiede zu fminsearch als auch untereinander hinzuweisen.
Ich habe versucht deinen Code auf meine Problemstellung zu übertragen. Leider bekomme ich keine Lösung.
Der Unterschied zu deinem Beispiel sind zum einen die Messdaten (Zeilen 10 / 11) und zum anderen die Zielfunktion (Zeile 88 ).
Falls Du oder jemand anders hier im Forum eine Idee hat, wo der Fehler liegt, bitte melden
Code:
function main
% Beispiel für eine Optimierung der Parameter für einen Fit an Messdaten mit fmisearch - ohne Abhängigkeit zu irgendeiner Toolbox.
%% ----- BEISPIELDATEN ERZEUGEN----- % % Wir erzeugen als Beispieldaten einige zufällige x-Werte und als % zugehörige y-Werte ein Polynom verrauscht mit Zufallszahlen. % In einem echten Fall hat man für gewöhnlich bereits Messdaten für x und y % vorliegen. Beispielsweise könnte x die Zeit beschreiben und y eine % gemessene Größe (Temperatur, Druck, Windgeschwindigkeit, ...).
x = [00.511.3]*1e-3;
y = [00.380.550.62]*1e-4;
%% ----- WEITERE OPTIMIERUNGSPARAMETER ----- % % Für bestimmte Fälle, insbesondere bei komplizierten fit-Funktionen will % man weitere Optimierungsparameter übergeben. Als Beispiel übergeben wir % noch den Grad des Polynoms, das wir an die Daten fitten wollen.
order = 3;
%% ----- FINDEN DER PARAMETER ----- % % In diesem Abschnitt werden die optimalen Parameter gesucht. Wenn wir % wirklich ein Polynom vom grad order fitten wollen, würde man das % selbstverständlich mit polyfit() erledigen: % % p_fit = polyfit(x, y, order); % % Wir wollen jedoch in diesem Beispiel mit fminsearch arbeiten. So haben % wir die Möglichkeit, ohne die Abhängigkeit von speziellen Toolboxen % beliebig komplexe fit-Funktionen zu fitten. Dafür rufen wir die eigens % geschriebene Funktion optimierung() mit unseren Daten x und y, sowie dem % zusätzlichen Argument "order" auf.
p_fit = optimierung(x, y, order);
%% ----- OPTIONAL: ERGEBNISSE DARSTELLEN ----- % % Anschließend bietet es sich an, das Ergebnis zu visualisieren. figure
x_fit = linspace(min(x), max(x), 100); % Mehr Stützstellen für weichen plot plot(x, y, 'ko')% Original gestörte Daten hold on
plot(x_fit, fit_function(x_fit, p_fit), 'g')% Ergebnis von fminsearch legend('location', 'northwest', 'Gestörte Datenpunkte', ...
'Lösung mit Optimierer') end
%% Hilfsfunktionen function p_fit = optimierung(x, y, order) % In dieser Methode wird die Optimierung durchgeführt, indem fminsearch % mit der zu minimierenden Funktion objFun und zusätzlich einem Startwert % p0 aufgerufen wird. % Wir legen der Einfachheit halbe einer Polynom oder Ordnung order durch x % und y. Dafür brauchen wir order+1 Koeffizienten.
p0 = ones(1, order+1); % wir raten: alle Koeffizienten sind 1
p_fit = fminsearch(@objFun, p0);
% ----- BEGIN: NESTED objFUN -----% function e = objFun(p) % Die Funktion objFun bekommt nur p als Parameter übergeben und liefert den % Wert e zurück, welcher dabei minimiert wird. e ist dabei das Maß, das % beschreibt, wie gut die aktuell gefundenen Parameter sind, um mit dem % Modell (fit_function) die Daten (x und y) zu beschreiben. Dafür % vergleicht man das Ergebnis der fit_function mit den Meswerten. Meistens % ist dafür "least squares" also die Summe der Fehlerquadrate angebracht. % In dieser "nested function" leben auch x und y (und ebfennalls auch p0, % und order), sodass dieser Vergleich gemacht werden kann.
yfitted = fit_function(x, p);
e = sum((yfitted-y).^2);
end % ----- END: NESTED objFUN -----%
function yfitted = fit_function(x, p) % In der fit_function steht die gesamte Magie beschrieben, wie aus den % x-Werten und dem Vektor "p" zu jedem x-Wert ein y-Wert aus unserem Modell % gemacht wird. Hierbei können beliebig andere Funktionen aufgerufen % werden, je nach Komplexität des Problems. % % Wir wollen in diesem Beispiel ein Polynom fitten (Fall 1). Dabei ist p % das Polynom in der Matlab-Repräsentation als Vektor der Koeffizienten, das % an den Stellen x ausgewertet wird. % % Zusätzlich sind noch andere Beispiele angeführt (Fall 2 bis 4), deren % Sinnhaftigkeit aber lieber nicht näher untersucht werden sollte, sie % dienen bloß zur Veranschaulichung des Prinzips für komplexere % fit_functions.
% Fall 1: Zielfunktion
yfitted = (p(1).*x+p(2).*x.^1.5)/sqrt(p(3)+p(4).*x.^2);
% Fall 2-4: Hier kann auch irgendetwas seltsames stehen. % Das versteckt man alles in dem p, das übergeben wurde. Die fit-function % muss bloß die Information in dem p richtig deuten.
% Fall 2: Summe aus Sinus-Welle und skalierter Standardparabel % A = p(1); % w = p(2); % f = p(3); % sin_wave = A*sin(2*pi*f*x + w); % parabola = p(4)*x.^2; % yfitted = sin_wave + parabola;
% Fall 3: Summe aus einer Exponentialfunktion und einer beliebigen Gerade % exp_fun = p(1)*exp(p(2)*x); % straight_line = polyval(p(3:4), x); % yfitted = exp_fun + straight_line;
% Fall 4: Produkt aus Wurzel und Exponentialfunktion + Verschiebung % yfitted = p(1)*sqrt(x).*exp(p(2)*x + p(3)) + p(4);
das Problem liegt in der Funktionsauswertung zum Plotten: du schickst einen x-Vektor rein, aber deine Funktion gibt da nur einen Skalar zurück. Wenn du / in der Zielfunktion durch ./ ersetzt, sieht es besser aus.
Davon abgesehen ist es nicht möglich, 4 Parameter sinnvoll durch 3 Datenpunkte zu bestimmen.
Grüße,
Harald
Einstellungen und Berechtigungen
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS