Hallo Uli,
was genau mit den Gewichten passiert ist in Mathematica bei Dir, kann
man mangels Beispiel nicht nachvollziehen.
Nur zum Sagen, auf eine Idee wie die folgende wirst Du als Physiker
sicher nicht gekommen sein:
In[1]:= <<Statistics`NonlinearFit`
In[2]:= <<Statistics`ContinuousDistributions`
In[8]:= With[{? = 0., ? = 0.75, ? = 0.25},
pic = Plot[1/(x^2 + (?/2)^2) + Random[NormalDistribution[?, ?]], {x, -2,
2},
PlotDivision -> 30, DisplayFunction -> Identity]; data = pic[[1,1,1,1]];
Print["There are ", Length[data], " datapoints."];
NonlinearRegress[data, 1/(x1^2 + (?1/2)^2) + Exp[-((x1 - ?1)^2/(2*?1^2))]/
(Sqrt[2*Pi]*?1), {x1}, {?1, ?1, ?1}, RegressionReport ->
BestFitParameters]]
"There are "860" datapoints."
Out[8]=
{BestFitParameters -> {?1 -> 0.0037703537868760337, ?1 ->
0.07695177010624978,
?1 -> 0.26111019517895107}}
Dieser Ansatz für die Fitfunktion ist falsch, weil (1)
NormalDistribution[my, sigma] positive und negative Werte abgibt, aber
die Normalverteilung im Ansatz bekanntermassen nicht; (2) die
funktionalen Abhängigkeiten verkehrt sind, das Rauschen ist unabhänhig
von x1. Mit dem Ansatz erzeugt man nur systematische Fehler.
Man kann testen, inwiefern die Lorentzdaten von einer Normalverteilung
verrauscht wurden, das ist ein Test von 2 Hypothesen; Fitter testen
immer Hypothesen in Form der Ansätze, in dem Fall habe ich keinen Weg
gefunden, beide Hypthesen mit einem Aufruf von NonlinearRegression zu
testen:
In[161]:= Clear[data, ?1, ?2, gauss];
With[{? = 0., ? = 0.75, ? = 0.25, dy = 0.1},
pic = Plot[1/(x^2 + (?/2)^2) + Random[NormalDistribution[?, ?]], {x, -2,
2},
PlotDivision -> 30, DisplayFunction -> Identity]; data = pic[[1,1,1,1]];
Print["There are ", Length[data], " data points."];
?2 = NonlinearRegress[data, 1/(x1^2 + (?1/2)^2), {x1}, {?1},
RegressionReport -> BestFitParameters][[1,2,1,2]];
Print["LorentzFit: ?2 = ", ?2];
dist = Sort[Last[Transpose[data]] - (1/(#1^2 + (?2/2)^2) & ) /@
First[Transpose[data]]]; val = First[dist] + dy/2.; oo = 0; gauss = {};
For[o = 1, o <= Length[dist], o++, If[dist[[o]] < val, oo++,
AppendTo[gauss, {val - dy/2., oo}]; val += dy; oo = 0; --o]];
gauss = Transpose[{First[Transpose[gauss]], Last[Transpose[gauss]]/
(dy*Length[dist])}];
Print["Gaussian Fit:"];
NonlinearRegress[gauss, Exp[-((x1 - ?1)^2/(2*?1^2))]/(Sqrt[2*Pi]*?1),
{x1}, {?1, ?1},
RegressionReport -> BestFitParameters]
]
From In[161]:= "There are "860" data points."
From In[161]:= "LorentzFit: ?2 = "0.25040192109887016
From In[161]:= "Gaussian Fit:"
Out[163]= {BestFitParameters -> {?1 -> -0.0067614107281022936, ?1 ->
0.7576946340056632}}
Also ganz gut getroffen diesmal. Andere Runs geben auch schon wesentlich
schlechtere Ergebnisse. Die Ansicht
ListPlot[gauss] ist instruktiv in solchen Fällen.
Schönen Sonntag
Udo.
Uli Jentschura wrote:
Lieber DEMUG-Mitglieder,
hier eine Frage betreffs
Statistik, insbesondere NonlinearRegress.
Bevor ich das Problem im Detail schildere, moechte ich mich vergewissern,
nicht von Vornherein einen Denkfehler gemacht zu haben:
Nehmen wir folgendes an:
1. Wir nehmen eine wohldefinierte Kurve, zum Beispiel eine Lorentzkurve
der Form
f(x) = 1/(x^2 + gamma^2/4)
(bis auf Normierung).
2. Wir addieren ein statistisches Rauschen, mit einer Gauss-schen
Verteilung ("Random" mit "NormalDistribution"), so dass die Linie
varrauscht wird, wobei sich die Standardabweichung der addierten
Zufallszahl ergibt zu
sigma(f) = Sqrt[f] + constant * f (*)
(man/frau koennte sich auch andere Abhaengigkeiten sigma(f)
vorstellen.)
3. Wir fitten mit "NonlinearRegress" die entstandene verrauschte Linie,
wobei wir die verrauschten Daten statistisch gewichten, und zwar genau
mit demjenigen Gewicht 1/sigma(f)^2, das wir vorher benutzt haben, um
die "Breite" des statistischen Rauschens festzulegen.
Ueberlegung (richtig oder falsch?):
Dann muesste doch im statistischen Mittel eine nichtlineare Anpassung an
die verrauschte Linie (unter Zuhilfenahme von "NonlinearRegress") mit
einem Modell f(x) = 1/(x^2 + gamma^2/4) unter Gewichtung der f_i mit
Gleichung (*) fuer den Parameter gamma den "Startwert" von gamma ergeben.
Beobachtung:
Die Breite gamma der Lorentzkurve wird nicht richtig gefittet. Vielmehr
ist diese Breite typischerweise kleiner als gamma. Aber vielleicht liegt
ein Denkfehler vor. (Hoffentlich liegt ein Denkfehler vor!)
Vielen Dank an alle Diskussionsteilnehmer im voraus!
Gruss,
Ulrich Jentschura
------------------------------------------------------------------------------
Ulrich Jentschura Theoretische Quantendynamik/ Physik
Albert-Ludwigs-Universitaet Freiburg, Hermann-Herder-Strasse 3, 79104 Freiburg
Tel. +49-761-203 5956 (fax extension 5883) jentschura@XXXXXXX.de
internet homepage: http://tqd1.physik.uni-freiburg.de/~ulj
DMUG-Archiv, http://www.mathematica.ch/archiv.html