Remove["Global`*"]

(* w1A *)

Nichtlineare Regression

Erster Datensatz, erzeugt mit  2 Exp[3 x] + 10 Exp[-x]  + “Fehler” (Random)

(* data3a *)

data3a={{0.05`,11.994333305587848`},{0.1`,11.65560900920308`},{0.15`,11.68107958893865`},{0.2`,12.180890713126544`},{0.25`,11.81781637246541`},{0.3`,11.958800259782324`},{0.35`,13.105524015844598`},{0.4`,13.039998993622376`},{0.45`,14.11706916628838`},{0.5`,14.70085639637875`},{0.55`,16.260041446614025`},{0.6`,17.628734779867973`},{0.65`,18.890158144959827`},{0.7`,21.171861988836092`},{0.75`,23.894269331951605`},{0.8`,26.748812570129424`},{0.85`,29.431339733648624`},{0.9`,33.54982103499556`},{0.95`,38.557965501995284`},{1.`,43.643313762836556`}};

(* Modell 1 *)

ymodell[x_ ]=  a Exp[b x] + c Exp[d x]

param=FindFit[data3a,ymodell[x] ,{a,b,c,d},x]

yfit1[x_] = ymodell[x]/.param

g1=ListPlot[data3a]

g2=Plot[yfit1[x],{x,0,1}]

Show[g1,g2]

(* Modell 2 *)

ymodell[x_ ]= a Exp[b x] + c

param=FindFit[data3a,ymodell[x],{a,b,c},x]

yfit2[x_] = ymodell[x] /.param

g2=Plot[yfit2[x],{x,0,1}]

Show[g1,g2]

(* Modell 3 *)

ymodell[x_ ]=a x^4 + b x^3 + c x^2 + d x + e

param=FindFit[data3a,ymodell[x ],{a,b,c,d,e},x]

yfit3[x_] = ymodell[x ]/.param

g2=Plot[yfit3[x],{x,0,1}]

Show[g1,g2]

(* Bewertung über die Maßzahlen: *)

(* die Anzahl der Messwerte-Paare *)

n = Length[data3a]

(* die x-Werte der Messung: *)

data3aX= Table[data3a[[i,1]],{i,1,n}]

(* die y-Werte der Messung: *)

data3aY= Table[data3a[[i,2]],{i,1,n}]

(* zunächst die jeweilige Summe der Schwankungsquadrate = Summe der quadrierten Residuen: *)

SQR1=Sum[(data3aY[[i]]- yfit1[data3aX[[i]]])^2,{i,1,n}]

SQR2=Sum[(data3aY[[i]]- yfit2[data3aX[[i]]])^2,{i,1,n}]

SQR3=Sum[(data3aY[[i]]- yfit3[data3aX[[i]]])^2,{i,1,n}]

(* jetzt die Varianzen *)

VarRes1=SQR1/(n-4)

VarRes2=SQR2/(n-3)

VarRes3=SQR3/(n-5)

Die Varianzen zeigen: Modell 1 ist am besten geeignet, aber die Varianzen der Fitfunktionen für Modell 1 und Modell 3 
unterscheiden sich nicht so sehr. Es sind also weitere Untersuchungen nötig: 
insbesondere müsste man ein weiteres Experiment vornehmen mit x > 1,  
siehe die beiden nächsten Abschnitte.

Vergleich der Fits für die Modelle 1 (rot) und 3 (grün) untereinander und mit dem Datensatz (blau):

g2=Plot[{yfit1[x],yfit3[x]},{x,0,0.2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]

g3=Plot[{yfit1[x],yfit3[x]},{x,0.2,1},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0],RGBColor[0,0,1]}]

Plot[{yfit1[x] - yfit3[x]},{x,0,1},PlotRange->All,PlotStyle->{Black}]

Plot[{(yfit1[x]-yfit3[x])/yfit1[x]},{x,0,1},PlotRange->All,PlotStyle->{Black}]

(* in der Grafik erkennt man: die Fits für die Modelle 1 und 3 unterscheiden sich offenbar nur am Anfang *)

data3a1= Take[data3a,4]

data3a2= Take[data3a,-17]

g1a= ListPlot[data3a1];

g1b= ListPlot[data3a2];

Show[g1a,g2]

Show[g1b,g3]

