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 *)
--------------
DMUG-Archiv, http://www.mathematica.ch/archiv.html