Hallo Udo,
tut mir leid, dass ich so spät antworte, aber ich bin in der Zwischenzeit
umgezogen!
Das iterative Fitten mit den Gewichten auf Grundlage des Funktionswertes der
Modellfunktion an der Stelle x des Datenpunkts ist nun wirklich des Rätsels
Lösung! Ich habe es geprüft mit der Konstanten und dem Lorentz. Keine
Systematik mehr und man hat die Information über die Statistik der
Messpunkte mit einbezogen!
Vielen Dank für die Diskussion, ich fand sie sehr fruchtbar.
Viele Grüße,
Martin Haas
> Hallo nocheinmal,
>
> mir ist noch ein Fehler in der funktionalen Programmierung aufgefallen:
<snip>
> Das sollte nicht sein. Deshalb mit der Bitte um Entschuldigung das Ganze
> nocheinmal und mit recht vernünftigem Abbruchprozedere:
>
> In[8]:=
> Clear[regressionTest];
> regressionTest[l0_List, l1_List, l2_List] :=
> (Last[l2] != Last[l1]) && If[Last[l2] == Last[l0], If[Last[l2] <
> Last[l1],
> Print["Period catched!\nPenultimate set: a = ", l1[[2, 1]], " c =
> ", l1[[2, 2]], " g = ", l1[[2, 3]],
> " wSOfSq = ", Last[l1]]; False, True], True]
>
> In[66]:=
> Clear[regressionStep];
> regressionStep[{data_List, {a0_, c0_, g0_}, x0_ (* unused *)}] :=
> Module[{x1, a2, c2, g2, gew, qdev},
> gew = {}; qdev = {};
> gew = ((Abs[#[[2]] - Lorentz[#[[1]], a0, c0, g0]]/Lorentz[#[[1]], a0,
> c0, g0]) & ) /@ data;
> {a2, c2, g2} = {a1, c1, g1} /. (BestFitParameters /.
> NonlinearRegress[data, Lorentz[x1, a1, c1, g1], {x1}, {a1, c1, g1},
> Weights -> gew, RegressionReport -> BestFitParameters]);
> qdev = (((#[[2]] - Lorentz[#[[1]], a2, c2, g2])^2) & ) /@ data;
> (*
> Print[" {a2, g2, c2}: ", {a2, c2, g2}, "wSOfSq: ", gew . qdev];
> *)
> {data, {a2, c2, g2}, gew . qdev}
> ]
>
> In[80]:=
> (* using NestWhile *)
> With[{a = 20., c = 0., g = 1., xL = -2., xR = 2., anZ = 101},
> data = ({#, Lorentz[#, a, c, g] + Random[NormalDistribution[0.,
> sigma[Lorentz[#, a, c, g]]]]} & ) /@ Sort[Table[Random[Real, {xL, xR},
> $MachinePrecision], {anZ}]];
> If[Length[Select[Last[Transpose[data]], Negative]] > 0,
> Print["Negative data point(s) generated."];
> Return[$Failed]
> ];
> pic = ListPlot[data, DisplayFunction -> Identity];
> (* StartWert *)
> {a2, c2, g2} = {a1, c1, g1} /. (BestFitParameters /.
> NonlinearRegress[data, Lorentz[x1, a1, c1, g1], {x1}, {a1, c1, g1},
> Weights -> Automatic, RegressionReport -> BestFitParameters]);
> unwSOfSqStart = Plus @@ ((((#[[2]] - Lorentz[#[[1]], a2, c2, g2])^2) &
> ) /@ data);
> Print["StartParameters: a = ", a2, " c = ", c2, " g = ", g2];
> (* NestWhile *)
> {{a3, c3, g3}, w} = Rest[NestWhile[regressionStep, {data, {a2, c2,
> g2}, -1.}, regressionTest, 3]];
> Print["FoundParameters: a = ", a3, " c = ", c3, " g = ", g3, " wSOfSq
> = ", w];
> Print["unweightedSOfSq: Start = ", unwSOfSqStart, " Found = ", Plus @@
> ((((#[[2]] - Lorentz[#[[1]], a3, c3, g3])^2) & ) /@ data)];
> Show[{pic,
> Plot[Lorentz[x, a, c, g], {x, xL, xR}, PlotStyle -> RGBColor[0, 0, 1],
> DisplayFunction -> Identity],
> Plot[Lorentz[x, a2, c2, g2], {x, xL, xR}, PlotStyle -> RGBColor[0,
> 1, 0],
> DisplayFunction -> Identity],
> Plot[Lorentz[x, a3, c3, g3], {x, xL, xR}, PlotStyle -> RGBColor[1,
> 0, 0],
> DisplayFunction -> Identity]}, DisplayFunction -> $DisplayFunction]
> ]
>
> From In[80]:=
> StartParameters: a = 18.8227 c = -0.00138881 g = 0.967654
>
> From In[80]:=
> Period catched!
> Penultimate set: a = 21.9753 c = -0.0123743 g = 1.07773 wSOfSq = 849.299
>
> From In[80]:=
> FoundParameters: a = 16.9547 c = 0.0108803 g = 0.898405 wSOfSq = 754.797
>
> From In[80]:=
> unweightedSOfSq: Start = 2854.31 Found = 3078.06
>
> Es ist klar, dass die ungewichtete Summe der quadratischen Abweichungen
> des Startwerts kleiner ist als die des Endwertes - der Startwert ist der
> beste ungewichtete Fit.
>
> Mit den besten Grüssen
> Udo.
>
>
--
Supergünstige DSL-Tarife + WLAN-Router für 0,- EUR*
Jetzt zu GMX wechseln und sparen http://www.gmx.net/de/go/dsl