(* Vergleich der Residuen: *)

residuals1=data3aY-Map[yfit1,data3aX]  (* Map wirkt auf die angegebenen Listen *)

l1=ListPlot[residuals1,Filling->Axis,DataRange->{Min[data3aX],Max[data3aX]},PlotRange->{{0,1.1},All},PlotStyle-> Red]

residuals3=data3aY-Map[yfit3,data3aX]  (* Map wirkt auf die angegebenen Listen *)

l3=ListPlot[residuals3,Filling->Axis,DataRange->{Min[data3aX],Max[data3aX]},PlotRange->{{0,1.1},All}]

Show[l1,l3]

(* nochmal die jeweilige SQR, direkt über die Listen berechnet:  *)

residuals1SumQ=Total[(residuals1)^2]

residuals3SumQ=Total[(residuals3)^2]

TEST: Taylor-Entwicklung der Fit-Funktion (Modell 1), um dieses mit dem Polynom 4. Grades (Modell 3) zu   
            vergleichen

yfit3[x]

yfit1[x]

s1=Series[yfit1[x],{x,0,4}]//Normal

(* Taylor von Modell 1: rot ; Modell 3:  grün *)

Plot[{s1,yfit3[x]},{x,0,1},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[{s1,yfit3[x]},{x,0,0.2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[0,1,0]}]

Plot[{s1-yfit3[x]},{x,0,1},PlotRange->All]

Plot[{(s1-yfit3[x])/s1},{x,0,1},PlotRange->All]

Zweiter Datensatz, erzeugt mit 2 Exp[3 x] + 0.1 Exp[5x] + “Fehler” (Random)

(* data3b *)

data3b={{0.05`,2.456288136607986`},{0.1`,2.5805689320636755`},{0.15`,3.8340528262900753`},{0.2`,3.5778376792647215`},{0.25`,4.39222865000867`},{0.3`,5.754414298307113`},{0.35`,5.976053884314916`},{0.4`,7.099298224811543`},{0.45`,8.633036594324585`},{0.5`,9.994712460040049`},{0.55`,12.337349499257485`},{0.6`,13.93169332834842`},{0.65`,16.38016360496567`},{0.7`,19.430757440466984`},{0.75`,23.353434846562408`},{0.8`,27.276428947837427`},{0.85`,32.89705823180665`},{0.9`,38.77167890133632`},{0.95`,46.35542857694972`},{1.`,54.60128444103972`}};

(* Modell 1 *)

ymodell[x_ ]=  a Exp[b x] + c Exp[d x]

param=FindFit[data3b,ymodell[x] ,{a,b,c,d},x]

yfit1[x_] = ymodell[x]/.param

(* diese Fit-Funktion stimmt NICHT gut überein mit der [im Experiment nicht bekannten!] Original-Fkt; 
   die Daten werden aber offenbar ganz gut gefittet:  *)

g1=ListPlot[data3b]

g2=Plot[yfit1[x],{x,0,1}]

Show[g1,g2]

(* Modell 2 *)

ymodell[x_ ]= a Exp[b x] + c

param=FindFit[data3b,ymodell[x],{a,b,c},x]

yfit2[x_] = ymodell[x] /.param

g2=Plot[yfit2[x],{x,0,1}]

Show[g1,g2]

(* Modell 3 *)

ymodell[x_ ]=a x^4 + b x^3 + c x^2 + d x + e

param=FindFit[data3b,ymodell[x ],{a,b,c,d,e},x]

yfit3[x_] = ymodell[x ]/.param

g2=Plot[yfit3[x],{x,0,1}]

Show[g1,g2]

(* Bewertung über die Maßzahlen: *)

(* die Anzahl der Messwerte-Paare *)

n = Length[data3b]

(* die x-Werte der Messung: *)

data3bX= Table[data3b[[i,1]],{i,1,n}]

(* die y-Werte der Messung: *)

data3bY= Table[data3b[[i,2]],{i,1,n}]

(* zunächst die jeweilige Summe der Schwankungsquadrate = Summe der quadrierten Residuen: *)

SQR1=Sum[(data3bY[[i]]- yfit1[data3bX[[i]]])^2,{i,1,n}]

SQR2=Sum[(data3bY[[i]]- yfit2[data3bX[[i]]])^2,{i,1,n}]

SQR3=Sum[(data3bY[[i]]- yfit3[data3bX[[i]]])^2,{i,1,n}]

(* jetzt die Varianzen *)

VarRes1=SQR1/(n-4)

VarRes2=SQR2/(n-3)

VarRes3=SQR3/(n-5)

Man findet kaum einen Unterschied! Ein weiteres Experiment ist erforderlich: 
für x > 1 unterscheiden sich zumindest die exponentiellen Fitfunktionen der 
Modelle 1 & 2 von  der polynomialen Fitfunktion 3:

Plot[{yfit1[x],yfit2[x],yfit3[x]},{x,0,2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

(* die beiden exponentiellen Modellfunktionen unterscheiden sich kaum: *)

Plot[{yfit1[x]-yfit2[x]},{x,0,2},PlotRange->All,PlotStyle->{Black}]

Plot[{(yfit1[x]-yfit2[x])/yfit1[x]},{x,0,2},PlotRange->All,PlotStyle->{Black}]

TEST: Taylor-Entwicklung der Fit-Funktionen

yfit3[x]

s1=Series[yfit1[x],{x,0,4}]//Normal

s2=Series[yfit2[x],{x,0,4}]//Normal

(* die Taylor-Polynome der beiden exponentiellen Fitfunktionen unterscheiden sich kaum: *)

Plot[{s1,s2,yfit3[x]},{x,0,1},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

Plot[{s1,s2,yfit3[x]},{x,0,2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

Weitere mathematische Tests für Modell 1

(* 1. Startwerte für die Parametersuche aus der Originalfkt:  2 Exp[3 x] + 0.1 Exp[5 x], 
       die man im Experiment aber nicht kennt, und die FindFit auch nicht braucht; 
       beachte ausserdem die Hilfe: “In the nonlinear case, it finds in general only a locally optimal fit.”  *)

ymodell[x_ ]=  a Exp[b x] + c Exp[d x]

param=FindFit[data3b,ymodell[x] ,{{a,2},{b,3},{c,0.1},{d,5}},x]

ymodell[x ]/.param

(* also das alte Ergebnis: *)

param=FindFit[data3b,ymodell[x] ,{a,b,c,d},x]

yfit1[x_] = ymodell[x]/.param

Bemerkung:
 
 Additive hat eine Funktion zur Verfügung gestellt, welche zufällige Startpunkte im vierdimensionalen Parameterraum auswählt 
 und von dort aus FindFit startet. Damit soll verhindert werden, dass der Algorithmus in ein lokales Minimum läuft,  
 mit anderen  Worten: man hofft, so das globale Minimum zu finden. 
 Es kommt bei vielen Startpunkten zu Konvergenzproblemen; der Fit mit der minimalen SQR ist derjenige mit obigen Parameterwerten.

(* 2. Berechnung der Parameter aus der Mininimalbedingung für die Summe der quadrierten Abweichungen, SQA;
       FindMinimum braucht keinen Startwert; beachte die Hilfe: 
       “searches for a local minimum in f, starting from an automatically selected point.” *)

SQA[a_,b_,c_,d_]=Sum[(data3bY[[i]]- ymodell[data3bX[[i]]])^2,{i,1,n}]

FindMinimum[SQA[a,b,c,d],{a,b,c,d}]

(* also das alte Ergebnis *)

(* 3. Parameterberechnung "von Hand": Nullstellen des Gradienten von SQA *)

gla= D[SQA[a,b,c,d],a]==0

glb= D[SQA[a,b,c,d],b]==0

glc= D[SQA[a,b,c,d],c]==0

gld= D[SQA[a,b,c,d],d]==0

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}}]

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},MaxIterations->1000]

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},MaxIterations->10000]

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},MaxIterations->1000,WorkingPrecision->100]

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},MaxIterations->1000,WorkingPrecision->20]

