Re: Statistisches

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

Hallo Martin,

Statt eines ungewichteten Fits könnte die Lösung auch ein iterativer Fit 
sein: Man kann FixedPoint[] in die Suche involvieren. Mit konstanten 
Gewichten (Weights -> Automatic) berechnet  man einen Startwert und dann 
geht es weiter in FixedPoint[] mit dem jeweils vorhergehenden Lorentz. 
Da ich nicht weiss, ob ein Fixpunkt vorliegen muss, missbrauche ich die 
Stichprobengrösse anZ als IterationLimit.

Schlecht sind Weights -> ((Abs[# - Lorentz[#, a0, c0, g0]]/#) &), welche 
den von Ihnen beschriebenen Einschnüreffekt zeigen. Recht gut sind dem 
Augenschein nach Weights -> ((Abs[# - Lorentz[#, a0, c0, g0]]/Lorentz[#, 
a0, c0, g0]) & ). Die Änderung kann leicht in regressionStep[] erfolgen.

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}, Weights -> ((Abs[# - 
Lorentz[#, a0, c0, g0]]/Lorentz[#, a0, c0, g0]) & ), RegressionReport -> 
BestFitParameters])}
   ]

In[128]:= (* using FixedPoint *)
With[{a = 20., c = 0., g = 1., xL = -2., xR = 2., anZ = 100, j = 0},
 data = ({#, Lorentz[#, a, c, g] + Random[NormalDistribution[0., 
sigma[Lorentz[#, a, c, g]]]]} & ) /@ Sort[Table[Random[Real, {xL, xR}, 
8], {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]);
 Print["StartParameters:\nBestFitParameters: a = ", a2, " c = ", c2, " 
g = ", g2];
 (* FixedPoint *)
 {a3, c3, g3} = Last[FixedPoint[regressionStep, {data, {a2, c2, g2}}, 
anZ]]; 
 Print["FixedParameters:\nBestFitParameters: a = ", a3, " c = ", c3, " 
g = ", g3];
 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]
]

In Ihrer Diskussion vermisst man ein wenig die Würdigung der 
experimentellen Daten. Wenn bei einer endlichen Stichprobe der Grösse 
anZ die Werte einseitig ausstreuen,  kann man dies auch nicht durch die 
Gewichte korrigieren. Ein korrekter Fitter folgt immer den Daten, im 
Ernstfall führt er den Ansatz ad absurdum. Nur wenn die Daten 
deterministisch und paritätisch streuen - mit anderen Worten, sie 
streuen im Mittel nicht, auch nicht in genügend grossen Teilstichproben 
- kann man den Rückerhalt der Startwerte fordern.

Bei wirklich streuenden Daten gibt es sowohl Verläufe wie:
From In[128]:=
StartParameters:
BestFitParameters: a = 18.0543  c = 0.0145595  g = 0.913402

From In[128]:=
FixedParameters:
BestFitParameters: a = 16.3942  c = 0.00917193  g = 0.849749

als auch nach einigem Wiederholen Verläufe wie diesen:

From In[134]:=
StartParameters:
BestFitParameters: a = 19.4353  c = 0.00802898  g = 0.980085

From In[134]:=
FixedParameters:
BestFitParameters: a = 20.3098  c = 0.00210746  g =  0.958748

Um es jetzt letzten Endes wirklich zu sehen, müsste man \[Chi]^2 
ausgeben. Genaugenommen muss \[Chi]^2 als Fixpunktkriterium verwendet 
werden, was ich Ihnen gern zum Programmieren und Experiemtieren überlasse.

Mit den besten Grüssen
Udo.

Martin Haas wrote:

> Hallo Udo,
>
> ich glaube, wir haben die Lösung. Man kann das Problem 
> herunterkochen auf das Fitten einer Konstanten, und 
> erhält den selben Effekt. Die Konstante wird immer zu 
> klein gefittet (s.u.). Kleine Datenpunkte haben eben ein 
> grösseres relatives Gewicht, wenn man die Gewichte
>
> > auf der Basis der gemessenen Werte <
>
> verteilt. Würde man sie auf der Basis der wahren Werte 
> berechnen, wäre alles in Butter, diese kennt man aber
> nicht, sonst bräuchte man gar nicht zu fitten. Unsere 
> Vermutung vom 02.08. trifft also zu und im Falle der 
> Lorentzfits äussert sich dies in einer kleineren Amplitude 
> und geringeren Breite, weil die Punkte, die weit von der 
> Linienmitte entfernt sind, die Linie sozusagen einschnüren.
>
> Die Lösung ist ein ungewichteter Fit, denn da wird kein 
> Punkt aufgrund seiner Lage (wahr oder real) ausgezeichnet.
> Entfernt man im Code unten die Gewichtung, erhält man 
> Fitergebnisse, die symmetrisch ("richtig") um den Startwert 
> verteilt sind und dies amplitudenunabgängig.
>
> Viele Grüße,
>
> Martin
>
> --------------
>
> << Statistics`ContinuousDistributions`;
> << Statistics`NonlinearFit`;
>
> sigma[n_] := Sqrt[n] + 0.025 n;
>
> Fitdiffs = {};
> amp = 50;
> For[k = 1, k <= 100, k++,
> {
>   thedata = 
>     Table[amp + Random[NormalDistribution[0, sigma[amp]]], {i, 1, 50}];
>   theweights = Table[1/sigma[thedata[[i]]]^2, {i, 1, 50}];
>   constfit = NonlinearRegress[thedata, c, x, c, Weights -> theweights, 
>          RegressionReport -> BestFitParameters];
>   AppendTo[Fitdiffs, amp - c /. (BestFitParameters /. constfit)];
> }
> ];
>
> ListPlot[Fitdiffs];
>
> (* fast alle positiv *)
>
> FitDifferences = {};
>
> For[amp = 50, amp <= 1000, amp *= Exp[Log[1000/50]/10],
> {
>   Fitdiffs = {};
>   For[k = 1, k <= 100, k++,
>    {
>      thedata = Table[amp + Random[NormalDistribution[0, sigma[amp]]],
>                      {i, 1, 50}];
>      theweights = Table[1/sigma[thedata[[i]]]^2, {i, 1, 50}];
>      constfit = NonlinearRegress[thedata, c, x, c, Weights ->
>                    theweights, RegressionReport -> BestFitParameters];
>      AppendTo[Fitdiffs, amp - c /. (BestFitParameters /. constfit)];
>      (* Abweichung des einzelnen Punkts *)
>   }
>  ];
>      
>  afit = NonlinearRegress[Fitdiffs, a, x, a, 
>         RegressionReport -> BestFitParameters];
>  AppendTo[FitDifferences, (a /. (BestFitParameters /. afit))/amp];
>  (* relativer Fehler des Fits *)
> }
> ];
>
> ListPlot[FitDifferences, Frame -> True, Axes -> False];
>
> (* immer positiv, relativer Fehler wird kleiner bei wachsender Amplitude *)
>
> --------------
>