Statistisches

[Udo und Susanne Krause <su.krause@XXXXXXX.ch>]

Hallo nocheinmal,

mir ist noch ein Fehler in der funktionalen Programmierung aufgefallen:

>(siehe (*HIER!*))
>
>> In[108]:= Clear[regressionStep];
>> regressionStep[{data_List, {a0_, c0_, g0_}}] :=
>> Module[{x1},
>>   {data, {a1, c1, g1} /. (BestFitParameters /. NonlinearRegress[data,
>> Lorentz[x1, a1, c1, g1], {x1}, {a1, c1, g1}, (*HIER: *)
>>  Weights -> ((Abs[# -
>> Lorentz[#, a0, c0, g0]]/Lorentz[#, a0, c0, g0]) & ), RegressionReport ->
>> BestFitParameters])}
>>     ]

dadurch bekommt der Lorentz[] den y-Wert oder die Response als Argument. 
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.