(* Da in SQA die MESSWERTE eingehen, kann man die Genauigkeit der "argument function" nicht erhöhen *)

(* laut Hilfe: Relax error tolerances for stopping: *)

FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},AccuracyGoal->4,PrecisionGoal->4]

(* laut Hilfe: Make estimated relative distance to the root the main criterion for stopping: *)

param=FindRoot[{gla,glb,glc,gld},{{a,2},{b,3},{c,0.1},{d,5}},AccuracyGoal->Infinity,PrecisionGoal->4]

yFindRoot[x_]=ymodell[x]/. param

Bemerkungen:

1. Da das Verfahren nicht konvergiert, sind die Ergebnisse OHNE BELANG!    

2. FindRoot benutzt zur Nullstellen-Bestimmung andere Algorithmen als FindMinimum zur Minimum-Bestimmung.
    Mit seinen Algorithmen kann FindRoot aber kein verlässliches Ergebnis erzeugen.

3. Obwohl die von FindRoot ermittelten Parameter-Werte den Werten der Originalfkt 2 Exp[3 x] + 0.1 Exp[5 x] sehr nahe kommen, 
    stellt diese Fitfunktion den Datensatz schlechter dar als die drei anderen (obigen) Fitfunktionen:

SQR=Sum[(data3bY[[i]]- yFindRoot[data3bX[[i]]])^2,{i,1,n}]

SQR/(n-4)

Auch die Original-Funktion, die man natürlich im Experiment NICHT kennt, stellt den Datensatz NOCH schlechter dar:

y[x_] = 2 Exp[3 x] + 0.1 Exp[5x]

SQRO=Sum[(data3bY[[i]]- y[data3bX[[i]]])^2,{i,1,n}]

Die Zufallszahlen, welche bei der Erzeugung des Datensatzes zu den Funktionswerten der  
Originalfunktion addiert werden, verändern diese Werte derart, 
dass die Fitfunktion  yfit1  den Datensatz besser darstellt als die Originalfunktion. 
Mit anderen Worten:  yfit1 "erklärt" die Daten besser als die Originalfunktion. 
 
ABER:
 
Mal angenommen, 2 Exp[3 x] + 0.1 Exp[5 x] wäre der "wahre" funktionale Zusammenhang, 
sozusagen das Naturgesetz, welches man finden will. Dann liefert dieses Experiment 
zwar eine Fitfunktion, welche die Daten gut erklärt, aber nicht dieses "Naturgesetz" wiedergibt.
Die gute Erklärung des experimentellen Befundes (des Datensatzes) durch diese Fitfunktion wäre 
nur eine auf diesen Fall passende phänomenologische Beschreibung.

Eine Theorie (hier: ein “Naturgesetz”) muss sich aber auch dahingehend bewähren, dass man Voraussagen machen kann, die wiederum einer experimentellen Überprüfungen standhalten (ansonsten wird die Theorie verworfen bzw. modifiziert). 
Wenn man mit der Fitfunktion eine Vorhersage macht (hier: für x>1), ergeben sich deutliche Abweichungen zum “Naturgesetz”:

Plot[{yfit1[x],y[x]},{x,0,2},PlotStyle->{RGBColor[1,0,0],RGBColor[0,0,1]},PlotRange->All](* Fitfunktion rot, Originalfkt blau *)

Plot[{yfit1[x]-y[x]},{x,0,2},PlotRange->All]

Ein weiteres Experiment ist erforderlich.

Dritter Datensatz: der zweite erweitert auf   0 < x < 2

(* data3c *)

(* ein neues Experiment *)

data3c={{0.05`,2.1239258453435177`},{0.1`,2.846062774283036`},{0.15`,3.436801895759832`},{0.2`,3.9425000732955713`},{0.25`,4.402529457974034`},{0.3`,5.117726798631371`},{0.35`,5.819753708982341`},{0.4`,7.725754266628263`},{0.45`,8.715124719358945`},{0.5`,10.66790078974186`},{0.55`,12.419074299076549`},{0.6`,14.455258407767813`},{0.65`,16.89758647232464`},{0.7`,19.416240955899347`},{0.75`,23.58468445646571`},{0.8`,27.60851225833905`},{0.85`,32.53003459349874`},{0.9`,38.482300370764776`},{0.95`,46.1684567488445`},{1.`,54.65502229206134`},{1.05`,65.79573219626845`},{1.1`,78.97053659553973`},{1.15`,94.8138206329427`},{1.2`,113.72247959408686`},{1.25`,136.35200353210237`},{1.3`,165.00701955902207`},{1.35`,200.41753816293993`},{1.4`,243.02362116218342`},{1.45`,295.88244719193403`},{1.5`,361.0182839614049`},{1.55`,441.50021609607256`},{1.6`,540.8729451985249`},{1.65`,665.0981054712286`},{1.7`,819.5501610496445`},{1.75`,1012.5810489946456`},{1.8`,1253.1262556041856`},{1.85`,1554.468048468654`},{1.9`,1934.1586774036734`},{1.95`,2409.9405197363108`},{2.`,3009.5557451683594`}};

(* Modell 1 *)

ymodell[x_ ]=  a Exp[b x] + c Exp[d x]

param=FindFit[data3c,ymodell[x] ,{a,b,c,d},x]

yfit1[x_] = ymodell[x]/.param

g1=ListPlot[data3c]

g2=Plot[yfit1[x],{x,0,2}]

Show[g1,g2]

(* Modell 2 *)

ymodell[x_ ]= a Exp[b x] + c

param=FindFit[data3c,ymodell[x],{a,b,c},x]

Obige Warnung beachten !! Das Ergebnis ist nicht vertrauenswürdig! 
Beobachtung:dieses Modell war im vorausgegangenen "Experiment" das BESTE! 
Die Parameterwerte waren aber sehr anders: 0.5184967857190231` +1.673954905751382` (e^(3.478079069037661` x)) 
Beachten Sie auch die völlig veränderten Bewertungszahlen (s.u.)!

yfit2[x_] = ymodell[x] /.param

g2=Plot[yfit2[x],{x,0,2}]

Show[g1,g2]

(* Modell 3 *)

ymodell[x_ ]=a x^4 + b x^3 + c x^2 + d x + e

param=FindFit[data3c,ymodell[x ],{a,b,c,d,e},x]

yfit3[x_] = ymodell[x ]/.param

g2=Plot[yfit3[x],{x,0,2},PlotRange->All]

Show[g1,g2]

(* Bewertung über die Maßzahlen: *)

(* die Anzahl der Messwerte-Paare *)

n = Length[data3c]

(* die x-Werte der Messung: *)

data3cX= Table[data3c[[i,1]],{i,1,n}]

(* die y-Werte der Messung: *)

data3cY= Table[data3c[[i,2]],{i,1,n}]

(* zunächst die jeweilige Summe der Schwankungsquadrate = Summe der quadrierten Residuen: *)

SQR1=Sum[(data3cY[[i]]- yfit1[data3cX[[i]]])^2,{i,1,n}]

SQR2=Sum[(data3cY[[i]]- yfit2[data3cX[[i]]])^2,{i,1,n}]

SQR3=Sum[(data3cY[[i]]- yfit3[data3cX[[i]]])^2,{i,1,n}]

(* jetzt die Varianzen *)

VarRes1=SQR1/(n-4)

VarRes2=SQR2/(n-3)

VarRes3=SQR3/(n-5)

Hier ist alles wieder, wie man es sich wünscht: das erste Modell ist zweifelsfrei das beste; 
graphisch lassen sich die Unterschiede aber nicht so gut ausmachen; 
dies liegt an der Auflösung: die Funktionswerte werden ja doch recht groß, 
die punktuelle Differenz dagegen eher klein.

g2=Plot[{yfit1[x],yfit2[x],yfit3[x]},{x,0,2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

Show[g1,g2,PlotRange->All]

Plot[{yfit1[x],yfit2[x],yfit3[x]},{x,1,2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

yfit1[x]

yfit2[x]

yfit3[x]

Plot[{yfit1[x],yfit2[x]},{x,0,2},PlotRange->All,PlotStyle->{RGBColor[1,0,0],RGBColor[1,1,0],RGBColor[0,0,1]}]

Plot[{yfit1[x]-yfit2[x]},{x,0,2},PlotRange->All,PlotStyle->{Black}]

Plot[{yfit1[x]-yfit3[x]},{x,0,2},PlotRange->All,PlotStyle->{Black}]

Abschließende Bemerkungen:

1. Bei einem echten (also in einem realen Experiment gemessenen) Datensatz würden diese 
    Unterschiede (data3b/data3c) möglicherweise nicht auftreten.

2. Aber auch dann gilt: traue keiner Statistik, auch nicht der eigenen, ohne gründliche Prüfung!

3. Traue keiner numerischen Berechnung ohne gründliche Prüfung!